さて、これまでヒトゲノム研究でRNAseqパイプラインをせっせと作っていたのだが、縁あって植物科学の世界に舞い戻ってきた。
そこで植物のRNAseqを解析することになったのだが、やはりと言うか、色々と環境が整備されていなくて、いちいち対応していく必要がある。
前任者、というかコラボでやってもらっていたデータをはいっと渡されはしたものの、どういうストラテジーで解析されたファイルなのかわからないまま眺めてみても、やはり今ひとつ状況が把握できない。
ということで、まずはシングルエンドのシーケンスリードをtrimmomaticで処理して、Hisat2にかけてマッピングまでやってみることにする。
trimmomaticのシングルエンドリード処理はやったことがなかったのでパイプラインを見直して、まずは対応させるところから。
基本的にはSEというオプションを付けて、入力するファイル、出力するファイル名を１つづつにしてやればそれで良いようだ。
次にHisat2でのマッピング。
そもそもリファレンスファイルから準備しなくてはならない。
http://plants.ensembl.org/info/website/ftp/index.html
こちらのサイトに植物のゲノム情報が集約されているようなので、そこから該当するファイルを貰ってくることにした。
とりあえずはarabiのDNAとアノテーションgtfファイルをダウンロードする。
更にHisatでマッピングするにはインデックスファイルが必要となるのだが、Hisat2のページで配布されているのは動物関係しかないようなので、これも自前で作成してやらねばならない。
ダウンロードしたfaファイルの階層にて
$ hisat2-build Arabidopsis_thaliana.TAIR10.dna.toplevel.fa Arabidopsis_thaliana.TAIR10.dna.toplevel
こういうふうに打つと、インデックスファイルが８分割して生成した。
Headers:
len: 119481543
gbwtLen: 119481544
nodes: 119481544
sz: 29870386
gbwtSz: 29870387
lineRate: 6
offRate: 4
offMask: 0xfffffff0
ftabChars: 10
eftabLen: 0
eftabSz: 0
ftabLen: 1048577
ftabSz: 4194308
offsLen: 7467597
offsSz: 29870388
lineSz: 64
sideSz: 64
sideGbwtSz: 48
sideGbwtLen: 192
numSides: 622300
numLines: 622300
gbwtTotLen: 39827200
gbwtTotSz: 39827200
reverse: 0
linearFM: Yes
Total time for call to driver() for forward index: 00:02:07まあこんなもんでしょう。
さてさて、次にHisat2にかけてみよう。
そもそもfastqファイルが４分割されているようなんで、これをどの段階でマージすればいいのかがわからず。
そのまま以前作ったパイプラインに流してみると、どうもBAMファイルに変換して、ソートをかけようとするとエラーが出て先に進まないようなので、ソートする前、もしくはBAMに変換する前にどうにかしないといけないのかもしれない。
やってみたところ、BAMに変換後、samtools mergeで結合したBAMファイルに対してsort-indexとかければちゃんと動作し、出来上がったBAMファイルをIGVブラウザで閲覧することができた。
つぎに、これをcufflinksにかけて発現量を計算してみる。
$ cufflinks -p 4 -q -g "$gtf" "$output_dir"/"$exp".bam -o "$output_dir2"
gtfファイルやインプット・アウトプットは前もって指定してあるのでこんな感じ。4コアなiMacでうりゃっとかけると
見事に4コア使い切っていて清々しい。1サンプル40分ってところだね。(2.7 GHz Intel Core i5 4core）
