BioPythonまとめ その2
その1が長くなったのでその２です。
～NGS処理ツール～
まずはおなじみSAM・BAMファイルからの抽出をpythonでもやってしまおうというもの。できるんですね。Linuxでやると調べるのに苦労するからpythonで実装してくれていてとてもありがたい。ちなみにこいつらもともとC⁺⁺言語用だったようで、これらを覚えているとサクッといけたらしい。pythonでも実装しようということで、python処理用の薄い皮で包み込んだpysamを誰かが作ってくれたのです。ありがたい。
と思ったら事件。Windows10ではpysam使えないらしい。インストールできない。はあ。。。
仮想Linux環境作っておいてほんと良かった。。。
これからSAMファイルとかBAMファイル使うときにpythonいじるときはこっちからになりそうだな。。
一応BAMファイルのパーサーとしてbamnosticがWindowsに対応しているらしいが。。。本当かな？
bamnosticはpipから無事インストールできました！使い方は下参照。
Welcome to bamnostic’s documentation! — bamnostic 1.1.4 documentation
これでSAMファイルもpythonからいじれるようになる。。。かな？
一応ubuntuでpysamもやっとこう。ダウンロードの仕方などはこちら。
pysam does not install on windows · Issue #575 · pysam-developers/pysam · GitHub
ubutuでpysam自体はインストールできたが、bamfile=pysam.AlignmentFile('bamファイル名’,'rb')とやっても、could not retrieve index file for 'ファイル名'と出てきてファイルの読み出しができない。普通こういう時はパスが通ってない時だが、そもそもbioinfoのディレクトリにいるのだからパスもくそもない。。。
むむむ。。。ファイルが大きすぎるのだろうか。import osからos.listdir()で現在のディレクトリ内のファイルを確認しても目的のファイルは列挙されるし。。。
わからん！いかんせん仮想環境内だから対策も練りようがない。
ここは大人しくbamnosticにいこう！さらばpysam！君には期待していたのだが！
ちなみに、覚書なのだが、samtoolsはもともとC/C⁺⁺ベースで作られたツールのため、samtoolsをいじる際に出てくる-oとか-uとか、こういうのは全部C/C⁺⁺にしたがっているぽい。今後調べるときの参考に。
お次はbamnostic。コマンドプロンプトでうまくできなかったので、一度諦めてjupyternotebookにアップロードしてから動作を実行した。パスは一応デスクトップに通してあるはずなのだが、なぜだろうか。とりあえず今沼にはまりかけているので、できるところからやる。
bamnosticの使い方はpysamと大体同じっぽい？データの読み出しの時に、ランダムアクセスはできませんうんちゃらかんちゃらとの警告が出たが、無視してそのままヘッダーを表示。(.header)よし、うまくいった。
ubuntuでできなかったのも、ランダムアクセスが制限されているからなのでは？と考えて、ubuntu でもやってみた。一応。
すると。。。なんだ、アクセスできるじゃん！
流れとしては、pysamを使ってAlignmentFileからbamファイルを開いてみたが、could not retrieveとかでファイルにはアクセスできない風だった。
bamfile = pysam.AlignmentFile('unmapped.bam','rb')
が、次の行で
print(bamfile.header)
を使うと下の行にズラーッとヘッダーが表示されてきた。つまり仮想Linux環境下で、カレントディレクトリのbamファイルにアクセスできた、というわけである！
ここからは推測だが、AlignmentFileからアクセスできなかったのは、たぶんbamnositcと同様にランダムアクセスをしてしまっていたから。headerなど指定してアクセスする分にはうまくいく。ということで今後はAlignmentFileを打ち込んだのちに起こるエラーは無視することにしよう。
またここでヘッダーとは、リファレンスとは、などと疑問に思ったので、SAMファイルの構造から理解することにした。
参考となるURLは↓
crusade1096.web.fc2.com/sam.html
要はSAMファイルには下のように11列の配列情報が入っているということ。
1列目：リードの名前
2列目：リードの状況（どんな風にマッピングされてるか判断できる，以下参照）
3列目：張り付いた染色体，コンティグの名前
4列目：張り付いた場所
5列目：マッピングスコア
6列目：マッピング状況（indel入ってるとか，どのくらいマッチしてるとか判断できる，以下参照）
7列目：paired-endの時の相方の名前
8列目：paired-endの時の相方の場所
9列目：paired-endの時のインサートの長さ
10列目：リードのシークエンス配列
11列目：リードのクオリティ
さて、構わずbamnosticやpysamで処理を続けていたが、ここで問題が生じてきた。ヘッダー情報やレングス情報など、ただファイルから情報を取り出すだけの処理ならば問題なくできるのだが、例えばpileupやsortなど、がうまくできないなあ。
調べてみるとbamnositicにはそもそもsort機能がついてない。sort以外の機能、読み出しなどについてはpysamと同じ感じなのか。ただbamnositicはランダムアクセスをするとエラーを返してくるようなので、for構文を使ってうまく順々にファイルアクセスするようにしてあげなければならないぽいなあ。。。bamnosticはNGS処理では使い勝手が悪いか。。。
お次pysam。これはまだ未知数だが、windowsでは仮想環境上に構築しているため、いかんせん使い勝手があまりよくないのでござる。。。
ただし、sortとindexは成功したようだ。
unmapped.bamに対して、pysam.sort('-o','unmapped.sorted.bam','unmapped.bam')を実行したらsortファイルであるunmapped.sorted.bamがカレントディレクトリ上に出現した。すなわち成功。さらにpysam.index(sortファイル)を行ってもエラーが出現しなかったため、これも成功ということでよかろう。
ちなみにindexが無事終わると、.bam.baiファイルが出現する。.baiファイルを作ることをbamファイルにindexを付与すると言うのだろうか。このbaiファイルはマッピングデータを視覚的に表示できるIGVファイルで用いるときに必要なようである。他の用途はわからん。。。下のサイトに色々とIGVの使い方が書かれている。
ただしpysamで問題なのはfetch(出現させる)である。一度やってダメだったが、indexを付与したファイルを用いて、
samfile = pysam.AlignmentFile('unmapped.sorted.bam','rb')
を打ち込んでみる。すると。。。
なんと何のエラーもなく実行された！！！
これらのことから、普通にbamファイルを上のコードを起動するとエラーが起こるのは、indexが付与されていなかったからであり、index付与後には問題なく扱えることがわかった。
bamnosticよりもpysamの方がよさそうである。これから実際の解析にはこれを使おう。（先ほどbamnosticにはソート機能がないと書いたが、bam.baiにできるところがありそうなので、詳しく調べればソートできるところもあるのだろうか？）
さて、pysamをだいぶ使いこなしてきました。ソート後インデックスを付与したファイルを対象（今回はunmapped.sorted.bam）に、
countread = samfile.count()
print(countread)
と打ち込むとカウント数がでた！これでBAMファイルにいくらリード数が入っているかわかるね！ついでにcountの()のなかにcontigや範囲を指定すればそこにおいてマッピングされたリードの数を数えられます。
とりあえずこれでpysamは終わり。samtoolsのviewは使えるぽいけど、bamファイル同士を結合させるmergeはないのか？そしたらsamtools使用に逆戻りだなあ。。。
と思ったらpysam.merge(""all.bam"",''入れたいファイル１","入れたいファイル２",,,,,,)とやれば入るっぽい。
unmappedも取り出せるらしいし、先生に言われた処理全部パイソンでできそう。多分表記方法さえ変えれば全部pysamで代用できるっぽいな。すごいぞpysam！今日一日かけた甲斐があった。pileupは必要なときにまた。
GO解析に移ったが、正直何がなにやら
とりあえずgoatoolsをインストールすればよいのかしら。GO解析の入出力は主にエクセル形式で行われるぽい？
また、jupyternotebookにサンプルコードを残しているわけだが、まあ必要になったら参照する程度でよいかと、、、かなりムズイ。
一つ収穫は今回のように、pydotなど、graphvizを使ったものはdot.exeという実行ファイルにパスが通ってないと使えないということ。今回もjupyternotebookで実行するとうまくパスが通らずに'File Not Found Error'、'dot' not foundと表示された。
jupyternotebookでは結局うまくできなかった。コマンドプロンプトでは実行してないが、dot.exeにもパスは通してあるのでうまくできる。はず。これにこだわってもしょうがないので一回スルーしよう。
