Tm based on Nearest Neighbor, monovalent ion, Mg2+, dNTP DNA

#advanced Tm calculation
#nearest-neighbour thermedynamics calculation of melting temperature
#corrected for monovalent ion, Mg+2, dNTP and even DMSO concentrations
#Credits:
#Main author: Ogan ABAAN <[email protected]>
#Overcount function: Greg Singer <[email protected]>

from math import log

def olcount(s, pattern):
    """Returns how many pattern on s, works for overlapping"""
    count= 0
    x= 0
    while 1:
        try:
            i= s.index(pattern,x)
        except ValueError:
            break
        count+= 1
        x= i+1
    return count



def Tm(s, dna=250, Na=50, K=0, Tris=0, Mg=1.5, dNTP=0.2, DMSO=0):

    """
    ### THE PCR BUFFER/REACTION CONDITIONS FOR CALCULATING ION CONCETRATIONS ###
    # dna      the oligo concentration (nM)
    # Na       the Na+ concentration (mM)
    # K        the K+ concetration (mM)
    # Tris     the Tris+ concentration (mM)
    # Mg       the Mg+2 concentration (mM)
    # dNTP     the dNTP concentration (mM)
    # DMSO     correction for DMSO (%)
    """

    # NN parameters in 1M NaCl (SantaLucia J, PNAS, 1998)
    # deltaH (1nd column unit in kcal/mol)
    # deltaS (2rd column unit in cal/mol K)
    deltaHS= {'AA': [-7.9, -22.2],\
                  'TT': [-7.9, -22.2],\
                  'AT': [-7.2, -20.4],\
                  'TA': [-7.2, -21.3], \
                  'CA': [-8.5, -22.7], \
                  'TG': [-8.5, -22.7], \
                  'GT': [-8.4, -22.4], \
                  'AC': [-8.4, -22.4], \
                  'CT': [-7.8, -21.0], \
                  'AG': [-7.8, -21.0], \
                  'GA': [-8.2, -22.2], \
                  'CG': [-10.6, -27.2], \
                  'GC': [-9.8, -24.4], \
                  'GG': [-8.0, -19.9], \
                  'CC': [-8.0, -19.9], \
                  'TC': [-8.2, -22.2]}


    s= s.upper()

    # Molar concentration of all monovalent ions
    # (Owczarzy R, Biochemistry, 2008)
    Mon= (Na + K + Tris / 2.0) / 1000.0

    # Molar concentration for non-self complementary oligonucleotide
    # (SantaLucia J, PNAS, 1998)
    Ct= (dna / 4.0) * 1e-9

    # Molar concentration of free Mg+2
    # (Owczarzy R, Biochemistry, 2008)
    fMg= (Mg - dNTP) / 1000.0

    if Mon > 0:
        ionR= fMg ** 0.5 / Mon
        #print 'Ion ratio is', ionR

    # calculate faction of GC
    fGC= float(s.count('G') + s.count('C')) / len(s)

    # DMSO correction
    # take the average of all four coefficients previously published
    # (von Ahsen N, Clin Chem, 2000)
    DMSO_coeff= (0.75 + 0.6 + 0.675 +0.5) / 4
    # 1% DMSO reduces Tm by DMSO_coeff
    dmso= DMSO_coeff*DMSO

    # double check for the 'fit-model' for criteria and range
    if Mon == 0 or Mg == 0:
        pass
    elif 1e-10 <= (dna*1e-9) <= 0.1 and 7 < len(s) < 61 and 0.015 <= Mon <= 1.2 and 0.01 <= Mg < 600 and dNTP <= (Mg*1.2):
        pass
    else:
        print 'buffer/DNA conditions out of range for calculation'
        return


    # collect the NN values for dh and ds
    dh = list()
    ds = list()
    # First, terminal base corrections
    if s[0]=='G' or s[0]=='C':
        dh.append(0.1)
        ds.append(-2.8)
    elif s[0]=='A' or s[0]=='T':
        dh.append(2.3)
        ds.append(4.1)
    if s[-1]=='G' or s[-1]=='C':
        dh.append(0.1)
        ds.append(-2.8)
    elif s[-1]=='A' or s[-1]=='T':
        dh.append(2.3)
        ds.append(4.1)

    # now the rest of the sequence
    for keys in deltaHS:
        H, S = deltaHS[keys]
        dh.append(olcount(s, keys) * H)
        ds.append(olcount(s, keys) * S)

    Rconst= 1.987 #gas constant in cal/mol K
    # this is the deltaG value calculated in Kelvin
    tmNN= 1000 * sum(dh) / (sum(ds) + (Rconst * log(Ct)))


    # now correct for ion composition
    # Decision tree based on Owczarzy R, Biochemistry, 2008
    if Mon == 0 or (Mon != 0 and ionR >= 6.0):
        # incorporate Mg correction
        # (Owczarzy R, Biochemistry, 2008)
        # various constants (1/K)
        Va= 3.92e-5
        Vb= -9.11e-6
        Vc= 6.26e-5
        Vd= 1.42e-5
        Ve= -4.82e-4
        Vf= 5.25e-4
        Vg= 8.31e-5

        # calculate tm in K
        Ktm= (1 / tmNN) + Va +\
             (Vb * log(fMg)) +\
             (fGC * (Vc + Vd * log(fMg))) +\
             ((Ve + Vf * log(fMg) + (Vg * (log(fMg)) ** 2)) / (2 * (len(s) - 1)))

    elif Mon != 0 and ionR < 0.22:
        # incorporate salt correction (for all monovalent ions)
        # (OwczarzyR, Biochemistry, 2004)
        # calculate tm in K
        Ktm= (1 / tmNN) +\
             (4.29 * fGC - 3.95) * 1e-5 * log(Mon) +\
             9.40 * 1e-6 * log(Mon) ** 2

    elif Mon != 0 and 0.22 <= ionR < 6.0:
        # incorporate Mg correction
        # (OwczarzyR, Biochemistry, 2008)
        # various constants (1/K) corrected for Mon+ as well
        Va= 3.92e-5 * (0.843 - (0.352 * Mon ** 0.5 * log(Mon)))
        Vb= -9.11e-6
        Vc= 6.26e-5
        Vd= 1.42e-5 * (1.279 - 4.03e-3 * log(Mon) - 8.03e-3 * log(Mon) ** 2)
        Ve= -4.82e-4
        Vf= 5.25e-4
        Vg= 8.31e-5 * (0.486 - 0.258 * log(Mon) + 5.25e-3 * log(Mon) ** 3)

        # calculate tm in K
        Ktm= (1 / tmNN) + Va +\
             (Vb * log(fMg)) +\
             (fGC * (Vc + Vd * log(fMg))) +\
             ((Ve + Vf * log(fMg) + (Vg * (log(fMg)) ** 2)) / (2 * (len(s) - 1)))

    else:
        print "error in Tm calculation"
        return

    # return Tm in degrees C
    return 1 / Ktm - 273.15 - dmso




if __name__ == '__main__':
    # a test case
    print Tm("ACGGAGAGGGATGGCATG", dna=10, Na=200, dNTP=0.1, Mg=10.0, Tris=10, K=10, DMSO=5)