前回の記事 「RNA-Seq発現差解析でedgeR・DESeq2・t検定を比較した事例（2）｜biopsy由来・中規模データでの検証」では、edgeR、DESeq2、t検定で抽出される遺伝子の違いを比較しました。その結果、手法間の違いの多くは低Count領域に集中しており、十分なシグナルを持つ領域では、独立2群比較における3手法の結果は大きくは変わりませんでした。 特に3手法の共通部分では、明確な発現差が確認できました。 ただし、non-lesional群とlesional群を独立した2群として扱いました。
一方、同一患者から得られたnon-lesionalとlesionalのペア構造を考慮することで、個体差を差し引き、 non-lesionalからlesionalへの変化をより直接的に評価できると期待されます。 この考え方は、一般的な統計解析ではpaired t-testに対応します。 一方、edgeRやDESeq2では、patientをdesign formulaに含めることで、 患者ごとの差をブロック因子として扱います。 直感的には、patientをブロック因子として使えば、 個体差が差し引かれ、non-lesionalからlesionalへの変化に注目できるため、 検出力が上がるように思えます。 しかし、edgeRやDESeq2の統計モデルは、必ずしもそう単純ではありません。
それでは、ペア構造を考慮して計算されたDEGリストを比較してみましょう。
本記事で見る範囲：シグナル領域に限定した比較
前回の記事で示したように、RNA-Seqの発現差解析では、 通常のp値ベースの検定結果を解釈しやすいのは、主にシグナル領域の遺伝子です。 低Count領域では、測定値の不安定性やON/OFF型の発現パターンが混在し、 通常の検定結果をそのまま解釈しにくくなります。
そこで本記事では、non-lesional群とlesional群のそれぞれで、半数以上のサンプルにおいてGene Countsが20以上である遺伝子をQC1リストとしました。このQC1リストを、シグナル領域の遺伝子群として扱います。
ただし、edgeRとDESeq2の計算時点では、全遺伝子を対象にp値を計算しました。その後、得られたp値リストをSubio Platformに取り込み、 QC1リストと組み合わせることで、 シグナル領域に属する有意遺伝子だけを取り出しました。
この手順を用いた理由は、edgeRやDESeq2のp値計算前に遺伝子集合を絞ると、 解析対象となる遺伝子の母集団が変わり、 分散推定や多重検定の前提も変わるためです。 低Count遺伝子のフィルターは、実用的には有効な前処理ですが、 同時に、解析対象となる遺伝子集合を変える処理でもあります。
本記事の目的は、edgeR/DESeq2の標準的な最終DEGリストを作成することではなく、 同じGSE121212データに対して、各手法がどのような遺伝子を有意と判定するかを比較することです。 そのため、edgeR/DESeq2のp値計算は全遺伝子を対象に行い、 解釈段階でシグナル領域に限定しました。
また、本記事では手法間の判定傾向を比較するため、 各手法のp<0.05のリストを用います。 FDRは最終的なDEG候補を絞り込む際には重要ですが、 ここではまず、手法ごとの有意判定の違いを可視化して比較することを目的としています。
paired t-testは、edgeR/DESeq2の有意遺伝子をほぼ包含した
paired t-test、patientを含むedgeR、patientを含むDESeq2の結果を比較しました。 edgeRおよびDESeq2でp<0.05となった遺伝子は、 paired t-testでもほぼ有意となっていました。 一方で、paired t-testは、edgeRやDESeq2では有意にならない多数の遺伝子も抽出しました。
Fig1: patientを考慮したedgeR、DESeq2、paired t-testの比較。 QC1シグナル領域に限定すると、 edgeRおよびDESeq2で有意となった遺伝子は、 paired t-testでもほぼ有意となっていました。 一方、paired t-testは、より小さな変動幅の遺伝子まで広く有意と判定しました。
散布図では、paired t-testで有意となった遺伝子が、 対角線に近い領域にも広く分布していました。 これは、paired t-testが大きなfold changeだけでなく、 患者内で一貫した小さな変化にも敏感であることを示しています。
この結果を、paired t-testがedgeR/DESeq2より優れている、 あるいはedgeR/DESeq2が遺伝子を見逃している、 と単純に解釈するのは適切ではありません。 paired t-testとedgeR/DESeq2は、どちらもペア構造を考慮できますが、 評価しているデータ形式とノイズモデルが異なります。
paired designにすると、t-testとedgeR/DESeq2で逆方向の変化が起きた
上記の違いを理解するために、前回記事の2 groups比較と、今回のpaired designの結果を比較しました。
t-testでは、2 groups t-testからpaired t-testに変更することで、 有意遺伝子数が大きく増加しました。 これは、患者間の基礎発現量の違いを差し引くことで、患者内のnon-lesionalからlesionalへの変化方向が揃った遺伝子への感度が高くなった結果と考えられます。
一方で、edgeRおよびDESeq2では、 patientをdesign formulaに含めたpaired designにすると、 2 groups比較よりも有意遺伝子数が減少しました。
Fig2: 2 groups比較とpaired designの比較。 t-testではpaired t-testにすることで有意遺伝子数が大きく増加しました。 一方、edgeRおよびDESeq2では、patientをdesign formulaに含めることで、 2 groups比較よりも有意遺伝子数が減少しました。
この結果は、paired designが単純に検出力を上げる処理ではないことを示しています。 edgeRやDESeq2では、Gene Countsを負の二項分布モデルで扱い、 平均Count、分散、library size、遺伝子ごとのばらつき、 およびpatient項を含むGLMモデルの不確実性を考慮して検定します。 そのため、患者内で変化方向が揃っていても、 Countデータとして十分に安定した差と評価されない遺伝子は、 有意になりにくくなります。
paired t-test onlyの遺伝子群は、単なるノイズなのか？
paired t-testのみで有意となった遺伝子群は、 edgeRやDESeq2では有意と判定されていません。 そのため、最終的な強いDEG候補として扱うには注意が必要です。
しかし、この遺伝子群を完全に捨ててよいかというと、 そう単純ではありません。 Fig1で見たように、これらの遺伝子は変動幅こそ小さいものの、 QC1シグナル領域に属しており、 paired t-testでは患者内の一貫した変化として検出されています。
そこで、同じ患者のnon-lesionalに対するlesionalのlog2 ratioを計算し、 ヒートマップで確認しました。 各列はpatientを表し、 各セルはそのpatientにおけるlesional/non-lesionalのlog2 ratioを示します。
Fig3: 同じpatientのnon-lesionalに対するlesionalのlog2 ratioを用いたヒートマップ。色は、各patientにおけるlesional/non-lesionalのlog2 ratioを表します。intersectionの遺伝子群では、明瞭な上昇・低下パターンが見られました。 paired t-test onlyの遺伝子群では変動幅は小さいものの、 non-lesionalからlesionalへの変化方向が一定していることが確認できました。
edgeR、DESeq2、paired t-testのすべてで有意となったintersection遺伝子群では、 患者内の変化が大きく、 ヒートマップ上でも明瞭な上昇・低下パターンとして確認できました。
一方、paired t-test onlyの遺伝子群では、 変動幅は小さいものの、 多くの遺伝子でnon-lesionalからlesionalへの変化方向が一定していました。 このことから、paired t-test onlyの遺伝子群は、 単なるランダムノイズの集合ではなく、 小さいが一貫した発現変化を含む遺伝子群として考えることができます。
paired t-test only遺伝子群のサブパターン
paired t-test onlyの遺伝子群をさらに詳しく見ると、 上昇方向と低下方向の大きなパターンに加えて、 複数のサブクラスターに分かれて見えました。
Fig4: paired t-test only遺伝子群のクラスターパターン。 変動幅は小さいものの、上昇方向、低下方向、 および患者間で変化量のばらつきが見られるサブパターンに分かれていました。
Cluster 1とCluster 2は、AD lesionalで上昇する遺伝子群として見えました。 一方、Cluster 3、4、5は、AD lesionalで低下する遺伝子群として見えました。 特に低下方向の遺伝子群は、単一の均質なパターンではなく、 複数のサブクラスターに分かれて見えました。
この分類は、厳密な細胞種同定を目的としたものではありません。 ヒートマップ上のクラスタリング結果をもとに、 paired t-test only遺伝子群に含まれる代表的な発現パターンを探索的に整理したものです。
paired t-test only遺伝子群にも、機能的なまとまりが見える
下の表は、各クラスターの代表遺伝子名と、 ADのnon-lesionalとlesionalの比較において Cluster 1、2はlesionalで上昇し、 Cluster 3、4、5はlesionalで低下するという情報をChatGPTに入力し、 クラスターごとの特徴を整理させた結果を、そのまま転載したものです。
この表は、代表遺伝子名から機能的傾向を整理したものであり、 細胞種を直接同定したものではありません。 bulk RNA-Seqでは、細胞種組成、細胞状態、発現プログラムの変化が重なって観察されるため、 ここではクラスターの解釈候補として扱います。
| Cluster | AD lesionalでの変化 | 主な特徴遺伝子 | 示唆される機能的傾向 | 関連が示唆される細胞成分・状態 |
|---|---|---|---|---|
| 1 | ↑ | HNRNPC, HNRNPK, U2AF1, SRSF2, PSMA/PSMB群, NUP群, ATP5F1群, NDUFS群 | RNA処理、翻訳、核輸送、タンパク質分解、エネルギー代謝の亢進 | 活性化ケラチノサイト、増殖性表皮 |
| 2 | ↑ | Cluster 1と近い遺伝子群 | Cluster 1と同様のハウスキーピング・代謝活性化プログラム | 活性化ケラチノサイト、増殖性表皮 |
| 3 | ↓ | IL18, KIT, SOX10, EPAS1, EFNB2, FILIP1L, SASH1, PLXDC2 | 血管関連、低酸素応答、免疫調節、特殊間葉系プログラム | 血管周囲細胞、特殊ストローマ、メラノサイト関連成分の可能性 |
| 4 | ↓ | YY1, APC, FBXW7, AGO1, SMAD2, LATS1, ATM, SETD2, CREBBP, PCDH群 | 分化維持、転写制御、細胞接着、恒常性維持に関わるプログラム | 分化した真皮細胞、神経関連シグナル、恒常性維持プログラム |
| 5 | ↓ | EGFL7, PROX1, AQP1, ACKR3, PDGFRA, HGF, VCAN, SVEP1, ABI3BP, PTGIS | 血管・リンパ管機能、ECM維持、組織恒常性 | リンパ管内皮、血管内皮、線維芽細胞、血管周囲ストローマ |
上昇方向のCluster 1およびCluster 2では、 RNA処理、翻訳、タンパク質分解、エネルギー代謝に関わる遺伝子が多く含まれていました。 これは、AD lesionalにおける活性化ケラチノサイトや増殖性表皮に関連する変化として解釈できます。
一方、低下方向のCluster 3、4、5では、 血管、リンパ管、ECM、組織恒常性、分化維持に関わる遺伝子が目立ちました。 これらは、lesional部位での組織構造や細胞状態の変化を反映している可能性があります。
このように、paired t-test onlyの遺伝子群は、 fold changeは小さいものの、 完全なノイズとして捨ててしまうには惜しい構造を持っていました。 edgeRやDESeq2では有意と判定されなかったとしても、 シグナル領域内で患者内の一貫した変化を示し、 クラスター単位では機能的まとまりも見えるため、 解析結果の解釈では考慮に入れる価値があります。