[BioRuby-ja] Fwd: Bio::Alignment::EnumerableExtension#consensus

[Naohisa GOTO]

後藤＠阪大です。

On Thu, 26 Jul 2007 20:50:46 +0900
"遠藤大二" <dendoh ＠ rakuno.ac.jp> wrote:

> degenerated nuclotide は、consensus_iupac メソッドを使ってください。
> 
>   ffconsensus = ffa.alignment.consensus_iupac
> 
> これで、今回障害になっていた問題は全て解決しました。
> .consensus_iupac　のメソッド一つを理解できるか否かで2日分の作業が違ってしまいます。
> 
> とりあえず、お礼申し上げます。

お役に立てて幸いです。
ドキュメントが、いまいち、わかりにくいのは、徐々に直していきたいと思います。

> 加えて、質問なのですが、mafftでは"-"で示されている
> ギャップがconsensusでは?にされています。この文字の選択については背景に意味があるのでしょうか。

BioPerlの実装に合わせただけで、深い意味はなかったと思います、たぶん。

実は、consensus_iupac は以下のようなギャップの取扱いを指定するための
オプションを取ります。consensus_string も同様のオプションを取りますが、
consensus_stringの第一引数は閾値である点は注意してください。

consensus_iupacでは、 :gap_mode => 0 と :gap_mode => 1 の
場合で、文字が?と-である点以外は違いがわからないのですが、
閾値の計算が絡む consensus_string では違いがわかると思います。

  a = Bio::Alignment.new([
      'aaaaaaaaccccgg',
      'aaccggttggtttt',
      'a-a-a-a-c-c-g-' ])

  ### consensus_string で閾値 1 と 0.5 の場合の例

  # ギャップも普通の文字として扱う(ギャップを閾値計算時の分母分子両方にカウント)
  # これがデフォルトの動作
  p a.consensus_string(1.0, {:gap_mode => 0})
  # "a?????????????"
  # 閾値0.5
  p a.consensus_string(0.5, {:gap_mode => 0})
  # "aaa?a?a?c?c?g?"

  # ギャップがあるサイトはギャップとして扱う(ギャップがあるサイトは必ず-になる)
  p a.consensus_string(1.0, {:gap_mode => 1})
  # "a-?-?-?-?-?-?-"
  p a.consensus_string(0.5, {:gap_mode => 1})
  # "a-a-a-a-c-c-g-"

  # ギャップは計算から除外する(ギャップは閾値計算時に分母分子両方から除外)
  p a.consensus_string(1.0, {:gap_mode => -1})
  # "aa????????????"
  p a.consensus_string(0.5, {:gap_mode => -1})
  # "aaaaaaaaccccgg"

  ### consensus_iupac の例

  # ギャップも普通の文字として扱う(よって、ギャップがあるサイトは?になる)
  # これがデフォルトの動作
  p a.consensus_iupac({:gap_mode => 0})
  # "a?m?r?w?s?y?k?"

  # ギャップがあるサイトはギャップとして扱う(ギャップがあるサイトは必ず-になる)
  p a.consensus_iupac({:gap_mode => 1})
  # "a-m-r-w-s-y-k-"

  # ギャップは計算から除外する
  p a.consensus_iupac({:gap_mode => -1})
  # "aammrrwwssyykk"


> degenerated nuclotide は、mafft alignment で一件でも　−　
> が有ると　? になってしまいます。alignment の状況によっては　−
> があっても善意に解釈したいときもあります。

という場合には、
  a.consensus_iupac(:gap_mode => -1)
を使えば、とりあえずは、大丈夫だと思います。

もちろん、ギャップが何個以下だったら救出するが、それ以上だったら
ギャップ扱い、のような複雑なことをしようと思えば、自分独自の
ルーチンを組む必要があります。

このとき、each_site というイテレータを使えば便利と思います。

  # さきほどのプログラムの続きとします
  a.each_site { |x| p x }

これを実行すると
["a", "a", "a"]
["a", "a", "-"]
["a", "c", "a"]
["a", "c", "-"]
["a", "g", "a"]
["a", "g", "-"]
["a", "t", "a"]
["a", "t", "-"]
["c", "g", "c"]
["c", "g", "-"]
["c", "t", "c"]
["c", "t", "-"]
["g", "t", "g"]
["g", "t", "-"]
と表示されます。つまり、各サイトの内容をArray (より正確には、Arrayを
継承した Bio::Alignment::Site クラスのオブジェクト)に入れた内容で、
繰り返し実行してくれます。

-- 
後藤 直久  ngoto ＠ gen-info.osaka-u.ac.jp
大阪大学微生物病研究所 遺伝情報実験センター ゲノム情報解析分野(安永研)
Phone: 06-6879-8365 / FAX: 06-6879-2047