RNA-seq による遺伝子発現解析

マッピング

シーケンシングされたリードをリファレンスゲノムにマッピングするとき、Bowtie、HISAT、STAR などのプログラムを利用する。リードをマッピングすることによって、どの転写領域に何リードマッピングされたかを計算することができ、間接的にその転写領域の発現量を測ることができる。

ヒト、マウス、シロイヌナズナなどのような全ゲノムが解読されているモデル生物のデータを処理する場合は、インターネットからリファレンスとなるゲノム(全ゲノムデータ)をダウンロードすればよい。一方、全ゲノムが解読されていない非モデル生物などの場合は、de novo アセンブリーにより、リファレンスゲノムを一度作成してから、解析に進む。両者の違いは、リファレンスとなるゲノムをインターネットからダウンロードするか、自作するかである。

モデル生物のリファレンスゲノムは Ensembl, RefSeq, UCSC などのデータベースからダウンロードすることができる。RefSeq ゲノムは NCBI から公開され、様々な研究に広く使われている。また、UCSC ゲノムは UniProt や GenBank mRNA などからデータを集めてアノテーションをつけている。Ensembl ゲノムは UCSC や Vega などデータを集めて自動アノテーションプログラムによってアノテーションを付けを行っている。アノテーション数としては Ensembl ゲノムが一番多い。三者において共通のアノテーションを持つ遺伝子が多く含まれているが、一方にあり他方にないような遺伝子も数多く含まれている。このようなアノテーションの違いが、トランスクリプトームの発現量を取得する際に大きく影響を与えてしまう。マッピング率として、アノテーション数の一番多い Ensembl の方が大きい。ただし、注意したいのはどのアノテーションを使おうと、不完全なアノテーションや完全に間違っているアノテーションが含まれていることを念頭におく必要がある。

Alignment-based transcript quantification

一般的なマッピングプログラムでは、リファレンス(ゲノム)配列からインデックスを作成して、インデックスを参照しながらリードをリファレンス(ゲノム)配列に効率よくマッピングしている。この種のツールは、マッピングだけで終わる。発現量を定量するには、マッピング結果を HTSeq や featureCountss に入力し、各転写領域にどれぐらいのリードがマッピングされたかを計数する必要がある。

Bowtie このバージョンの Bowtie は、マッピングするときにギャップを許容しない。リードが 50 bp 未満と比較的に短いときによく使われていた。
> 実行例: A. thaliana single-end
Bowtie2 短いリード(50-100 bp 前後) をマッピングするときによく使われていた。1 つのリードが 2 つのエクソンに跨る場合に、対応できない。
> 実行例: R. norvegicus single-end
> 実行例: H. sapiens paired-end
TopHat 1 つのリードが複数のエクソンに跨る場合でも、マッピングできる splice-aware aligner である。TopHat の内部では Bowtie2 を利用して、1 回目のマッピングを行う。続いて、1 回目のマッピングでマップされなかったリードを分断して、再度 Bowtie2 でマッピングを試みる。TopHat と TopHat2 の 2 つのバージョンが存在するが、とくに何かこだわりでもない限り最新版の TopHat2(あるいはその後継版である HISAT2)を使った方がいい。
> 実行例: D. melanogaster paired-end
> 実行例: A. thaliana single-end (TopHat2)
HISAT HISAT は TopHat を改良し、高速化したマッピングツールである。TopHat と同様に 1 つのリードが複数のエクソンに跨る場合にも対応できる splice-aware aligner である。HISAT と HISAT2 の 2 つのバージョンが存在する。HISAT2 は HISAT と TopHat2 の後継バージョンとしてリリースされている。HISAT と TopHat2 よりも HISAT2 を使用することが推奨されている。
BWA ギャップを許容してマッピングを行うことができる。1000 bp のような長いリードにも対応可能。塩基の挿入や欠損に強い。
> 実行例: H. sapiens single-end
STAR リード先端の塩基から、リファレンス配列へのマッピングを試みる。リードの途中までしかマッピングできなかった場合、中断された場所をリードの先端とみなして再度リファレンス配列へのマッピングを試みる。このように STAR も、 1 つのリードが複数のエクソンに跨る場合に対応している。TopHat に比べ非常に高速であるが、メモリ使用量も非常に大きい。

1 つの遺伝子から複数の転写産物が転写される。これら複数の転写産物は同じエクソンを持つ場合もある。遺伝子レベルの発現量解析では問題にならないが、転写産物レベルの発現量解析の場合にリードの分配方法が問題にある。例えば、ある遺伝子に 4 つのエクソン領域(E1、E2、E3、E4)があり、そこから 3 の転写産物(T1:E1-E2、 T2:E1-E2-E4、T3:E1-E3-E4)が作られるものとする。このとき、あるリードが E1 にマッピングされたとき、このリードはどの転写産物に分配されるかによって、あとあとの解析に影響を及ぼす。このようなリードの分配方法を考慮して、転写産物レベルの発現量解析を可能にしたプログラムも開発されてきている。

Cufflinks Bowtie/TopHat/HISAT などのマッピング結果である BAM ファイルを入力として、転写産物の構造を推定する。1 つのリードが同一遺伝子由来の複数転写産物にマッピングされる場合、推定された転写産物の構造に基づいて、この 1 リードを分割して、各転写産物に配分する。
RSEM Bowtie/TopHat/HISAT などのマッピング結果である BAM ファイルを入力として、転写産物の発現量推定を行う。1 つのリードに複数箇所のマップ候補がある場合、RSEM は、推定される転写産物の構造に基づいて、尤度などを計算して 1 つのリードをこれらの複数箇所に適切に配分する。

Alignment-free transcript quantification

リードをリファレンスゲノムにマッピングしてから発現量を定量する戦略の他に、マッピングせずに発現量を定量する戦略もある。後者の場合、マッピングのステップがないため、発現量の定量を非常に高速に行うことができる。

Sailfish Sailfish は、既知の転写産物の配列に対して k-mer のインデックスをあらかじめ作成しておく。次に、各 k-mer がリードに出現する頻度に基づいて、転写産物の定量を行う。
Salmon Salmon は、既知の転写産物の配列に対してインデックスを作り、インデックスを介してリードを転写産物配列へ擬似的にマッピングして、転写産物の定量を行う。Salmon は、Sailfish の後継バージョンとしての位置付けである。
> 実行例: A. thaliana paired-end
kallisto kallisto は、既知の転写産物の配列を k-mer に分け、de Bruijn グラフを作り、リードを de Bruijn グラフに当てはめて、転写産物の定量を行う。
> 実行例: A. thaliana paired-end

References

  • Zhao S, Zhang B. A comprehensive evaluation of ensembl, RefSeq, and UCSC annotations in the context of RNA-seq read mapping and gene quantification. BMC Genomics. 2015, 16:97. PubMed Abstract