Virtual electrophotogram

CairoNumberWidth.py

#!/home/zhangcg/local/python/bin/python
# -*- coding:utf-8 -*-

#-----------------------------------------------------
# File: CairoNumberWidth.py
# Date: 2008-07-04
# Description: 
#-----------------------------------------------------

__version__ = '1.0'
__author__  = 'Wubin Qu <quwubin@gmail.com> @ZCGLAB @BMI @CHINA'
__blog__    = 'http://quwubin.cnblogs.com'
__license__ = 'GPL'


font50 = {
"a" : 23.0,
"b" : 26.0,
"c" : 22.0,
"d" : 25.0,
"e" : 22.0,
"f" : 17.0,
"g" : 24.0,
"h" : 26.0,
"i" : 14.0,
"j" : 14.0,
"k" : 25.0,
"l" : 14.0,
"m" : 39.0,
"n" : 25.0,
"o" : 24.0,
"p" : 25.0,
"q" : 25.0,
"r" : 17.0,
"s" : 20.0,
"t" : 15.0,
"u" : 25.0,
"v" : 25.0,
"w" : 36.0,
"x" : 25.0,
"y" : 25.0,
"z" : 22.0,
"0" : 25.0,
"1" : 26.0,
"2" : 25.0,
"3" : 25.0,
"4" : 25.0,
"5" : 25.0,
"6" : 25.0,
"7" : 25.0,
"8" : 26.0,
"9" : 25.0,
"A" : 36.0,
"B" : 34.0,
"C" : 33.0,
"D" : 36.0,
"E" : 30.0,
"F" : 28.0,
"G" : 36.0,
"H" : 35.0,
"I" : 16.0,
"J" : 19.0,
"K" : 36.0,
"L" : 30.0,
"M" : 44.0,
"N" : 36.0,
"O" : 37.0,
"P" : 28.0,
"Q" : 37.0,
"R" : 33.0,
"S" : 28.0,
"T" : 30.0,
"U" : 36.0,
"V" : 36.0,
"W" : 47.0,
"X" : 36.0,
"Y" : 36.0,
"Z" : 31.0,
"_" : 25.0,
"-" : 16.0,
"|" : 9.0,
"!" : 16.0,
"@" : 47.0,
"#" : 27.0,
"$" : 25.0,
"%" : 41.0,
"^" : 23.0,
"&" : 39.0,
"*" : 26.0,
"(" : 16.0,
")" : 17.0,
"+" : 28.0,
"[" : 17.0,
"]" : 16.0,
"{" : 24.0,
"}" : 24.0,
":" : 14.0,
'"' : 21.0,
";" : 14.0,
"'" : 9.0,
"\\" : 14.0,
"," : 13.0,
"." : 12.0,
"/" : 14.0,
"?" : 22.0,
">" : 28.0,
"<" : 28.0,
}


font30 = {
" " : 14.0,
"a" : 14.0,
"b" : 15.0,
"c" : 13.0,
"d" : 15.0,
"e" : 14.0,
"f" : 10.0,
"g" : 15.0,
"h" : 15.0,
"i" : 8.0,
"j" : 9.0,
"k" : 16.0,
"l" : 9.0,
"m" : 24.0,
"n" : 15.0,
"o" : 15.0,
"p" : 15.0,
"q" : 15.0,
"r" : 10.0,
"s" : 12.0,
"t" : 8.0,
"u" : 15.0,
"v" : 15.0,
"w" : 22.0,
"x" : 15.0,
"y" : 15.0,
"z" : 13.0,
"0" : 15.0,
"1" : 15.0,
"2" : 15.0,
"3" : 15.0,
"4" : 16.0,
"5" : 14.0,
"6" : 15.0,
"7" : 15.0,
"8" : 15.0,
"9" : 15.0,
"A" : 22.0,
"B" : 19.0,
"C" : 20.0,
"D" : 21.0,
"E" : 19.0,
"F" : 17.0,
"G" : 21.0,
"H" : 21.0,
"I" : 9.0,
"J" : 12.0,
"K" : 22.0,
"L" : 18.0,
"M" : 26.0,
"N" : 21.0,
"O" : 21.0,
"P" : 17.0,
"Q" : 21.0,
"R" : 20.0,
"S" : 17.0,
"T" : 19.0,
"U" : 21.0,
"V" : 22.0,
"W" : 28.0,
"X" : 22.0,
"Y" : 21.0,
"Z" : 19.0,
"_" : 15.0,
"-" : 9.0,
"|" : 6.0,
"!" : 10.0,
"@" : 27.0,
"#" : 17.0,
"$" : 15.0,
"%" : 25.0,
"^" : 14.0,
"&" : 23.0,
"*" : 15.0,
"(" : 10.0,
")" : 10.0,
"+" : 17.0,
"[" : 11.0,
"]" : 10.0,
"{" : 15.0,
"}" : 15.0,
":" : 8.0,
'"' : 12.0,
";" : 8.0,
"'" : 5.0,
"\\" : 8.0,
"," : 8.0,
"." : 7.0,
"/" : 8.0,
"?" : 14.0,
">" : 17.0,
"<" : 17.0,
}


