Package Bio :: Package SeqUtils :: Module MeltingTemp
[hide private]
[frames] | no frames]

Source Code for Module Bio.SeqUtils.MeltingTemp

  1  import math 
2 -def Tm_staluc(s,dnac=50,saltc=50,rna=0):
3 """Returns DNA/DNA tm using nearest neighbor thermodynamics. 4 5 dnac is DNA concentration [nM] 6 saltc is salt concentration [mM]. 7 rna=0 is for DNA/DNA (default), for RNA, rna should be 1. 8 9 Sebastian Bassi <sbassi@genesdigitales.com>""" 10 11 #Credits: 12 #Main author: Sebastian Bassi <sbassi@genesdigitales.com> 13 #Overcount function: Greg Singer <singerg@tcd.ie> 14 #Based on the work of Nicolas Le Novere <lenov@ebi.ac.uk> Bioinformatics. 15 #17:1226-1227(2001) 16 17 #This function returns better results than EMBOSS DAN because it uses 18 #updated thermodynamics values and takes into account inicialization 19 #parameters from the work of SantaLucia (1998). 20 21 #Things to do: 22 #+Detect complementary sequences. Change K according to result. 23 #+Add support for heteroduplex (see Sugimoto et al. 1995). 24 #+Correction for Mg2+. Now supports only monovalent ions. 25 #+Put thermodinamics table in a external file for users to change at will 26 #+Add support for danglings ends (see Le Novele. 2001) and mismatches. 27 28 dh=0 #DeltaH. Enthalpy 29 ds=0 #deltaS Entropy 30 31 def tercorr(stri): 32 deltah=0 33 deltas=0 34 if rna==0: 35 #DNA/DNA 36 #Allawi and SantaLucia (1997). Biochemistry 36 : 10581-10594 37 if stri.startswith('G') or stri.startswith('C'): 38 deltah=deltah-0.1 39 deltas=deltas+2.8 40 elif stri.startswith('A') or stri.startswith('T'): 41 deltah=deltah-2.3 42 deltas=deltas-4.1 43 if stri.endswith('G') or stri.endswith('C'): 44 deltah=deltah-0.1 45 deltas=deltas+2.8 46 elif stri.endswith('A') or stri.endswith('T'): 47 deltah=deltah-2.3 48 deltas=deltas-4.1 49 dhL=dh+deltah 50 dsL=ds+deltas 51 return dsL,dhL 52 elif rna==1: 53 #RNA 54 if stri.startswith('G') or stri.startswith('C'): 55 deltah=deltah-3.61 56 deltas=deltas-1.5 57 elif stri.startswith('A') or stri.startswith('T') or \ 58 stri.startswith('U'): 59 deltah=deltah-3.72 60 deltas=deltas+10.5 61 if stri.endswith('G') or stri.endswith('C'): 62 deltah=deltah-3.61 63 deltas=deltas-1.5 64 elif stri.endswith('A') or stri.endswith('T') or \ 65 stri.endswith('U'): 66 deltah=deltah-3.72 67 deltas=deltas+10.5 68 dhL=dh+deltah 69 dsL=ds+deltas 70 # print "delta h=",dhL 71 return dsL,dhL
72 73 def overcount(st,p): 74 """Returns how many p are on st, works even for overlapping""" 75 ocu=0 76 x=0 77 while 1: 78 try: 79 i=st.index(p,x) 80 except ValueError: 81 break 82 ocu=ocu+1 83 x=i+1 84 return ocu 85 86 R=1.987 # universal gas constant in Cal/degrees C*Mol 87 sup=s.upper() 88 vsTC,vh=tercorr(sup) 89 vs=vsTC 90 91 k=(dnac/4.0)*1e-9 92 #With complementary check on, the 4.0 should be changed to a variable. 93 94 if rna==0: 95 #DNA/DNA 96 #Allawi and SantaLucia (1997). Biochemistry 36 : 10581-10594 97 vh=vh+(overcount(sup,"AA"))*7.9+(overcount(sup,"TT"))*7.9+(overcount(sup,"AT"))*7.2+(overcount(sup,"TA"))*7.2+(overcount(sup,"CA"))*8.5+(overcount(sup,"TG"))*8.5+(overcount(sup,"GT"))*8.4+(overcount(sup,"AC"))*8.4 98 vh=vh+(overcount(sup,"CT"))*7.8+(overcount(sup,"AG"))*7.8+(overcount(sup,"GA"))*8.2+(overcount(sup,"TC"))*8.2 99 vh=vh+(overcount(sup,"CG"))*10.6+(overcount(sup,"GC"))*9.8+(overcount(sup,"GG"))*8+(overcount(sup,"CC"))*8 100 vs=vs+(overcount(sup,"AA"))*22.2+(overcount(sup,"TT"))*22.2+(overcount(sup,"AT"))*20.4+(overcount(sup,"TA"))*21.3 101 vs=vs+(overcount(sup,"CA"))*22.7+(overcount(sup,"TG"))*22.7+(overcount(sup,"GT"))*22.4+(overcount(sup,"AC"))*22.4 102 vs=vs+(overcount(sup,"CT"))*21.0+(overcount(sup,"AG"))*21.0+(overcount(sup,"GA"))*22.2+(overcount(sup,"TC"))*22.2 103 vs=vs+(overcount(sup,"CG"))*27.2+(overcount(sup,"GC"))*24.4+(overcount(sup,"GG"))*19.9+(overcount(sup,"CC"))*19.9 104 ds=vs 105 dh=vh 106 107 else: 108 #RNA/RNA hybridisation of Xia et al (1998) 109 #Biochemistry 37: 14719-14735 110 vh=vh+(overcount(sup,"AA"))*6.82+(overcount(sup,"TT"))*6.6+(overcount(sup,"AT"))*9.38+(overcount(sup,"TA"))*7.69+(overcount(sup,"CA"))*10.44+(overcount(sup,"TG"))*10.5+(overcount(sup,"GT"))*11.4+(overcount(sup,"AC"))*10.2 111 vh=vh+(overcount(sup,"CT"))*10.48+(overcount(sup,"AG"))*7.6+(overcount(sup,"GA"))*12.44+(overcount(sup,"TC"))*13.3 112 vh=vh+(overcount(sup,"CG"))*10.64+(overcount(sup,"GC"))*14.88+(overcount(sup,"GG"))*13.39+(overcount(sup,"CC"))*12.2 113 vs=vs+(overcount(sup,"AA"))*19.0+(overcount(sup,"TT"))*18.4+(overcount(sup,"AT"))*26.7+(overcount(sup,"TA"))*20.5 114 vs=vs+(overcount(sup,"CA"))*26.9+(overcount(sup,"TG"))*27.8+(overcount(sup,"GT"))*29.5+(overcount(sup,"AC"))*26.2 115 vs=vs+(overcount(sup,"CT"))*27.1+(overcount(sup,"AG"))*19.2+(overcount(sup,"GA"))*32.5+(overcount(sup,"TC"))*35.5 116 vs=vs+(overcount(sup,"CG"))*26.7+(overcount(sup,"GC"))*36.9+(overcount(sup,"GG"))*32.7+(overcount(sup,"CC"))*29.7 117 ds=vs 118 dh=vh 119 120 ds=ds-0.368*(len(s)-1)*math.log(saltc/1e3) 121 tm=((1000* (-dh))/(-ds+(R * (math.log(k)))))-273.15 122 # print "ds="+str(ds) 123 # print "dh="+str(dh) 124 return tm 125