[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の書き方が正しいのかどうか良く解らんので、ちゃんと勉強する。