font20 = {
" " : 9.0,
"a" : 9.0,
"b" : 10.0,
"c" : 9.0,
"d" : 11.0,
"e" : 9.0,
"f" : 7.0,
"g" : 11.0,
"h" : 9.0,
"i" : 6.0,
"j" : 6.0,
"k" : 10.0,
"l" : 6.0,
"m" : 16.0,
"n" : 11.0,
"o" : 10.0,
"p" : 11.0,
"q" : 10.0,
"r" : 7.0,
"s" : 8.0,
"t" : 5.0,
"u" : 10.0,
"v" : 10.0,
"w" : 14.0,
"x" : 10.0,
"y" : 10.0,
"z" : 9.0,
"0" : 11.0,
"1" : 10.0,
"2" : 10.0,
"3" : 10.0,
"4" : 11.0,
"5" : 10.0,
"6" : 10.0,
"7" : 10.0,
"8" : 10.0,
"9" : 11.0,
"A" : 14.0,
"B" : 12.0,
"C" : 14.0,
"D" : 15.0,
"E" : 12.0,
"F" : 12.0,
"G" : 15.0,
"H" : 14.0,
"I" : 6.0,
"J" : 8.0,
"K" : 14.0,
"L" : 12.0,
"M" : 17.0,
"N" : 14.0,
"O" : 14.0,
"P" : 11.0,
"Q" : 14.0,
"R" : 13.0,
"S" : 10.0,
"T" : 12.0,
"U" : 14.0,
"V" : 14.0,
"W" : 19.0,
"X" : 14.0,
"Y" : 14.0,
"Z" : 13.0,
"_" : 10.0,
"-" : 7.0,
"|" : 4.0,
"!" : 7.0,
"@" : 18.0,
"#" : 12.0,
"$" : 10.0,
"%" : 16.0,
"^" : 9.0,
"&" : 16.0,
"*" : 10.0,
"(" : 7.0,
")" : 7.0,
"+" : 12.0,
"[" : 7.0,
"]" : 7.0,
"{" : 10.0,
"}" : 10.0,
":" : 6.0,
'"' : 8.0,
";" : 6.0,
"'" : 4.0,
"\\" : 6.0,
"," : 5.0,
"." : 4.0,
"/" : 6.0,
"?" : 9.0,
">" : 12.0,
"<" : 12.0,
}

GelMobility.py

