[Ruby][BioInformatics] BioRubyで核酸配列のTmを計算する

ちょっと、必要があって、Rubyでmethodとかを書く練習も兼ねて、書いてみた。
BioRuby-1.4.1のBio::Sequence::NAを継承して、tmを計算するmethodを定義する。

#!/opt/local/bin/ruby

require 'bio'

class MYNA < Bio::Sequence::NA
 # Calculate the melting temperature
 #  reference: http://www.sigmaaldrich.com/japan/lifescience/custom-products/custom-dna/oligo-learning-center.html#o04
 # 
 # s = Bio::Sequence::NA.new('CGACGTTGTAAAACGACGGCCAGT')
 # puts s.tm                            #=> 73.15
 # 
 # * (optional) _ct : conc of oligo nucleotide(mM) (default 0.5)
 # * (optional) _na : conc of Na+ (mM) (default 50)
 # * (optional) _fa : conc of Formamide (mol/L) (default 0)
 # *Returns*:: Float
 def tm(_ct=0.5, _na=50, _fa=0)
   naseq = self.dna
   len = naseq.length.prec_f # base length
   ct = _ct.prec_f/1000000   # concentration of oligo (mol/L)
   na = _na.prec_f/1000      # concentration of Na+ (mol/L)
   fa = _fa.prec_f           # concentration of Formamide (mol/L)
   @@r=1.987                   # gas constant R
   
   # to_upper
   naseq.upcase!

   # thermodynamic parameter 
   #  _h ...delta enthalpy
   #  _s ...delta entropy
   _h = {"AA" => -9.1, "TT" => -9.1, "AT" => -8.6, "TA" => -6.0,
            "CA" => -5.8, "TG" => -5.8, "GT" => -6.5, "AC" => -6.5,
            "CT" => -7.8, "AG" => -7.8, "GA" => -5.6, "TC" => -5.6,
            "CG" => -11.9, "GC" => -11.1, "GG" => -11.0, "CC" => -11.0}
   _s={"AA" => -24.0, "TT" => -24.0, "AT" => -23.9, "TA" => -16.9,
          "CA" => -12.9, "TG" => -12.9, "GT" => -17.3, "AC" => -17.3,
          "CT" => -20.8, "AG" => -20.8, "GA" => -13.5, "TC" => -13.5,
          "CG" => -27.8, "GC" => -26.7, "GG" => -26.6, "CC" => -26.6}
   tot_h=0.0
   tot_s=0.0
   for i in 0..(naseq.length-2)
     tot_h += _h[ naseq.slice(i,2) ]
     tot_s += _s[ naseq.slice(i,2) ]
   end

   count = self.composition
   at = count['a'] + count['t']
   gc = count['g'] + count['c']

   _tm = ((1000*tot_h)/(-10.8+tot_s+@@r*Math::log(ct/4)))-273.15+16.6*Math::log10(na)
   if  _tm > 80
     # gc% method
     tm = 81.5+16.6 * Math::log10(na) + 41*((gc.prec_f)/len)-500/len-0.62*fa
   elsif _tm < 20
     # wallace method
     tm = 2*(at)+4*(gc)
   else
     tm = _tm
   end

   return tm
 end
end

dna = MYNA.new("CGACGTTGTAAAACGACGGCCAGT")
puts dna.tm  # =>73.1531598741961

Tmの計算方法自体は、http://www.sigmaaldrich.com/japan/lifescience/custom-products/custom-dna/oligo-learning-center.html#o04 に紹介されている方法をそのまま書いてみた。

次にやること

Tmの計算方法は、このmethodにも書いたように一つではないので、methodに引数を指定することで、どの方法で計算するのかを指定できるようにする。
MeltCalc( http://www.meltcalc.com/ )の論文に書かれているパラメータと式も使えるようにする。
そもそも、BioRuby的に、こういうmethodの書き方が正しいのかどうか良く解らんので、ちゃんと勉強する。