さて、BRAKER2による遺伝子予測が終わったら、今度は近縁種からの遺伝子アノテーションを移植・・・する前に、RNA-seq のデータからノイズを除去する操作を行う。

GeMoMa では最初に RNA-seq のマッピングデータに基づいてイントロンを予測するのだが、実は RNA-seq のデータというのはかなり厄介である。クオリティや厚みが足りないのは問題外だが、どんなにきれいな RNA を使って丁寧にライブラリを作っても、なんだこりゃ、というリードが一定量含まれてしまう。そもそも細胞内における RNA の転写というのがかなりいい加減なシステムだからだ。転写開始点じゃないところから転写が始まったり、逆に転写終了点を超えてリードスルーを起こしたり、スプライシングをミスったり、、、RNA-seq という技術は、そういう失敗作まで拾い上げてしまうから厄介なのだ。

下の図を例に取ると、２つのORFをまたいで４本のリードがマッピングされている。これ使ってそのまま遺伝子予測をすると、２つの遺伝子座がフュージョンした「それ絶対違うだろ」的な遺伝子座が予測されてしまう。人間の目にはこの４つのリードは無視すべきなのは明らかでも、これをうまいこと除去したり、無視して遺伝子予測を進めてくれるツールは現時点では存在しない。機械学習を使えばできるのかも知れないが、私も、そして恐らくこれを参考ししようとしている読者のみなさんも、自らの手でプログラムを作ることはできない人だろう。

だが、私やあなたのようなインフォマ素人でも、下図のようにノイジーなリードを除去することはできる。まずRNA-seq のデータをゲノムにマップする前に一旦アセンブルし、次にできた転写産物の配列にリードをマッピングし、そしてマップできたリードのみを回収するのだ。上の図に示されたような少数派のリードが、de novo assembly においては無視されることを利用するのである。

1. ツールのインストール

mamba install hisat2 samtools
mamba create -n trinity trinity

ここでは私が使い慣れているTRINITY (https://github.com/trinityrnaseq/trinityrnaseq) を使うが、この論文の検証結果を見る限り SPAdes (http://cab.spbu.ru/software/spades/) でもいいと思う。好きな方を使おう。
あと、hisat2 (http://daehwankimlab.github.io/hisat2/) の代わりに STAR (https://github.com/alexdobin/STAR) を使ってももちろん構わない。だが WSL 環境で STAR を動かすのはトリッキーなので hisat2 を標準としておく。

2. アセンブル

ショートリードを扱うときにデータの Pre-processing が推奨されていることが多いが、あまり気にしなくてよい。そもそもクオリティの低いリードが結果に貢献することはほとんどないから。

conda activate trinity
Trinity --seqType fq --left READ1.fastq.gz --right READ2.fastq.gz --max_memory 48G --CPU 24

オプションは色々ある（https://github.com/trinityrnaseq/trinityrnaseq/wiki/Running-Trinity) が、特に使う必要もないと思う。メモリは 100万リードペアあたり 1 GB ほどあればよいらしい。Trinity.fasta がアウトプット。途中で止まった場合はほぼメモリ不足なのでメモリを増やすかCPUを減らすかしてやり直そう。

3. リードマッピング

TRINITYでできたTrinity.fastaに元のリードをマッピングする。trinity 環境を activate した状態のままだったら、まずは conda deactivate で base 環境に戻ろう。マッピングの実行手順は以下の通り。

hisat2 用 index ファイルの作成

hisat2-build -p 8 Trinity.fasta INDEX

-p で複数のCPUを使った並列処理ができるが、あまり欲張るとメモリ不足を起こして失敗するので注意。しかもINDEXファイルがちゃんとできてしまっているように見えるので、気付かずに次のコマンドを打って永遠に何も起きない（走っているフリをされる）という罠に陥る。せいぜい一桁にしておこう。

リードマッピングの実施

hisat2 -x INDEX -1 READ1.fastq.gz -2 READ2.fastq.gz -p 24 | samtools view -f 0x2 -bh | samtools sort -@24 -o Trinity.bam && samtools index -@24 Trinity.bam

あるいは

hisat2 -x INDEX -1 READ1.fastq.gz -2 READ2.fastq.gz -p 24 -S Trinity.sam
samtools view -f 0x2 -bh Trinity.sam > temp.bam
samtools sort -@24 -o Trinity.bam temp.bam
samtools index -@24 Trinity.bam

-x で INDEX ファイルのprefixを指定、-1 で Forward Read のファイルを、-2 で Reverse Read のファイルを指定する。
-p は CPU の数、-S はアウトプットとなる.sam ファイルの名前の指定だ。 -S を付けなければ「標準出力（知らない人はググろう）」になるので、パイプ ( | ) を使って samtools に流し込むことができる。
samtools view はマッピングクオリティ等でフィルタリングを掛けるときに使うコマンドで、-f 0x2 でペアエンドの両方のリードが正しい向きにマッピングされたリードだけを抽出する。-q でマッピングクオリティによるフィルタリングもできるが、あるリードが複数の転写産物にマッピングされることはあり得るのでここでは使わない。
samtools sort はリードのポジションソートに使う。IGVでマッピングの結果をブラウジングするときには染色体やコンティグの配列に沿ってアラインメントを並び替えておかなければいけないので、リードマッピングから BAM ファイルを作るときにはほぼ自動的にで使う。なお、samtools ではCPUを複数使って並列化するときには -@ を使う。
samtools index で .bam ファイルをインデックス化して .bai ファイルを作る。IGV でのブラウジングには必須だし、これがあると色んな処理が高速化するので、これもほぼ自動的に使う。

さすがに.sam や .bam についての説明は省略するが、知らないと苦労するので知らない人は自習すること。

3. BAM ファイルからのリードの回収

samtools bam2fq -n -1 CLEAN_1.fastq.gz -2 CLEAN_2.fastq.gz Trinity.bam

-n は回収したリードのリード名に変更を加えないためのオプション。忘れると Forward Read のリード名の最後に "_1" が、 Reverse Read のリード名の最後に "_2" が追加されてしまい、ツールによっては不都合になることも。

4. 回収したリードをゲノムにマッピング

hisat2-build -p 8 GENOME.fasta GENOME
hisat2 -x GENOME -1 CLEAN_1.fastq.gz -2 CLEAN_2.fasta.gz -p 24 --max-intronlen 15000 | samtools sort -@24 -o CLEAN.bam && samtools index -@24 CLEAN.bam

--max-intronlen で最大イントロン長（gapped alignment で許すギャップの大きさ) を指定。植物のイントロンはかなり小さいので15 kbくらいにしておけばよい。動物だと数百 kbp のイントロンもザラにあるので、妥当な数値を自分で調べて適宜調整。
これで GeMoMa を実行する準備が整った。

関連

ゲノムアノテーション（目次）
RepeatModeler/RepeatMaskerを使ってゲノムにSoftmaskを掛ける
BRAKER2のインストール
BRAKER2を使ってタンパク質コード領域を予測する
GeMoMa で近縁種のアノテーションを移植する
Gffcompareを使ってGeMoMaとBrakerの差分を取り、アノテーションを補完する
EnTAP のセットアップ
EnTAP の実行と結果の解釈