#!/home/zhangcg/local/python/bin/python
# -*- coding: utf-8 -*- #   
'''Mobility of the PCR amplicons in agarose gel electrophoresis
was calculated following the semi-log formula reported by Helling
et al [Ref. 1] in 1974.

    Y = a - b * ln (X + k)

  Where Y is the mobility and X is the size of the PCR amplicons.

  To get an reliable prediction of mobility of PCR amplicons, we
determined these parameters (a, b and k) by statistically analyzing
the commercially available DNA markers of known size (100 bp, 250 bp,
500 bp, 750 bp, 1000 bp and 2000 bp, from TaKaRa Inc.) in agarose gel
with different concentrations. The DNA fragment of 100 bp was taken as
a reference band in the agarose gel electrophoresis, and the relative
mobility was determined according to the 100 bp DNA band. When the 
mobility of the reference band reached 5 cm in running the agarose gel,
the electrophoresis was stoped and the mobility of other DNA bands
were measure to determine the relative mobility. Briefly, the 
parameters in the abovementioned formula were determined on the 
0.5% ~ 2% agarose gel electrophoresis with triplicate measurements
in out lab and shown below for producing the virtual electrophoretogram
in MPprimer program.

For 0.5% agarose gel concentration:

    a = 2.7094
    b = 0.2691
    k = 464.4412
    R^2 = 0.9989

For 1.0% agarose gel concentration:

    a = 2.3977
    b = 0.2700
    k = 73.9788
    R^2 = 0.9984

For 1.5% agarose gel concentration:

    a = 2.3221
    b = 0.2634
    k = 48.0873
    R^2 = 0.9977

For 2.0% agarose gel concentration:

    a = 2.1333
    b = 0.2561
    k = 18.5417
    R^2 = 0.9948

Reference

[1] Helling, R. B., Goodman, H. M., & Boyer, H. W. (1974). Analysis of 
endonuclease R-EcoRI fragments of DNA from lambdoid bacteriophages and 
other viruses by agarose-gel electrophoresis. Journal of Virology, 
14(5): 1235-1244.
'''

# Author: Wubin Qu [quwubin@gmail.com], BIRM, China
# Date: 2009-10-12
# License: GPL v3

import sys

def load_gel_para_dict(gel_conc=1.0, formula='Helling'):
    gel_para_dict = {}
    gel_para_dict['Helling'] = {}

    # Parameter for Helling's formula
    gel_para = {}
    gel_para[0.5] = {
	'a' : 2.7094,
	'b' : 0.2691,
	'k' : 464.4412,
    }

    gel_para[1.0] = {
	'a' : 2.3977,
	'b' : 0.2700,
	'k' : 73.9788,
    }

    gel_para[1.5] = {
	'a' : 2.3221,
	'b' : 0.2634,
	'k' : 48.0873,
    }

    gel_para[2.0] = {
	'a' : 2.1333,
	'b' : 0.2561,
	'k' : 18.5417,
    }

    gel_para_dict['Helling'] = gel_para

    err_msg = 'Gel concentration or formula is illegal, Currently, only 0.5, 1, 1.5, 2.0 for concentration and "Helling" for formula are allowed.'

    try: 
        gel_conc = float(gel_conc)
        a = gel_para_dict[formula][gel_conc]['a']
        b = gel_para_dict[formula][gel_conc]['b']
        k = gel_para_dict[formula][gel_conc]['k']
    except:
	print >> sys.stderr, err_msg
	exit()

    return gel_para_dict, a, b, k

def get_size_range(size, gel_conc=1.0, ref_mobility=50, offset=2, formula='Helling'):
    Y = cal_mobility(size, gel_conc)
    offset = int(offset)
    # Set 2 mm as the distance which the bands can be \
            #seperated by naked eyes
    Xmin = cal_size(Y + offset, gel_conc=gel_conc, ref_mobility=ref_mobility, formula='Helling')
    Xmax = cal_size(Y - offset, gel_conc=gel_conc, ref_mobility=ref_mobility, formula='Helling')

    return Xmin, Xmax

def cal_mobility(X, gel_conc=1.0, ref_mobility=50, formula='Helling'):
    '''Cal mobility based on size'''
    import math
    gel_para_dict, a, b, k = load_gel_para_dict(gel_conc=gel_conc, formula=formula)

    X = float(X)
    gel_conc = float(gel_conc)

    # X: size (bp)
    # ref_mobility: the mobility distance of the fastest DNA segment

    if formula == 'Helling':
        Y = a - b * math.log(X + k)
    else:
        pass
        #Y = math.exp(a - b * math.log(X + k))

    # Y: the relative mobility = mobility distance / ref_mobility
    Y = Y * ref_mobility
    # Y: the mobility distance
    return round(Y, 1)

