12
$\begingroup$

I experience big discrepancies when calculating melting temperature of RNA 7-mers with Biopython over values generated by a popular algorithm.

I tried the nearest neighbour algorithm with RNA and salt concentrations as described in a respective paper (thermodynamic table used as in paper below from: Freier et al 1986). Yet, the values largely differ (execute code below to see). I tried all seven salt correction methods provided by Biopython, still I never get close to the values generated by siRNA design algorithm for the same 7-mers.

Can someone tell me how accurate Biopython's melting temperature nearest neighbour algorithm is? Especially for short oligomers like my 7-mers? Is there maybe something I am implementing wrong? Any suggestions?

Values derived from executing sample input: http://sidirect2.rnai.jp/ Tm is given for the seed duplex of the guide strand: bases 2-7

Literature: "Thermodynamic stability and Watson–Crick base pairing in the seed duplex are major determinants of the efficiency of the siRNA-based off-target effect" http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2602766/pdf/gkn902.pdf

from Bio.Seq import Seq  
from Bio.SeqUtils import MeltingTemp

test_list = [
('GGAUUUG', 21.5),
('CUCAUUG', 18.1),
('CAUAUUC', 8.7),
('UUUGAGU', 19.2),
('UUUUGAG', 12.2),
('GUUUCAA', 14.9),
('AGUUUCG', 19.7),
('GAAGUUU', 13.3)
]


for t in test_list:
    myseq = Seq(t[0])
    tm = MeltingTemp.Tm_NN(myseq, dnac1=100, Na=100,   nn_table=MeltingTemp.RNA_NN1, saltcorr=7)  # NN1 = Freier et al (1986)
    tm = round(tm, 1)  # round to one decimal
    print 'BioPython Tm: ' + str(tm) + '  siDirect Tm: ' + str(t[1])


# Output:
>>> BioPython Tm: 10.1  siDirect Tm: 21.5  # GGAUUUG
>>> BioPython Tm: 8.2  siDirect Tm: 18.1   # CUCAUUG
>>> BioPython Tm: -2.0  siDirect Tm: 8.7   # CAUAUUC
>>> BioPython Tm: 8.2  siDirect Tm: 19.2   # UUUGAGU
>>> BioPython Tm: 1.3  siDirect Tm: 12.2   # UUUUGAG
>>> BioPython Tm: 4.4  siDirect Tm: 14.9   # GUUUCAA
>>> BioPython Tm: 9.0  siDirect Tm: 19.7   # AGUUUCG
>>> BioPython Tm: 2.1  siDirect Tm: 13.3   # GAAGUUU

Question first asked at stack-exchange: https://stackoverflow.com/questions/30037939/biopython-big-discrepancy-calculating-rna-melting-temperature-over-literature

$\endgroup$
  • $\begingroup$ Also asked at Biostars: biostars.org/p/140873 $\endgroup$ – pykong May 5 '15 at 18:22
  • $\begingroup$ What version is your Bio? $\endgroup$ – dustin May 5 '15 at 19:19
  • 1
    $\begingroup$ The source code for that function is available on Github, there are also some comments in there explaining what methods Biopython uses. The issue might also be that you're comparing apples and oranges, and that your siRNA Tm are not for 7mer parts of an siRNA but something else. $\endgroup$ – Mad Scientist May 5 '15 at 19:28
  • $\begingroup$ @ dustin: My Biopython is up to date V 1.65. @ WYSIWYG: added outputs for saltcorrection method 7. @ Mad Scientist: Thanks for suggesting the doc, yet I feel drawing conclusions from the included algorithms on my issue is beyond me. I need to think about whether I maybe misunderstood something fundamental, but this I deem quite unlikely. $\endgroup$ – pykong May 6 '15 at 18:11
  • $\begingroup$ I was forwarded a PERL script by the researchers behind siDIRECT 2.0. I still do not know why the big discrepancy as both the physical chemistry as well as PERL is beyond me at this point. But at least I was able to integrate that script into the larger frame of my Python program doing melting point calculations. For this concludes the topic without providing an actual answer. Thanks for your contributions. $\endgroup$ – pykong May 24 '15 at 20:25
4
$\begingroup$

Comparing Biopython MetlingTemp to other calculators.

I have written the recent version of MeltingTemp in Biopython's SeqUtils. I have extensively tested the Tm calculations against other programs like MELTING and Primer3Plus and other online Tm calculators with consistent results, thus I'm pretty confident that there is no gross error in the module. The simple answer in this case is: the calculation of siDirect is wrong.

Sources of discrepancy.

One minor thing: Many programs calculate k as k = total Oligo/4. MeltingTemp uses k = Oligo1 - (Oligo2/2). Thus to mimic this behaviour, you have to use

