RNA発現解析
次にマッピングしたデータを処理する。
参照配列の塩基にどれだけのリードがあったかを集計・整理する。リードを集計・整理できたら（リードカウント数として確認出来たら）そのカウント数をTPMなどで補正してそれぞれの出現量とする。
さらに出現量の違いを距離として入力したサンプル間で系統樹を書いたりなどの処理が可能となる。
また発現量の違いが大きい遺伝子を抽出すればそれらの生物的な関連をジーンオントロジーやパスウェイに当てはめて考察できる。
まず必要なことは参照配列のgffデータに対して遺伝子ごとに通番のgene_idを挿入する。これはのちの作業のインデックス化に必要となる元データになる。
pip install bcbio-gffでbcbio-gffをインストールしたのちに元のgffファイルから作る。
with open('gffファイル名',"w") as (name):
でgffファイル名を開き、GFF.parseメソッドを用いてfor構文で一行ずつ読み出す。
そのさいにfeaturesをfor構文で読みだしてfeature.typeからgeneとpseudogeneに一致するものの、gene_idを0001から連番でつけていく。f.qualifiersで名づけたら次にsub_featuresも同様に、さらにそのsub_featuresのgene_idも同様に名付けていく。
最後にGFF.writeメソッドで出力して終了。
これでgene_idを連番にした参照ファイルを作ることができた。
次にこれからリードカウント数を数えるためにsbreadパッケージをダウンロードしてfeatureCountsを使って作業を行った。
The Subread package (sourceforge.net)
のLatest Version2.0.1のリンクからダウンロードした。
ダウンロード後はsudo apt install tarで
tar -xvcf
でtar.gzファイルを解凍して、
目的ファイルのあるディレクトリまで下りたらmake -f Make.Linuxと打ち込んで完了。あとはsudo apt install subreadでfeatureCountsが使える状態となった。
featureCountsはオプションが色々とあり、-pでペアエンドの指定、-T （数字）でスレッド数の指定、-t exonでexonを数えることを指定、-g （アトリビュート）でグループ化の際のアトリビュートを指定できる。また-a （アノテーションファイル）でアノテーションファイルを指定する。これでアトリビュートごとのカウント数を得ることができる。
使い方の詳細は公式サイト
http://bioinf.wehi.edu.au/subread-package/SubreadUsersGuide.pdf
を参照のこと。
一通り処理が終わったらリードカウント数のデータがまとめて入ったファイルがtsv形式で出てきた。RNA発現量解析はtsvデータをいじって行うらしい。
とまあこんな感じでいったんデータの集計は終わり。
次からこのtsvファイルを使ってTPM値を使ってデータの補正を行う。
あと今参考にしている本がさらに参考にしてたサイト見つけたのでのっけとく。
これさえ見ればpython使ったNGS解析は一応できるんじゃ、てほど丁寧。
動画付き。
情報解析講習会ビデオ＜2019年度第1回PAGS・DDBJ・DBCLS合同情報解析講習会＞ ｜ 先進ゲノム支援（先進ゲノム解析研究推進プラットフォーム） (genome-sci.jp)
ここで話は変わるが宿主と寄生生物の遺伝子の分割、SVMを用いたwebツールを使えば一発ととある論文に書いてあったので楽観視してこれまで調べなかったのだが、今日そのサイトをみたらサービス終了してた。。。
そもそもその論文自体が2005年に発行されたものだから仕方がないっちゃないけども、また論文探しから始めないといかん。
一応使えそうなツールとして、まだしっかりと論文を読んではいないがblobtoolsなるものが(運がいいことに)pythonのツールとして2018年にGithubで出されてた。
これ使えばなんとかなるかもしれないが、さすがに機械学習をさせたいので、Webツールに頼りたい。。。
論文をまた探してみるか。
ubuntuでAnaconda3を入れた話
はい、とうとうubuntuでAnaconda3を入れました。
別にubuntuではAnacondaの縮小バージョンminicondaでいいかな、と思ってたけど、minicondaからだとどうもpythonがバージョン2しか稼働しないぽく。。。
するとpython3以降からのみ対応しているパッケージが全部minicondaによって邪魔されてインストールできないor稼働しない。。。
新しいバージョンのminicondaをインストールする手もあったんだろうけど、linux置いてるところからjupyternotebook使えるのも便利だし、そこもう統一しようかということで入れちゃいました。
割と四苦八苦したので書いておく。
まずネットに公式サイトからLinux用のAnacondaをインストールしてくださいって書いてたんだけど、そのままインストールするだけじゃ多分無理。
ubuntuの置いてあるファイルがWIndowsの環境から独立しているので、インストールしたところでubuntuのとこまで正常に持っていけない。
なのでubuntuのコマンドから
wget (ダウンロードのurl)
でダウンロードするのが無難。
無事ダウンロードできたら
bash Anaconda3-2020.11-Linux-x86_64.sh
でダウンロードしたAnacondaファイルを開く。このときの2020.11は今後Anacondaのバージョンアップに伴い変化するので、そこは適宜変更する。
ファイルを開いたらあとはubuntuの画面に出てくる通りにエンターとかyesとか押してけばインストールは完了。
でも大体はここでパスが通ってない。
Thank you for installing
みたいな文字が出て打ち込み画面が現れたら
conda -V
でcondaにパス通っているか確認。パス通ってなかったらエラーが出てくる。
なのでubuntuのユーザーディレクトリにある.bashrcのファイルを開いて自分で直接パスを書き込む。ubuntuなので直接書かなきゃいけないのがめんどくさいね。WIndowsなら環境変数でいけるのに。ちなみに.bashcとかそういうファイルを開くためにAtomをインストールしとくとよい。このアプリ見やすくてお勧め。
.bashrcの最後の行にexport PATH=~/anaconda3/bin:$PATH
と打ち込んでsaveしたらubuntu に戻って
source ~/.bashrc
と入力してbashrcファイルを読み込んで終了。
これでパスも通ったので
conda -V
と打てばバージョン情報が出てくるはず。
あとubuntuからjupyter notebookも使えて、
jupyter notebook --no-browser
と打てばリンクを返してくれるからそれをweb検索エンジンに打ち込めばhomeディレクトリ内のフォルダが既に入ってる状態のjupyternotebookを返してくれる。優秀。
anacondaが無事入ったのですかさずpip install biopython、pip install bcbio-gff、conda install hisat2でインストール。無事入りましたとさ。
こういうもので壁にぶち当たったり沼にハマったときのひとつの解決策として、確かにそのエラーなりの解決策を考えるのも重要だけど、そんなに時間もないわけで。
もう沼から抜け出せねえ、てなったらいっそそのパッケージなりごと消してもう一回インストールしなおす、というのも重要かもしれん。根本からぶっ壊す、みたいな。
ついにubuntuに移行した
これまでjupyternotebook便利すぎてwindows環境にこだわってたけど明らかにLinuxのが使いやすいしwindows対応してないパッケージ多すぎるしで結局RNAの発現解析調べるのはubuntuでやることにした。
まずBiocondaをインストールなどもろもろ環境整備してからfastqcをダウンロード。これで自由にfastqcファイルを見れるぞ！（これまではテキスト形式でみてた）
さて、バイオインフォマティクスを学んできて一か月、やっとお目当てのRNA発現量解析のところまできた。（長かった。。。）
今回はマッピングから。
ubuntuでBiocondaをインストールして、そこからhisat2を使って行う。
マッピングにはまず参照配列をインデックス化する必要がある。インデックスとはソフトウェアのマッピングプログラムから参照しやすくなるように整理・加工したもの。参照配列は必ずインデックス化する必要がある。
hisat-2を用いてインデックス化する際にはhisat2-build (ファイル名).fna (新しいファイル名).fnaでよい。
次にhisat2を用いてマッピングをする。
hisat2 (パラメータ指定) -x（参照配列） -1（入力配列） -2 （入力配列その２） -s（出力配列）
のようにマッピングする。マッピングのパラメーターには様々なものがある。パラメータ指定については
Manual | HISAT2 (daehwankimlab.github.io)
を参考にすればよい。
んでいざマッピング、と思ったら
Error, fewer reads in file specified with -2 than in file specified with -1
terminate called after throwing an instance of 'int'
Aborted (core dumped)
(ERR): hisat2-align exited with value 134
というエラー。どうやら2つ目のファイルがうまくトリミングできていないっぽい？
のでもう一度trimmomaticでトリミングしてから同じコードでマッピング。したらうまくいった。どうやらさっきやったトリミングはうまくできていないらしい。
ちなみにこのサイト覗いてみると、どうしてもread数がそろわないけどマッピングしたいときにはhisat2ではエラーが出てしまうのでSTARを使えばいいらしい。
Paired-endのデータが不揃いの場合に起こるエラーについて - Qiita
無事samファイルを出力しようと思ったのでbamファイルを作ろうと思い、
samtools view -@ 32 -b (ファイル名sam) -o (ファイル名bam)
で実行したら事件。
samtools: error while loading shared libraries: libcrypto.so.1.0.0: cannot open shared object file: No such file or directory
とまたもやエラー。Biocondaでsamtoolsをインストールするとlibcrypto.so.1.0.0を認識できずにこうしたエラーが起こるらしい。
色々ネットに解決策が転がっていたが、Ubuntuのアップデートによりopensslがアップグレードされてしまい、かなりめんどくさい状況。opensslをダウングレードすればいいっぽいが、すると他の色々もダウングレードされるので
conda uninstall samtools
でsamtoolsをアンインストールして、
sudo apt install samtools
でsamtoolsを再インストール。逃げの一手。笑。で、結局うまくbamファイルを作れました。めでたし。
samtools sort -@ 32 （もとファイルbam） > （ソートするファイル.sorted.bam）
以上のようにsortしてbamファイルの準備は完了。このソート作業でリードはマップされたリード開始位置の順になる。
またsortの際にmerging from 0 files and 32 in-memory blocks...と出るがどうやらエラーではないっぽい。sortするときにreadを一時的に保存した場所を知らせているだけらしい。
numpy np.r_/np.c_の使い方
numpyのnp.r_とnp.c_が便利そうだったのでまとめます。
下の記事の内容まんまです。
多次元配列の結合を行うオブジェクトnp.c_とnp.r_の使い方 - DeepAge
配列を結合したい際、次元数が違う時がある。そこでnp.r_を使えばちがう次元でも結合可能。
数値の文字列で表すことができます(a, b, cはそれぞれ整数)。デフォルトでは'0, 0, -1'となっています。
aは、どの軸(axis)の方向に沿って配列を結合するのかを指定します。bは、できあがる配列の次元数の最小値を指定します。cは、次元数の少ない配列の次元数の拡張を行った際、形状(shape)の表記として、どこに配列の最後の次元が置かれるべきかを指定します。
cについては、例えば (2,3) の形状(shape)を持つ二次元配列があったとします。この配列を他の(2,2,3)の形状(shape)を持つ三次元配列と結合したいとします。このとき、二次元配列は三次元に拡張しなければ結合できません。
c=-1またはc=1であれば、先頭に1が追加され、(1, 2, 3)の形状で結合することができます。c=0であれば、元の配列の形状が先頭に詰められるので、(2,3,1)となります。
上のリンクに書いてあったもの。np.r_に上の三種類の条件値を指定すれば結合の方法を指定できる。
それでnp.c_はこのうちa=-1,b=2,c=0を指定して入れる。
a=-1というのは最も低い次元のもとで結合していく。c=0は次元の違うものを結合させる際、末尾に1を追加するというもの。
ちなみにaxisとは座標軸のもの。行と列のみの行列なら列がaxis=0,行がaxis=1となる。
つまり、任意のx,yについてnp.gridmeshで行列をそろえ、np.ravel( )でリストにし、np.c_[ ]で結合すれば[x,y]の要素が入ったサイズが列数の行列を手に入れることができる。
WindowsからGraphviz利用
Graphvizを使うのに手間取ったので覚え書き。
決定木を書こうと思い、pydotplusを使いデータの読み込み、pydotplus内のgraphvizでデータの書き出しをしようとしたがここでエラー。
graphに決定木の学習データclfを入れてgraph.write_png('DecisionTree.png')でグラフを書こうと思ったがエラーが。
pydotplus.gaphviz.InvocationExceptionと表示され、pngファイルを認識できません、Use one of : でそのまま何も出ず。
graphvizはpip install graphvizでインストール済み。
これはgraphvizがただインストールしただけじゃダメな奴かな、と色々調べたところどうやらWindowsではコンパイラに何か問題があるぽい。やっぱMacほしいなあ。。。
しょうがないのでpip uninstall graphvizでgraphvizをuninstallしてgraphvizの公式サイトから直接ダウンロード。zipファイルを解凍してPathを通してリトライ。
コマンドプロンプトで dot -vと打ち込んでバージョン情報が出てくればgraphvizがインストールされている状態。無事バージョン情報が出てきたので、よっしゃこれでいけるわ！と思い再びpngファイルを書き出してみる。
エラー。
また同じエラー。嘘じゃん。
エラー内容にdot -cを入れてみたらどうですか、とあったのでコマンドプロンプトに戻りdot -cを入れる。
これが入らん。問題はdot -cが無事実行されれば解決されるだろうと思い、公式サイトのWindowsのインストール方法に飛ぶ。
New simplified installation procedure on Windows - Announcements - Graphviz
このサイトを読んでみると「管理者として実行せよ」との文字が。なるほど。コマンドプロンプトを管理者として実行してもう一度dot -cを打ち込んでみる。
Trueが返ってきた。どうやら管理者からしかアクセスできないようになっているらしい。
これでいけるやろと思いもう一度pythonからpngファイルの書き出し。ドンピシャ。無事決定木の画像データが表示された。
Graphvizは系統樹作成にも使えるから、graphvizをインストールできたのはでかい。
ちなみにdot言語を学ぶときは下のウェブサイトを参考にしよう
IGVの使い方
IGVの使い方について少し。
pythonのigv Jupyter Extensionを使えばどうやらpythonでも使えるらしいけど、リファレンスゲノムが多分サーバー上にあるやつしか使えないんじゃないかな。ということでサーバー上にない生物種を使うために結局アプリをインストール。
ちなみにigvでbamファイルを扱うためにはインデックス化してbam.baiの形にすることが必要です。
ここからはアプリでリファレンスゲノムを設定するための覚え書き。以下のサイト参照。リファレンスゲノムの読み込み | IGV ゲノムブラウザ (biopapyrus.jp)
要はリファレンスゲノムのFastaファイルを用意し、Genomeの欄からCreate.fileを選択してリファレンスゲノムをインポートすれば用意すればいいらしい。
OptionでGFFファイルを打ち込めるらしいが、必須ではないだろう。これは遺伝子名を入れたいときに必要なのか。
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にもパスは通してあるのでうまくできる。はず。これにこだわってもしょうがないので一回スルーしよう。
