bowtie2によるローカルアライメントを使い、Poly-Tプライマーで作成したcDNAライブラリのRNA-seq解析を行う。
まずはbowtie2をインストールする。
$ git clone https://github.com/BenLangmead/bowtie2.git $ cd bowtie2/ $ make $ sudo make install
次にリファレンスのインデックスを作成する。
$ cd PATH/TO/reference $ mkdir target_ref_bowtie2index $ bowtie2-build --threads 40 ./genome_assembly.fasta ./target_ref_bowtie2index/genome_assembly
ゲノムのサイズにもよるがマルチスレッドを使ってもかなりの時間を要する。
ローカルアライメントオプションを入れてマッピングしてみる
bowtie2 -p 40 -x [リファレンスゲノムのインデックスへのパス] -U [RNA-seqデータのファイル.fastq.gz] --local -k 1 -S [出力SAMファイル名.sam] 2> mapping_report.txt
標準ではエラー出力にしかレポートが出力されないため、バッチ処理等で大量にまとめて処理するとレポートが取れない。なので2>でファイルにリダイレクトしている。
また-k 1オプションにより、1つのリードが複数箇所にマッピングされるときは最良の1箇所だけを残す選択をしている。（これは必須ではない）
このあと生成したSAMをsamtoolsでBAMに変換、ソート、インデックス作成までは問題なく実施できた。
これをRSEMで発現解析させたいのだが、まだ上手く行っていない。
追記
どうもうまく行かない。
結局bowtie2で前もってゲノムにマッピングしておいて、それをRSEMで作成したtranscriptデータベースを使って発現解析をするという手順に問題がありそうだ。最初からbowtie2のtranscriptインデックスを作成してそちらにマッピング〜発現解析までをRSEMで一括して行うのが良さそう。
まずRSEMでtranscriptデータベースを作成する。
$ rsem-prepare-reference -p 40 --bowtie2 --bowtie2-path /usr/local/bin --gff3 genes.gff3 genome.fa rsem_db/transcript_db
ポイントとしては--bowtie2-pathはあくまでbowtie2のあるディレクトリのパスであり、bowtie2まで入れないこと。このコマンドでRSEM用の発現解析用データベースと一緒にbowtie2のマッピング用インデックスも生成される。
そしてbowtie2でのマッピングもRSEMの中でやってしまう。
$ rem-calculate-expression -p 40 --bowtie2 --bowtie2-path /usr/local/bin input_reads_1.fastq.gz --bowtie2-options "--local" rsem_db/transcript_db output_prefix
と思ったが、RSEMでは--localオプションを扱えないらしい。上のようなオプションの付け方は無効だ。そのままend-to-endで解析すると結果としてマップ率が非常に低くなってしまうため、トリミングの時点で3’末端を削っておく必要がある。(以下の例では50base削った）
fastp --trim_tail1 50 -w 40 -i input_reads_1.fastq.gz -o input_reads_1_trim50.fastq.gz -q 20 -l 20 -j input_reads_1_trim_report50.json -h input_reads_1_trim_report50.html
そのうえで、
$ rsem-calculate-expression -p 40 --bowtie2 --bowtie2-path /usr/local/bin input_reads_1_trim50.fastq.gz rsem_db/transcript_db output_prefix 2> mapping_report.txt
レポートのリダイレクトも有効だ。
マッピングレートは--localオプションを付けると90%以上であったが、トリミングしてのend-to-endアライメントではせいぜい60％くらいにしかならなかった。
さて、これでとりあえずリードカウントができることがわかったが、少しまだ問題が残る。
上記パイプラインをコマンドラインで順番に処理する上では問題がないのだが、クラスタサーバでslurmを使ったジョブスケジューリングで自動運転した場合、RSEM内部で呼び出しているsamtoolsがうまく見つけてもらえない。bowtie2は明示的にパスを与えているのだが、samtoolsにはそのオプションがなく、処理が止まってしまうようだ。
なので、やはりbowtie2によるマッピングはRSEM内部で処理せず、個別に実施した上で、RSEMではリードカウントのみ実施してもらうのが良さそうだ。
そこで、
$ bowtie2 -p 40 -x rsem_db/transcript_db -U input_reads_1_trim50.fastq.gz -S output_prefix.sam 2> output_prefix_mapping_report.txt $ samtools view -@ 40 -bS output_prefix.sam > output_prefix.bam $ samtools sort-@ 40 output_prefix.bam -o output_prefix_sort.bam $ samtools index output_prefix_sort.bam $ rsem-calculate-expression -p 40 --alignments output_prefix_sort.bam rsem_db/transcript_db output_prefix
というパイプラインにしてみたのだが、
どういうわけか
rsem-parse-alignments rsem_db/transcript_db output_prefix.temp/output_prefix output_prefix.stat/output_prefix output_prefix.transcript_sort.bam 1 -tag XM Read LH00220:78:22H5YFLT3:7:1284:39149:5518: RSEM currently does not support gapped alignments, sorry! "rsem-parse-alignments rsem_db/transcript_db output_prefix.temp/output_prefix output_prefix.stat/output_prefix output_prefix.transcript_sort.bam 1 -tag XM" failed! Plase check if you provide correct parameters/options for the pipeline!
というエラーが出る。どうもbowtie2は厳密にギャップアライメントを排除することが難しく、素のbowtie2のデフォルトオプションとRSEM内部でやっているオプションが違うらしい。
RSEMのヘルプの‐‐bowtie2オプションの説明にヒントが有った。
--bowtie2
Use Bowtie 2 instead of Bowtie to align reads. Since currently RSEM
does not handle indel, local and discordant alignments, the Bowtie2
parameters are set in a way to avoid those alignments. In
particular, we use options '--sensitive --dpad 0 --gbar 99999999
--mp 1,1 --np 1 --score-min L,0,-0.1' by default. The last parameter
of '--score-min', '-0.1', is the negative of maximum mismatch rate.
This rate can be set by option '--bowtie2-mismatch-rate'. If reads
are paired-end, we additionally use options '--no-mixed' and
'--no-discordant'. (Default: off)つまり
--sensitive --dpad 0 --gbar 99999999 --mp 1,1 --np 1 --score-min L,0,-0.1
というオプションがデフォルトで入っているらしい。
なので上のパイプラインの1行目を
$ bowtie2 --sensitive --dpad 0 --gbar 99999999 --mp 1,1 --np 1 --score-min L,0,-0.1 -p 40 -x rsem_db/transcript_db -U input_reads_1_trim50.fastq.gz -S output_prefix.sam 2> output_prefix_mapping_report.txt
としてやることでRSEM内部での動作と同等にしてやることが可能となった。
今後ローカルアライメントを受け付けるリードカウントについては検討の余地があると思う。そもそもbowtie2を使いたかったのもローカルアライメントが使えるからであり、それが使えないならSTARでマッピングしても同じようなものだから。bowtie2を内部で実施しないなら、そもそもRSEMを使う必然性も低く、これにかわるリードカウント方法を考えてもいいように思う。RSEMのカウントが比較的正確だ、というのもあるけど。
