SeqCap Epi連載[5]｜BSMAP methratio.pyでメチル化検出

先日鈴鹿サーキットに初めて行ってきた久保（kubor）です。レーシングカーのスピードもさることながら、時の流れも早いもので、僕がアメリエフに入社して2ヶ月が経ちました。これからもどうぞ、よろしくお願いいたします。

さて、SeqCap Epi連載第5回目の今回は、BSMAP^*1でマッピングした結果から、メチル化塩基を検出するPythonスクリプト（methratio.py）を紹介します。BSMAPは次世代シーケンサー（NGS）でバイサルファイトシーケンスを行った結果得られるシーケンスリードをマッピングするためのツールソフトウェアです。このBSMAPについては、前回「SeqCap Epi連載[4]｜BSMAPでバイサルファイトシーケンスのマッピング」こちらの記事で紹介致しました。ご興味のある方は併せてお読みいただくと、より解析手順を把握しやすいと思います。

｜methratio.pyはBSMAPディレクトリに含まれています

BSMAPには、あらかじめ、メチル化塩基検出用のmethratio.pyが同梱されています。こちらを使用することにより、メチル化塩基検出結果をBSMAP独自のフォーマットとして得ることが可能です。 methratio.pyの格納場所は、BSAMPのルートディレクトリです。実行にはPython 2.x系が必要です。

｜使い方

使い方は、非常にシンプルです。必要なのは、オプションと、入力ファイル。入力ファイルは、BSMAPでマッピングしたBAMファイルを用います。弊社のSeqCap Epi解析パイプラインにおいては、オーバーラップリードおよび重複リードの削除を行ったBAMファイルから、メチル化検出を行います。こちらがコマンド例です。

$ python methratio.py [オプション] BSMAPの出力.bam

｜フィルタリングを設定する

このmethratio.pyには複数のフィルタリングオプションが用意されています。そのうち特徴的なものが以下のオプションです。

-m &ltint&gt : メチル化を検出する領域について、指定した値以上のdepthがある場合に結果を出力する
-z : メチル化率が0になる領域のシトシン（C）についても結果を出力する
-c &ltCHR&gt : 染色体を指定してメチル化コールする
-i [skip, correct, no-action] : SNPを考慮するかどうか設定する

skip - C->TのSNPの可能性（逆鎖でAが検出された場合）がある領域を結果に出力しない
correct - C->TのSNPの可能性がある場合でも、その領域を結果に出力し、同時に逆鎖のGの数およびG+Aの数を出力する
no-action - SNPの可能性を考慮せずにメチル化コールをする（逆鎖をチェックしない）

これらの中でも-iオプションが特徴的です。例えば弊社のパイプラインでは、SNPの可能性がある領域を結果に出力しないよう、-i skipを指定しています。この理由については、次回、詳しく解説したいと考えています。 methratio.pyでSNPを考慮しないように設定したことから、本パイプラインにおけるSNP検出には、この後紹介予定のBisSNP^*2を用いています。BisSNPは既知のSNP情報も利用したSNP検出とメチル化コールが可能なため、より正確なSNP検出結果を期待できます。

｜出力フォーマットを知る

methratio.pyの出力フォーマット例です。

chr     pos     strand  context ratio   eff_CT_count    C_count CT_count        rev_G_count     rev_GA_count    CI_lower        CI_upper
NC_001416       4       +       GGCGG   0.000   4.00    0       4       0       0       0.000   0.490
NC_001416       7       +       GGCGA   0.000   5.00    0       5       0       0       -0.000  0.434
NC_001416       10      +       GACCT   0.000   5.00    0       5       0       0       -0.000  0.434

項目が多いですが、上述の-iオプションでskipを指定している場合、注目するべきは塩基のポジションの他に、2点です。まず1つがC_countです。メチル化Cの数はC_countに相当します。次に、CT_countはそのポジションで検出されたすべてのC+Tの数を示します。したがって、5列目のratio（メチル化率）は、C_count/CT_countで表されます。これらの項目の詳細およびその他の項目については、READMEで確認出来ます。注意したいのは、このフォーマットはBSMAP独自のものであるため、下流の解析で用いるソフトウェアがこのフォーマットに対応しているかどうか確認が必要であるという点です。

｜編集後記

解析結果がオリジナルのフォーマットというだけで少々厄介ですが、なによりもフォーマットの説明がどこに明記されているのかということが問題になることが多いです。BSMAPのように、配布サイトでの説明がなく、READMEに記載されているものも多いです。そういう僕も、過去に何度も説明文を探すことに骨を折りました。そうした経緯もあって、最近では、READMEを見つけた際には、とにかく読んでみる癖がつきました。今となっては当たり前ですが、Linuxに慣れるためには大切なことですね。

SeqCap Epi連載シリーズ過去記事一覧

*1:bsmap - Bisulfite Sequence Mapping Program, BSMAP配布サイトにて公開中 *2:Bis-SNP, Bis-SNP配布サイトにて公開中

アメリエフの技術ブログ

Amelieff Staff Blog