def cal_size(Y, gel_conc=1.0, ref_mobility=50, formula='Helling'):
    '''Predict size based on the relative mobility'''
    import math

    gel_para_dict, a, b, k = load_gel_para_dict(gel_conc=gel_conc, formula=formula)

    # Y: the mobility distance
    Y = Y / ref_mobility
    # ref_mobility: the mobility distance of the fastest DNA segment
    if formula == 'Helling':
        #Y = a - b * math.log(X + k)
        X = math.exp((a - Y) / b) - k
    else:
        pass

    return int(round(X, 0))

def main ():
    '''Test'''
    import sys
    #X = 100
    #gel_conc = sys.argv[1]
    #mobility = cal_mobility(X, gel_conc, ref_mobility=50, formula='Helling')
    #print mobility
    #mobility = cal_mobility(X, gel_conc, ref_mobility=50)
    #print mobility
    min, max = get_size_range(sys.argv[2], gel_conc=1.0, ref_mobility=50, offset=sys.argv[1], formula='Helling')
    print min, max


if __name__ == '__main__':
    main()

ve.py

#!/home/zhangcg/local/python/bin/python
# -*- coding:utf-8 -*-
from __future__ import division

__version__ = '1.0'
__author__  = 'Wubin Qu <quwubin@gmail.com>'
__license__ = 'GPL v3'

import os
import time
import CairoNumberWidth as CNW
import cairo
import GelMobility

def draw_title(ctx, height, width):
    '''Draw title'''
    title = 'Virtual Electrophotogram'

    title_length = 0
    for i in range(len(title)):
        title_length = title_length + CNW.font30[title[i]]

    ctx.move_to(width/2 - title_length/2, 50)
    ctx.set_source_rgb(0, 0, 237)
    ctx.select_font_face("Times", cairo.FONT_SLANT_NORMAL)
    ctx.set_font_size(30)
    ctx.show_text(title)

    ISOTIMEFORMAT = '%Y-%m-%d %X'
    title = 'Created by VE @ %s' % time.strftime(ISOTIMEFORMAT, time.localtime(time.time()))

    title_length = 0
    for i in range(len(title)):
        title_length = title_length + CNW.font20[title[i]]

    ctx.move_to(width/2 - title_length/2, 80)
    ctx.set_source_rgb(0, 0, 237)
    ctx.select_font_face("Times", cairo.FONT_SLANT_NORMAL)
    ctx.set_font_size(20)
    ctx.show_text(title)

    ctx.stroke()

    return ctx

def draw_foot(ctx, height, width, agarose, product_number):
    '''Draw title'''
    title = 'M: Marker DL2000, %s agarose gel' % agarose

    title_length = 0
    for i in range(len(title)):
        title_length = title_length + CNW.font20[title[i]]

    ctx.move_to(width/2 - title_length/2, 715)
    ctx.set_source_rgb(1, 1, 1)
    ctx.select_font_face("Times", cairo.FONT_SLANT_NORMAL)
    ctx.set_font_size(20)
    ctx.show_text(title)

    ctx.stroke()

    title = '1-%s: PCR amplicons' % (product_number + 1)
    title_length = 0
    for i in range(len(title)):
        title_length = title_length + CNW.font20[title[i]]

    ctx.move_to(width/2 - title_length/2, 738)
    ctx.set_source_rgb(1, 1, 1)
    ctx.select_font_face("Times", cairo.FONT_SLANT_NORMAL)
    ctx.set_font_size(20)
    ctx.show_text(title)

    ctx.stroke()

    return ctx