dnac1=50, dnac2=50

instead of

dnac1=100

And now the major thing: Oligo concentration is usually given in nanomolar. However, if you look at the calculation in the Perl script, siDirect seems to use micromolar!

my \$tm = ( 1000 * $dH / ( -10.8 + $dS + 1.987 * log(0.0001/4) ) - 273.15 + 16.6 * log(0.1)/log(10) )

The critical part here is:

log(0.0001/4)

which should be:

log(0.0000001/4)

Thus, you have to pass 1000-fold higher oligo concentrations to Biopython's MeltingTemp to get the same result as siDirect:

>>> from Bio.SeqUtils import MeltingTemp as mt
>>> print mt.Tm_NN('GGAUUUG', dnac1=50000, dnac2=50000, Na=100,
                  nn_table=mt.RNA_NN1, saltcorr=1)
20.1472140567
>>> print mt.Tm_NN('CUCAUUG', dnac1=50000, dnac2=50000, Na=100,
                   nn_table=mt.RNA_NN1, saltcorr=1)
18.1074422939

This is of course wrong if your primer concentration is 100 nM, but that's the result you get from siDirect.

| improve this answer | |
$\endgroup$
  • $\begingroup$ A very insightful and thorough answer, thanks. I can't help but wonder why siDirect went for micro molar. $\endgroup$ – James Dec 1 '15 at 2:05
  • $\begingroup$ Thank You for your in-depth answer Markus. Indeed striking that siDirect - one of the most used oligo tool on the web could be so fundamentally wrong. This saves SeqUtils's honour! $\endgroup$ – pykong Dec 1 '15 at 21:14
1
$\begingroup$

I just received a PERL code snippet by the researcher behind sidirect2 which allows to calculate the Tm. I still do not know the source of the calculation discrepancy as the thermodynamics is beyond me.

I post the code to conclude the question here. The can below can be saved into a .pl file and run from shell with the respective RNA sequence as a string argument.

Code:

#!/usr/bin/perl


# ====================
sub tm_RNA {

# Calculating Tm for RNA-RNA hybrid using nearest neighbour method
#
# usage: $tm = tm_RNA('GGCUGCCAAGAACCUGCAGG') ;

my $seq = lc ($_[0] // '') ;
$seq =~ /^[augc]+$/ or return '' ;

my $dH  = deltaH_RNA($seq) or return '' ;
my $dS  = deltaS_RNA($seq) or return '' ;

my $tm = ( 1000 * $dH / ( -10.8 + $dS + 1.987 * log(0.0001/4) ) - 273.15 + 16.6 * log(0.1)/log(10) ) ;

return sprintf("%.2f", $tm) ;
} ;
    # ====================
    sub deltaH_RNA {  # delta H parameters
    my $seq = lc ($_[0] // '') ;
    $seq =~ /^[augc]+$/ or return '' ;

my %rna_dH = (
    'aa' =>  -6.6,
    'uu' =>  -6.6,
    'au' =>  -5.7,
    'ua' =>  -8.1,
    'ca' => -10.5,
    'ug' => -10.5,
    'cu' =>  -7.6,
    'ag' =>  -7.6,
    'ga' => -13.3,
    'uc' => -13.3,
    'gu' => -10.2,
    'ac' => -10.2,
    'cg' =>  -8.0,
    'gc' => -14.2,
    'gg' => -12.2,
    'cc' => -12.2
) ;

my $sum_dH = 0 ;
    foreach (0..length($seq) - 2){
    my $dinucleotide = substr($seq, $_, 2) ;
	    $sum_dH += $rna_dH{$dinucleotide} ;
}

return $sum_dH ;
    } ;
    # ====================
    sub deltaS_RNA {  # delta S parameters
    my $seq = lc ($_[0] // '') ;
    $seq =~ /^[augc]+$/ or return '' ;

my %rna_dS = (
    'aa' => -18.4,
    'uu' => -18.4,
    'au' => -15.5,
    'ua' => -22.6,
    'ca' => -27.8,
    'ug' => -27.8,
    'cu' => -19.2,
    'ag' => -19.2,
    'ga' => -35.5,
    'uc' => -35.5,
    'gu' => -26.2,
    'ac' => -26.2,
    'cg' => -19.4,
    'gc' => -34.9,
    'gg' => -29.7,
    'cc' => -29.7
) ;

my $sum_dS = 0 ;
    foreach (0..length($seq) - 2){
    my $dinucleotide = substr($seq, $_, 2) ;
	    $sum_dS += $rna_dS{$dinucleotide} ;
}

return $sum_dS ;
} ;
# ====================
# Calling main method:

my $seq = shift;
    $tm = tm_RNA($seq) ;
    print $tm  # returns value to stdout
| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.