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

マッピング

シークエンシングされたリードは、トランスクリプトの一部となっている。これらのリードの塩基配列とゲノム(リファレンス)配列を比べることで、各々のリードがゲノム上のどの遺伝子領域に由来するのかを明らかにすることができる。このように、リードとゲノム配列を比較して、リードの由来箇所を決めることをマッピングあるいはアラインメントとよぶ。

RNASeq マッピングはリードとゲノム配列を比較して、リードの由来箇所を決めることである。

マッピング用のプログラムとして、STARHISAT2BWA などが有名である。これらのプログラムは、FM index とよばれる文字列検索アルゴリズムを使用してマッピングを行なっている。したがって、これらのプログラムを使用するには、まずゲノム配列を FM index で検索できるように Burrows-Wheeler 変換しておく必要がある。この作業をインデックスの作成などという。インデックスを作成してから、FM index でリードと完全一致するゲノム配列上の位置を検索する。FM index では完全一致による検索アルゴリズムであるため、ギャップやミスマッチを許容できない。そこで、これらのプログラムは、まずリードの一部分だけをゲノム配列に対して完全一致で検索し、ヒット箇所があればその近くで Smith-Waterman アルゴリズムでアラインメントを伸ばすこととしている。

マッピングの際に使用するリファレンス配列は、モデル生物であれば、データベースで公開されていることが多い。有名なデータベースとして、RefSeqUCSCEnsembl などのデータベースがある。RefSeq ゲノムは NCBI から公開され、様々な研究に広く使われている。また、UCSC ゲノムは UniProt や GenBank mRNA などからデータを集めてアノテーションをつけている。Ensembl ゲノムは UCSC や Vega などデータを集めて自動アノテーションプログラムによってアノテーションを付けを行っている。アノテーション数としては Ensembl ゲノムが一番多い。三者において共通のアノテーションを持つ遺伝子が多く含まれているが、一方にあり他方にないような遺伝子も数多く含まれている。このようなアノテーションの違いが、トランスクリプトームの発現量を取得する際に大きく影響を与えてしまう(Zhao et al, 2015)。マッピング率として、アノテーション数の一番多い Ensembl の方が大きい。ただし、注意したいのはどのアノテーションを使おうと、不完全なアノテーションや完全に間違っているアノテーションが含まれていることを念頭におく必要がある。非モデル生物のシークエンサーデータを解析する場合は、リファレンス配列を作ることから始める必要がある。シークエンサーデータ、すなわちリードの配列を利用して、元の転写領域の配列を復元することをアセンブリという。非モデル生物を対象する場合に、アセンブリを行なってからマッピングを行うが一般的である。

シークエンサーのデータ量が増えるにしたがって、マッピングに要する時間も長くなりつつある。そこで、マッピングを行うのではなく、リファレンス配列をハッシュ化して、リードの配列をハッシュのキーとして、そのキーがどの転写産物の配列から作られたハッシュに存在するのかを調べることで、擬似的にそのリードの由来元を推定する プログラムが開発されている。kallisto や salmon などがハッシュを利用した擬似アラインメントプログラムである。この種のプログラムは、速度が速いものの、ミスマッチや indel に対して弱いことが知られている。解析対象の生物種のゲノムの質とそのアノテーションの質に応じて使うかどうかを決めると良い。

Alignment-based transcript quantification

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

Bowtie このバージョンの Bowtie は、マッピングするときにギャップを許容しない。リードが 50 bp 未満と比較的に短いときによく使われていた。
> 実行例: A. thaliana single-end RNA-Seq
Bowtie2 短いリード(50-100 bp 前後) をマッピングするときによく使われていた。1 つのリードが 2 つのエクソンに跨る場合に、対応できない。
> 実行例: A. thaliana single-end RNA-Seq
> 実行例: A. thaliana paired-end RNA-Seq
TopHat2 1 つのリードが複数のエクソンに跨る場合でも、マッピングできる splice-aware aligner である。TopHat2 の内部では Bowtie2 を利用して、1 回目のマッピングを行う。続いて、1 回目のマッピングでマップされなかったリードを分断して、再度 Bowtie2 でマッピングを試みる。現在では TopHat2 をダウンロードして使うことは可能だが、とくに何かこだわりでもない限りその後継版である HISAT2 を使用することが推奨されている。
> 実行例: A. thaliana single-end RNA-Seq
> 実行例: A. thaliana paired-end RNA-Seq
HISAT2 HISAT は TopHat を改良し、高速化したマッピングツールである。TopHat と同様に 1 つのリードが複数のエクソンに跨る場合にも対応できる splice-aware aligner である。HISAT と HISAT2 の 2 つのバージョンが存在する。HISAT2 は HISAT と TopHat2 の後継バージョンとしてリリースされている。HISAT と TopHat2 よりも HISAT2 を使用することが推奨されている。
> 実行例: パン小麦(T. aestivum) paired-end RNA-Seq
> 実行例: A. thaliana paired-end RNA-Seq
> 実行例: A. thaliana single-end RNA-Seq
BWA PacBio などが出力 DNA-Seq ロングリードのマッピングによく使われている。ギャップを許容してマッピングを行うことができ、塩基の挿入や欠損に強い。1000 bp のような長いリードにも対応可能。ただし、塩基の挿入には対応しているものの、イントロンのような長い挿入がある場合に失敗する。そのため、BWA はイントロンを持つ真核生物の RNA-Seq マッピングに使われない傾向がある。
> 実行例: トウモロコシ(Z. mays) paired-end DNA-Seq
> 実行例: トウモロコシ(Z. mays) single-end DNA-Seq
STAR リード先端の塩基から、リファレンス配列へのマッピングを試みる。リードの途中までしかマッピングできなかった場合、中断された場所をリードの先端とみなして再度リファレンス配列へのマッピングを試みる。このように STAR も、 1 つのリードが複数のエクソンに跨る場合に対応している。TopHat に比べ非常に高速であるが、メモリ使用量も非常に大きい。
> 実行例: A. thaliana single-end RNA-Seq
> 実行例: A. thaliana paired-end RNA-Seq

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 リードを分割して、各転写産物に配分する。
StringTie HISAT などのマッピング結果である BAM ファイルを入力として、転写産物の構造を推定する。Cufflinks の後継プログラムとして使われる。
RSEM Bowtie/Bowtie2/STAR を内部的に呼び出してマッピングを行い、マッピング結果を解析して、転写産物の発現量推定を行う。1 つのリードに複数箇所のマップ候補がある場合、RSEM は、推定される転写産物の構造に基づいて、尤度などを計算して 1 つのリードをこれらの複数箇所に適切に配分する。

Alignment-free transcript quantification

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

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

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