def draw_virtual_elec(ctx, height, width, idname_size, agarose, product_number):
    '''Draw'''
    # Draw slot
    ctx.select_font_face("Times", cairo.FONT_SLANT_NORMAL)
    x = 30 + 80
    y = 160
    for i in range(product_number + 2):
	ctx.set_source_rgba(0, 0, 0, 1.0)
	ctx.rectangle(x, y, 70, 500)
	ctx.fill()
	x = x + 90

    x = 30 + 80
    y = 160
    for i in range(product_number + 2):

	if i == 0:
	    line_title = 'M'
	else:
	    line_title = str(i)

	line_title_length = 0
	for j in range(len(line_title)):
	    line_title_length = line_title_length + CNW.font20[line_title[j]]

	ctx.move_to(x + 35 - line_title_length/2, y-10)
	ctx.set_source_rgb(1, 1, 1)
	ctx.set_font_size(20)
	ctx.show_text(line_title)

	if i != 0 and i <= product_number:
	    id, size = idname_size[i-1]
	    y1 = GelMobility.cal_mobility(size, gel_conc=agarose, ref_mobility=50, formula='Helling')
	    y1 = y1*13 - 40# 13 is just a amplificationfactor
	    ctx.move_to(x+10, y1)
	    ctx.set_line_width(6)
	    ctx.line_to(x + 50, y1)

	    x2 = 30 + 80 + 90 * (product_number + 1) + 10
	    ctx.move_to(x2, y1)
	    ctx.set_line_width(6)
	    ctx.line_to(x2 + 50, y1)
	    ctx.stroke()

	    x3 = 30 + 80 + 90 * (product_number + 2) 
	    ctx.move_to(x3, y1+4)
	    bp_num = str(size).rjust(4) + ' bp'
	    #bp_num = size + ' bp'
	    ctx.show_text(bp_num)
	    ctx.stroke()

	    y4 = y + 525
	    id_length = 0
	    for i in range(len(id)):
		id_length = id_length + CNW.font20[id[i]]
	    ctx.move_to(x + 35 - id_length/2, y4)
	    ctx.set_source_rgb(1, 1, 1)
	    ctx.set_font_size(20)
	    ctx.show_text(id)


	x = x + 90

    
    # Mark 2000 kb
    mark_size_list = [2000, 1000, 750, 500, 250, 100]
    mark_size_dict = {
	2000 : '2000 bp',
	1000 : '1000 bp',
	750  : ' 750 bp',
	500  : ' 500 bp',
	250  : ' 250 bp',
	100  : ' 100 bp',
    }
    
    ctx.set_source_rgb(1, 1, 1)
    for mark_size in mark_size_list:
	x = 30 + 80 + 10
	y = GelMobility.cal_mobility(mark_size, gel_conc=agarose, ref_mobility=50, formula='Helling')
	y = y*13 - 40# 13 is just a amplificationfactor

	ctx.move_to(x, y)
        #x, y = ctx.get_current_point()
        ctx.set_line_width(6)
        ctx.line_to(x + 50, y)
        ctx.stroke()

        ctx.move_to(x - 90, y+4)
	ctx.show_text(mark_size_dict[mark_size])
        ctx.stroke()	
	
    return ctx

def paint(size_list, agarose, pic_path):
    '''Paint the virtual electrophotogram '''
    size_list = list(set(size_list))
    size_list.sort()
    idname_size = [(str(size), size) for size in size_list]

    height = 750
    # Marker line width: 80
    product_number = len(idname_size)
    width = (product_number + 2) * 90 + 80 + 60 + 80

    # Set cairo
    surface = cairo.ImageSurface(cairo.FORMAT_ARGB32, width, height)
    ctx = cairo.Context(surface)

    # Mask
    ctx.set_source_rgba(0, 0, 0, 0.9)
    ctx.rectangle(0, 100, width, height)
    ctx.fill()

    ctx = draw_title(ctx, height, width)
    ctx = draw_virtual_elec(ctx, height, width, idname_size, agarose, product_number)
    ctx = draw_foot(ctx, height, width, agarose, product_number)

    #output a PNG file
    surface.write_to_png(pic_path)
    
    return pic_path 


def main ():
    size_list = [100, 200, 300, 400, 500]
    agarose = 1
    pic_path = 've.png'
    ele_picture = paint(size_list, agarose, pic_path)


if __name__ == '__main__':
    main()

Download these three script and put them into a directory.

Run ./ve.py to get an example virtual electrophotogram called "ve.png".