Biopython - Big Discrepancy Calculating RNA melting Temperature over Literature

Question

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

Answer 1

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.

Answer 2

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