kuroの覚え書き

96の個人的覚え書き

arabidopsisでhisat2

さて、これまでヒトゲノム研究でRNAseqパイプラインをせっせと作っていたのだが、縁あって植物科学の世界に舞い戻ってきた。

そこで植物のRNAseqを解析することになったのだが、やはりと言うか、色々と環境が整備されていなくて、いちいち対応していく必要がある。

前任者、というかコラボでやってもらっていたデータをはいっと渡されはしたものの、どういうストラテジーで解析されたファイルなのかわからないまま眺めてみても、やはり今ひとつ状況が把握できない。
ということで、まずはシングルエンドのシーケンスリードをtrimmomaticで処理して、Hisat2にかけてマッピングまでやってみることにする。

trimmomaticのシングルエンドリード処理はやったことがなかったのでパイプラインを見直して、まずは対応させるところから。
基本的にはSEというオプションを付けて、入力するファイル、出力するファイル名を1つづつにしてやればそれで良いようだ。

次に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

こういうふうに打つと、インデックスファイルが8分割して生成した。

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ファイルが4分割されているようなんで、これをどの段階でマージすればいいのかがわからず。
そのまま以前作ったパイプラインに流してみると、どうも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)