RNA-Seqでは、readを参照ゲノムや転写産物配列に対応づけ、 遺伝子ごとの発現量を計算します。
しかし、配列のよく似た遺伝子ファミリー、pseudogeneを持つ遺伝子、 互いに重なる遺伝子などでは、 同じreadが複数の遺伝子に対応することがあります。
ケーススタディ No.414 「RNA-Seqはマイクロアレイの単純な上位互換ではない」 では、次のように説明しました。
相同性の高い領域を共有する遺伝子ファミリーでは、 short-read RNA-Seqのreadが複数の遺伝子に対応することがあります。 解析方法によっては、multi-mapping readsを除外する場合もあれば、 複数の候補遺伝子に確率的に分配して発現量を推定する場合もあります。
このようなreadが多い遺伝子では、 個々の遺伝子へのreadの割り当てが解析手法やアノテーションに依存し、 本来の遺伝子ごとの発現量や発現差を正確に推定できない場合があります。 つまり、配列類似性の高い遺伝子群では、 RNA-Seqの測定値が、実際の遺伝子ごとの発現差を十分に反映しない場合があります。
この問題がどのように検証されているのかを確認するため、 関連する研究についてAIに質問し、紹介された文献を調べました。
曖昧なreadはどのように処理されるのか
複数の遺伝子に対応するreadの扱いは、 大きく次のように分けられます。
| 処理方法 | 扱い | 注意点 |
|---|---|---|
| 除外する | 複数の遺伝子に対応するreadを数えない | 一意に対応する配列が少ない遺伝子ほど、発現量が低く見えやすい |
| すべての候補に数える | 対応するすべての遺伝子にreadを加える | 同じreadが重複して数えられ、発現量が水増しされる |
| 均等に分配する | 候補遺伝子へ同じ割合で配分する | 実際の発現比率が異なっていても均等に分けられる |
| モデルで分配する | unique readsや推定発現量などを使って確率的に配分する | 区別する情報が乏しい場合、結果がモデルやアノテーションに依存する |
処理方法によって結果は変わる
CoCoの研究では、 nested genesを対象に、 HTSeq-count、featureCounts、RSEMなどの定量結果が比較されています。
複数の遺伝子に割り当てられるreadを除外する標準的な設定では、 実際にはreadが存在する遺伝子が検出されない、 または過小評価される場合がありました。 一方、RSEMや設定を変更したfeatureCountsでは改善が見られましたが、 readが正しい比率で各遺伝子に割り当てられるとは限りませんでした。
single-cell RNA-Seqの定量法である alevinの研究でも、 multi-mapping readsを捨てると、 影響を受ける遺伝子のcountsが系統的に低くなることが示されています。
つまり、曖昧なreadを除外すれば発現量を失う可能性があり、 均等に配分すれば本来の発現比率を崩す可能性があります。
モデルで分配すれば正しくなるのか
Salmonや RSEMは、 一意に対応するread、推定発現量、fragment lengthなどを利用して、 曖昧なreadを候補となる転写産物へ確率的に配分します。
遺伝子を区別するreadが十分にあれば、 単純な除外や均等分配より妥当な推定ができる可能性があります。
しかし、二つの遺伝子を区別できる配列がほとんど存在しなければ、 RNA-Seqデータだけから本当の発現比率を決めることはできません。 それでも解析ソフトウェアは、モデルに基づいた発現量を出力します。
曖昧なreadを統計モデルで分配する方法であっても、 使用するプログラムによって同じ結果になるとは限りません。 readと候補転写産物の対応づけ方、 発現量の推定アルゴリズム、 bias補正、アノテーションの扱いなどが異なるためです。
配列類似性の高い多重遺伝子ファミリーをシミュレーションデータで比較した研究 では、 反復配列の多いTrypanosoma cruziを対象として、 使用する定量方法によって、 個々の遺伝子へのreadの割り当てや推定誤差が大きく異なり、 配列類似性が100％に近づくほど 正確な定量が難しくなることが示されています。
遺伝子を区別する情報がRNA-Seqデータに十分含まれていない場合、 出力される発現量は、 使用するプログラムやアノテーションによって大きく異なる可能性があります。 また、いずれかの推定値が、 実際の遺伝子ごとの発現比率を正確に表しているとは限りません。
この点は、1細胞あたりに得られるread数が少なく、 3′末端または5′末端の限られた領域を測定することが多い single-cell RNA-Seqでは、さらに重要になります。
遺伝子を一意に区別できるreadが個々の細胞で得られない場合、 解析手法によっては、readやUMIがどの候補遺伝子に対応するかという情報を利用して、 曖昧なreadの由来を推定します。 alevinの研究では、 遺伝子間で曖昧なreadを除外せずに利用することで、 gene-level countsの過小評価を抑える方法が検討されています。
このような推定はcountsを改善する可能性がありますが、 個々の細胞でその遺伝子を固有に識別して測定したことを 意味するわけではありません。 相同性の高い遺伝子を細胞型のmarkerとして使用する場合には、 その発現値がどのようなreadの割り当て方法によって 計算されたものかにも注意が必要です。
発現量だけでなく、遺伝子間の関係にも影響する
遺伝子間のcross-mappabilityを調べた研究では、 配列類似性によるreadの誤割り当てが、 見かけ上の遺伝子間の共発現や、 偽陽性のtrans-eQTLシグナルを生じさせる可能性が示されています。
multi-mappingの影響は、 個々の遺伝子の発現量だけでなく、 相関解析やネットワーク解析にも及ぶ可能性があります。
通常のbulk RNA-Seqではどの程度問題になるのか
150 bp paired-endのbulk RNA-Seqでは、 短いsingle-end readより多くの配列情報を利用できるため、 readの由来を区別しやすくなります。 ただし、配列類似性の高い遺伝子では、 paired-endであっても曖昧さが残る場合があります。
また、複数のtranscriptに対応するreadであっても、 それらが同じgeneに属するなら、 gene-level解析では大きな問題にならない場合があります。
注意が必要なのは、同じfragmentが 異なる遺伝子に対応する場合です。 特に、paralog、pseudogeneを持つ遺伝子、 重複領域にある遺伝子では、 gene-levelの発現量でも解析方法の影響を受ける可能性があります。
3′ RNA-Seqやマイクロアレイなら区別できる可能性もある
遺伝子全体では相同性が高くても、 3′側にはその遺伝子に固有の配列が存在する場合があります。
その場合、3′末端付近を測定する3′ RNA-Seqや、 固有配列にprobeが設計されたマイクロアレイでは、 通常のRNA-Seqより遺伝子を区別しやすい可能性があります。 3′ UTRは配列の違いが比較的大きいことがあり、 マイクロアレイのprobe設計に利用されてきました。
ただし、3′側にも相同性があれば改善しません。 マイクロアレイでも、probeが類似配列に反応すれば cross-hybridizationが起こります。
本来であれば、定量が難しい遺伝子について、 3′末端配列や実際のmicroarray probeまで確認し、 測定方法ごとの識別可能性を比較する必要があります。 しかし、今回の調査では、 そこまでを包括的に検証した研究は確認できませんでした。
計算結果をそのまま信じない
今回確認した文献から、 相同性の高い遺伝子では、 multi-mapping readsの扱い方によって 発現量の推定結果が変わり得ることが分かりました。
解析ソフトウェアから数値が出力されたとしても、 その数値が遺伝子ごとの実際の発現量を 正確に反映しているとは限りません。
特に、配列のよく似た遺伝子を比較する場合には、 使用した定量方法、multi-mapping readsの扱い、 一意に対応する配列の有無を確認しながら 結果を解釈する必要があります。
RNA-Seqで重要な候補遺伝子が見つかった場合には、 BAMファイルを IGVなどで確認し、 readがどのexonやsplice junctionに対応しているかを 調べることも有効です。
より詳しく確認する場合には、 STARや RSEMから transcript座標のBAMファイルを出力し、 exonが連結された成熟mRNA配列上で readの対応位置を確認する方法もあります。
さらに、readが集中して見える領域や、 サンプル群間でcoverageの差が大きい領域の配列を取得して NCBI BLASTで検索し、 他の遺伝子やpseudogeneにも高い配列類似性を持つ領域がないかを 確認できます。
特に、発現量や発現差が研究の結論に大きく関わる遺伝子については、 readが対応する領域にその遺伝子固有の配列が含まれているかを確認し、 必要に応じて固有領域を対象としたqPCRなどで 検証することが望まれます。
