RNA-Seqのdeconvolution解析には、 xCell、CIBERSORTx、MCP-counter、EPIC、quanTIseqなど、 複数の解析ツールがあります。
しかし、実際に解析を始めようとすると、 まずどのツールを選ぶべきか迷います。 判断基準としては、
- インプットする値が何か
- ノイズやバッチエフェクトなどの処理は含まれるか
- ペア構造を統計解析に使えるか
- どのような細胞種の推定ができるか
- 結果として得られる値はどのように使えるか
- どのような図を出力できるか
といったことになるでしょう。
今回は、アトピー性皮膚炎患者のペアRNA-Seqデータ GSE121212を題材として、 これらをAIと相談しながら一つずつ決めました。 解析ツールの比較から始め、 入力データの前処理を設計し、 その条件に合わせたRスクリプトをAIに生成してもらいました。 さらに、実行時のエラー修正、 統計結果の整理、 生物学的解釈、 作図方法についてもAIと対話しながら調整しました。 この一連の作業を、 従来では考えられないほど短時間で、 解析方針の決定から結果の可視化まで進めることができました。
RNA-Seq deconvolution解析ツールを比較する
RNA-Seq deconvolution解析には複数のツールがありますが、 すべてが同じ結果を出すわけではありません。
推定対象、出力値、必要な入力データ、 実行方法が異なるため、 解析目的に合ったツールを選ぶ必要があります。
| ツール | 主な対象 | 出力 | 実行環境 | 主な入力 | 長所 | 注意点 |
|---|---|---|---|---|---|---|
| xCell | 免疫細胞・間質細胞など64細胞種 | Enrichment score | R、Web | Gene Symbol × Sampleの発現行列。RNA-SeqではTPMなどの遺伝子長・ライブラリーサイズ補正値が推奨されるが、今回は解析目的に応じて正規化Gene Countsを使用 | 細胞種が広く、免疫細胞以外も含む。各サンプルの相対的な変化を探索しやすい | 細胞割合ではない。近縁細胞間のsignatureの重複や、低発現値の扱いに注意が必要 |
| CIBERSORTx | 主に免疫細胞。カスタムsignatureも作成可能 | 相対的または絶対的な細胞量推定、細胞種別発現推定 | Web、Docker | Mixture expression matrixとsignature matrix | 広く利用されており、細胞割合に近い出力を得られる。scRNA-Seqから独自signatureを作成できる | 既定のLM22は主に白血球を対象とする。referenceと解析組織の違いが結果に影響する |
| MCP-counter | 免疫細胞8群、endothelial cells、fibroblasts | Abundance score | R | Gene Symbol × Sampleの発現行列 | 対象となる細胞集団が限定されている分、主要な免疫・間質細胞のサンプル間比較に使いやすい | 細胞割合ではない。対象となる細胞集団が限定され、細かな免疫細胞サブタイプまでは区別できない |
| EPIC | 免疫細胞、fibroblasts、endothelial cells、その他の細胞 | 細胞割合またはmRNA fraction | R、Web | Bulk RNA-Seq発現行列 | 免疫細胞以外を含めた構成比を推定できる | 主に腫瘍微小環境を想定して設計されている。皮膚など別組織への適用ではreferenceとの違いに注意が必要 |
| quanTIseq | 10種類の免疫細胞 | 細胞割合と未分類細胞の割合 | R、コマンドライン、Web | FASTQまたはRNA-Seq発現行列 | FASTQから解析可能で、細胞割合として解釈しやすい | 主に血液・腫瘍免疫を対象とし、非免疫細胞の分類は限定的 |
| immunedeconv | 複数のdeconvolution手法を統一して実行 | 使用する手法に依存 | R | 通常はGene Symbol × Sampleの発現行列 | xCell、MCP-counter、EPIC、quanTIseqなどを同じ形式で比較できる | ツールそのものではなく、複数手法へのインターフェース。各手法の出力の違いを理解する必要がある |
今回は、皮膚組織に含まれる免疫細胞と間質細胞について、 まず広く変化の候補を探索し、 得られたスコアを患者内で比較することを目的としました。 そこで、対象細胞種が広く、 Rを使ってローカル環境で実行できることから、 xCellを選びました。 ただし、xCellの出力は細胞割合ではなく、 細胞種に関連する遺伝子群のenrichment scoreです。 この点は、統計結果の解釈と作図方法を決めるうえで注意が必要です。
xCellに何を入力するかをAIと検討する
ツールを選んだ後は、 xCellにどの発現値を入力するかを検討しました。
まずAIは、xCellではRNA-Seq入力として標準的に使われる TPMを提案しました。 しかし、TPMからは元のread数が分からないため、 低Counts領域の測定精度を評価できません。
そこで、AIと相談しながら、遺伝子長補正よりも、 低Counts領域に生じる不安定なダイナミックレンジの影響を確認することを優先し、 正規化Gene Countsから二つの入力を作成することにしました。
Input 1は、Gene Countsに正規化とLog2変換を適用したものです。 対数変換により0は欠損値となります。 正規化とLog2変換を施したGene Counts表の最小値が-0.6436だったので、 Counts 0ともともとの欠損値だったセルを-0.7で埋めました。
ところで、正規化後の低Counts側の下限値は、 サンプルごとのダイナミックレンジによって異なります。 そのため、Input 1には、生物学的な違いだけでなく、 低Counts領域の測定の不安定さや、 サンプルごとのダイナミックレンジの違いを反映した 非生物学的な変動も含まれています。
Input 2では、この問題に対処するため、 低Counts領域の不安定さと非生物学由来の変動を xCellの計算へそのまま反映させないことにしました。
Low Signal Cutoffを4.3219（線形値で20）に設定し、 1以上20未満に相当する値を4.3219へ置き換えました。 また、元からの欠損値と、 0 CountsがLog2変換によって欠損値となった部分は4で補完しました。 ここで示している閾値は、散布図を見ながら決めたもので、 どのようなデータセットに無条件に適用可能なものではありません。 これにより、 低Counts領域に含まれる測定の不安定さや、 サンプルごとのダイナミックレンジの違いによって生じる 非生物学由来の変動を解析へ反映させず、 0または欠損値と、 1以上のCountsが得られた低発現値は区別しました。
今回使用した正規化Gene Countsのインプットファイルはこちらからダウンロードできます。
解析条件に合わせたRスクリプトをAIに生成してもらう
解析方針が決まった後、 入力ファイルの形式と実行したい処理をAIへ伝え、 Rスクリプトを生成してもらいました。
生成したRスクリプトでは、次の処理を自動化しました。
- メタデータと発現値の読み込み
- Entrez Gene IDからGene Symbolへの変換
- 重複Gene Symbolの整理
- Input 1とInput 2に対するxCell解析
- ペアを無視したWelchのt検定
- 患者内のpaired t-test
- 多重検定補正
- 患者ごとの増加・減少方向の集計
- Input 1とInput 2の結果比較
最初の実行では、エラーが発生しました。 エラーメッセージをAIへ提示すると修正案が示され、 解析を継続できました。
このようにAIは、コードを一度生成して終わるものではなく、 実行結果やエラーを確認しながら、 スクリプトを修正する相手として利用できます。
統計結果をAIと一緒に確認する
解析後、Input 1とInput 2の統計結果ファイルをAIへ渡し、 主要な結果が前処理によって変わるかを確認しました。 paired p-valueが0.05未満となった項目は、 Input 1で28項目、Input 2で27項目であり、 27項目は両方に共通していました。 Th2 cellsやImmuneScoreなどの主要な結果は、 Low Signal Cutoffの有無にかかわらずほぼ再現しました。
一方、B-cellsのpaired p-valueは、 Input 1では0.0115、 Input 2では0.0609でした。 この結果から、 低Counts領域の値の扱いに影響されるシグナルと、 前処理を変えても維持される比較的安定したシグナルを区別できました。
もしTPMだけをインプットしていた場合、 B-cellsを有意とする判断が 低Counts領域の不安定さに依存していないかを 検証できなかった可能性があります。
作図方法もAIと相談しながら修正する
最初に、各サンプルのxCell scoreを積み上げ棒グラフとして表示し、 この図をAIに提示して問題がないか尋ねました。 すると、AIは、xCell scoreは細胞割合ではなく、 複数の細胞種のスコアを合計しても 組織全体の構成比にはならないので、積み上げ棒グラフは不適切だと指摘しました。 AIは、代わりに 患者ごとにnon-lesional皮膚とlesional皮膚を隣接させた ヒートマップを提案しました。
さらに、AIは、細胞種ごとのZ-score化を提案しました。 しかし、Z-score化すると、大きなスコアを持つ細胞種と、 ほぼ0付近で小さく変動する細胞種の差が見えなくなります。 そこで、この提案は却下しました。
一方、表に残す主要な細胞種と、 それらを並べる順番についての提案は採用しました。 そのうえで、xCell scoreの実値を共通の色スケールで表示し、 極端に大きな値だけで色の濃さが決まらないよう、 全表示値の95th percentileをカラースケールの上限に設定しました。 できたのが下の図です。
GSE121212の各サンプルについて、 Input 2から算出した主要なxCell enrichment scoreを示した。 各患者のnon-lesional皮膚とlesional皮膚を隣接して表示している。 セル内の数値はxCell scoreの実値であり、 全表示値の95th percentileをカラースケールの上限に設定した。 ただし、xCell scoreは細胞割合ではなく、 異なる細胞種間の細胞数や存在割合を直接比較できる値ではない。
AIによって解析は速くなるが、学ぶべきことは判断へ移る
今回の解析では、deconvolutionツールの比較、 入力データの設計、Rスクリプトの生成と修正、 統計結果の確認、作図方法の改善にAIを活用しました。 これにより、従来では考えられないほど短時間で、 一連の解析を進めることができました。
しかし、AIの最初の提案をすべて採用したわけではありません。 TPMだけを使用することは避け、 低Counts領域の影響を確認するために二つの入力を作成し、 最終的な図ではZ-score化も採用しませんでした。
AIによってコーディングや実装が容易になるにつれて、 データ解析を学ぶ際の優先度も、 コードを書くことに大半の労力を費やすことから、 解析プロセス全体を理解し、 どこを精査すべきかを見極め、 AIの最初の提案を論理的かつ批評的に評価することへ 移していく必要があります。 そのためには、 インプットする値の性質、 統計手法の前提、 解析ツールの限界、 得られた出力値の意味など、 データ解析に関する幅広い知識が必要です。
AIの主な価値は、 正しい解析方法を自動的に決めることではなく、 選択肢を検討し、実際に解析を試し、 対話を重ねながら結果を改善できることにあります。
AIによってコードや結果を短時間で得られるようになるほど、 提案された解析が適切かどうかを判断する能力は、 これまで以上に重要になります。
