ある遺伝子の発現量が、細胞の状態遷移の初期には低く、中盤で大きく上がり、終盤には再び低くなったとする。初期と終盤だけを取り出せば、両者の発現量はほとんど変わらない。二つの時点だけを比較する検定では、この遺伝子は「大きく変化していない」側へ分類される可能性がある。しかし細胞が初期状態から終盤の状態へ移る途中では、発現量が大きく上昇し、その後に元の水準まで低下している。始点と終点の差が小さいことと、途中で変化が起きなかったことは同じではない。
この違いは、細胞の分化や状態遷移を調べるときに実際の解析結果へ影響する。途中だけ発現する遺伝子を拾える検定と、始点と終点の差を強く評価する検定では、同じデータを入力しても上位に現れる遺伝子が変わり得るからである。しかも、途中の一過性変化が観測されても、その遺伝子が状態遷移を引き起こした原因だと直ちに言えるわけではない。まず必要なのは、状態遷移に沿ってどの遺伝子がどのような発現パターンを示したかを、途中の形を失わずに検出することである。
単一細胞 RNA-seq は、多数の細胞について遺伝子発現を細胞単位で測定する。一方で、一般的な単一細胞 RNA-seq から得られるのは、同じ一個の細胞を初期から終盤まで連続して追跡した時系列ではない。ある時点で回収された多数の細胞が、それぞれ異なる転写状態を持っている。そのため、細胞集団の中から状態の連続性を推定し、各細胞を発生、分化、応答などの進行に沿って並べる必要がある。
この推定された進行軸が擬似時間である。たとえば、発現状態が未分化な細胞に近い細胞を前方に、分化後の状態に近い細胞を後方に配置すれば、多数の細胞から状態遷移の順序を再構成できる。Trapnell らは、単一細胞を擬似時間上に並べることで、細胞運命決定に伴う遺伝子発現動態を再構成する枠組みを示した[1]。ここで得られる時間は実験装置で測った経過時間ではなく、細胞の発現状態から推定された位置である。1 時間、2 時間という物理時間ではなく、「この細胞はこの遷移のどのあたりにいると考えられるか」を表す軸として読む必要がある。
擬似時間が得られると、解析対象は単なる細胞間の発現量差から、擬似時間に沿った発現パターンへ広がる。ところが、遺伝子発現は右肩上がりや右肩下がりのような単純な形だけを取るとは限らない。途中だけ上がって元へ戻る、一度下がった後に上がる、複数の山を持つ、ある条件では早く上がり別の条件では遅れて上がる、といった形があり得る。擬似時間を推定できても、その上で何を「変化した遺伝子」と判定するかは別の統計問題として残る。
擬似時間上の変化をどのような量として表すかによって、強く検出される発現パターンは変わる。まず既存の検定がどの差を問うのかを確認し、そのうえで、発現曲線を別の表現へ変換する方法が一過性・非単調な変化をどこまで拾えるのかを、実データと統計的な限界の両方から見ていく。
1. 細胞の変化は、最初と最後だけではわからない
擬似時間に沿った遺伝子発現を考えるとき、最初に区別すべきなのは「状態の差」と「過程の差」である。状態の差は、たとえば初期細胞と終盤の細胞で発現量がどれだけ違うかを比較すれば測れる。過程の差を知るには、その二点の間で発現量がどのように動いたかまで見る必要がある。二つの遺伝子が同じ初期値から始まっても、途中の動きが異なれば、状態遷移との関わり方も同じとは限らない。
| 発現パターン | 初期 | 中期 | 後期 | 始点と終点だけを比較した場合 | 過程全体で観察できること |
|---|---|---|---|---|---|
| 遺伝子 A | 低い | 中程度 | 高い | 初期と後期の差が大きく、変化を検出しやすい。 | 状態遷移に沿って発現量が持続的に上昇している。 |
| 遺伝子 B | 低い | 高い | 低い | 初期と後期の差が小さく、変化が弱く見える。 | 中期に一過性の大きな発現ピークが存在する。 |
遺伝子 A は、初期から中期、後期へ進むにつれて発現量が増えている。始点と終点の二点だけでも低い状態から高い状態への移行が残るため、終点差を評価する方法でも変化を拾いやすい。遺伝子 B は事情が異なる。中期には遺伝子 A より大きな変化を示していても、後期には初期の水準へ戻っている。終盤の細胞だけを取り出せば、途中で起きた発現上昇に関する情報は失われる。
この違いは、単にグラフを何点表示するかという可視化上の話ではない。始点と終点の差を検定統計量として使えば、遺伝子 A のような変化には強く反応する一方、遺伝子 B の中期ピークは評価へ入りにくくなる。逆に、擬似時間全体の曲線形状を評価対象にすれば、終点では元へ戻った発現変化も検出対象に含められる。どの情報を統計量へ残すかが変わるため、最終的に「変化した遺伝子」として選ばれる集合も変わる。
この解析を実時間の時系列と同じものとして扱うことはできない。単一細胞 RNA-seq では、一つの細胞について時刻 0、時刻 1、時刻 2 の発現量を順番に測っているのではなく、それぞれ別の細胞から一度ずつ転写状態を取得する。その集団の中に初期に近い細胞、中間状態の細胞、終盤に近い細胞が含まれていると考え、発現状態の連続性から順序を推定することで、初めて一本の進行軸が得られる。
擬似時間は、この推定された順序へ連続的な位置を割り当てたものである。実時間の測定なら、10 分後と 20 分後の間隔が 10 分であることは観測条件から決まる。擬似時間では、隣り合う二つの細胞の距離も、どの細胞が先に位置するかも、観測された遺伝子発現から推定される。そのため、擬似時間上で見つかった発現変化には、元の RNA-seq の測定変動だけでなく、細胞をどの軌跡へ置いたかという推定も影響する。
軌跡推定の方法が一つに定まっていないのは、この推定対象が単純ではないためである。Saelens らは 45 種類の軌跡推定法を実データとシミュレーションデータで比較し、直線的な進行、分岐、複雑な軌跡など、対象とする構造によって各手法の適性が異なることを示した[2]。Street らの Slingshot は、その中でも複数の細胞系譜を推定し、それぞれに沿った擬似時間を与える方法である[3]。
たとえば、一つの前駆状態から二種類の細胞へ分化する場合、細胞を一本の直線上に並べるだけでは二つの進路を表現できない。まず共通する経路があり、途中で二本へ枝分かれするという軌跡を推定し、それぞれの枝に沿って擬似時間を与える必要がある。軌跡の形を誤れば、その後に並ぶ細胞の順序も変わる。順序が変われば、同じ遺伝子発現値でも「途中で上昇した」「後半で低下した」といった曲線の形が変わり、下流の発現差検定にも影響する。
細胞を擬似時間上へ配置できても、解析はそこで完了しない。次に決める必要があるのは、擬似時間上の何を差として検定するかである。初期と終盤の発現量を比べれば、持続的に上昇または低下する遺伝子を捉えやすい。曲線全体を推定すれば、中間のピークや谷を扱える。分岐ごとに曲線を分ければ、一方の細胞系譜だけで起きる変化を探せる。野生型とノックアウトの軌跡を比較すれば、ある一点の発現量だけでなく、発現が上がる時期や下がる時期まで条件差として扱える。
これらはすべて「擬似時間に沿った発現変化」を調べる解析だが、同じ問いではない。遺伝子 A と遺伝子 B の例で見たように、終点差を問う検定と途中の形を問う検定では、強く反応する発現パターンが違う。軌跡推定によって細胞の進行順序を作る工程と、その軸の上でどの発現変化を有意と判定するかを決める工程は分離して考える必要がある。以下で扱う方法が変更するのは後者であり、次に見るべきなのは、擬似時間上の発現曲線をどのような量へ変換して検定するかである。
2. 擬似時間を作った後にも、「どの変化を拾うか」という問題が残る
擬似時間が得られても、「どの遺伝子が状態遷移に沿って変化したか」は自動的には決まらない。下流解析では、同じ細胞、同じ発現量、同じ擬似時間を使っていても、何を差として検定するかによって結果が変わる。途中で一度だけ発現が上がる遺伝子と、終盤まで発現が上がり続ける遺伝子では、強く反応する検定が異なるからである。
tradeSeq は、この違いを具体的に示す既存法の一つである。各遺伝子のカウント値を負の二項分布で扱い、一般化加法モデルによって系譜ごとの発現量を擬似時間上の滑らかな曲線として推定する。そのうえで、擬似時間に沿った発現変化を見る associationTest、始点と終点を比較する startVsEndTest、複数の系譜で発現パターン全体を比較する patternTest、系譜の終盤を比較する diffEndTest など、異なる問いに対応する検定を用意している[4]。
| 検定で問う差 | 強く捉えやすい発現変化 | 弱く見える可能性がある発現変化 |
|---|---|---|
| 始点と終点の差 | 初期から終盤まで発現量の差が残る変化を捉えやすい。 | 中盤だけ上昇し、終盤には初期水準へ戻る変化は差が小さく見える。 |
| 擬似時間全体との関連 | 単調増減だけでなく、中間のピークや谷を含む変化も検出対象になる。 | 変化の振幅が小さい場合や、擬似時間との対応が不安定な場合には信号が弱くなる。 |
| 系譜間の曲線差 | 同じ遺伝子が二つの分岐で異なる時期や形で発現する場合を比較できる。 | 両系譜の曲線形状が近ければ、特定時点の小さな差だけでは強く検出されないことがある。 |
| 終盤の系譜差 | 異なる系譜が最終状態で異なる発現水準へ到達する場合を捉えやすい。 | 途中だけ系譜差が生じ、終盤で再び近づく場合には情報を取りこぼし得る。 |
この表で区別しているのは、同じ目的に対する四つの精度競争ではない。それぞれが別の帰無仮説、つまり「その種類の差はない」と置いた出発点に対して、別の発現差を問う検定である。たとえば発現量が「低い → 高い → 低い」と変化する遺伝子では、始点と終点の差は小さいままだが、擬似時間全体との関連は強くなり得る。逆に「低い → 中程度 → 高い」と持続的に上昇する遺伝子なら、途中の形を詳しく見なくても終点差が明瞭に残る。どちらの検定が正しいかではなく、何を検出したいかによって適切な問いが違う。
そのため、二つの方法から異なる遺伝子集合が得られても、それだけで一方を誤検出と判断することはできない。検定が違えば、同じ発現曲線のどの特徴を統計量へ残すかが変わる。始点と終点だけを使えば中間部分の情報は統計量から落ちる。曲線全体を使えば中間の形まで評価できる代わりに、曲線をどう推定するかという別の仮定が入る。P 値は、設定した帰無仮説の下で観測されたもの以上に極端な検定統計量が得られる確率として計算されるため、その大小だけを遺伝子固有の重要度のように読むことはできない。下流解析の結果を読むには、その P 値がどの帰無仮説に対して計算されたものなのかを確認する必要がある。
さらに、擬似時間解析では検定に入力する横軸そのものにも不確実性がある。実験で「培養開始から 6 時間後」と記録された時間なら、その値は観測条件として与えられる。擬似時間はそうではない。遺伝子発現データから細胞の軌跡と順序を推定し、その結果を時間軸として下流解析へ渡す。軌跡推定で細胞の順序が変われば、同じ遺伝子について描かれる発現曲線も変わり、その曲線に基づく検定結果まで変化し得る。
Song と Li が提案した PseudotimeDE は、この推定された時間軸を固定値として無条件に扱うことによって生じる統計的な問題を正面から扱っている。擬似時間を再推定する過程を含めて帰無分布を構成し、遺伝子発現と擬似時間の関連について P 値を適切に較正することを重視した方法である[5]。これは、擬似時間解析では「発現曲線をどの形で表現するか」と「その横軸をどこまで確かなものとして扱うか」が別々の論点になることを示している。
下流解析には、少なくとも二つの推定が連続して入る。まず多数の細胞から軌跡と擬似時間を推定し、その後で、その推定軸に沿って各遺伝子の変化を検定する。前段の順序推定がずれれば後段の発現パターンも変わり、後段で採用する検定が違えば同じ曲線から選ばれる遺伝子も変わる。擬似時間を一度作れば客観的な時系列が完成するわけではなく、軌跡推定と発現差検定の両方を含めて初めて「変化した遺伝子」という結果が成立する。
次章で扱う方法が変更するのは、このうち後段で発現パターンをどのように表現して検定するかという部分である。既存の回帰モデルへさらに別の曲線形状を追加するのではなく、擬似時間に沿った発現値を別の表現へ変換する。次に見るべきなのは、この表現変換によって、時間領域の曲線を直接比較する場合とは何が変わるのかである。
3. scLS は発現曲線を周波数領域へ移して調べる
Iuchi と Hamada が提案した scLS は、擬似時間に沿った遺伝子発現を Lomb–Scargle periodogram へ変換し、その結果を使って二種類の検定を行う[6]。第一の dynamic expression test は、一つの軌跡の中で、ある遺伝子の発現が擬似時間と無関係にばらついているのか、それとも擬似時間に沿った構造的な変化を持つのかを調べる。第二の shifted expression test は、野生型とノックアウトのような二条件について、それぞれの発現量が単純に高いか低いかではなく、擬似時間に沿った発現パターン全体が異なるかを調べる。前者が「この遺伝子は状態遷移の途中で変化するか」を問い、後者が「条件を変えると、その変わり方まで変化するか」を問う構成になっている。
この方法を理解するには、まず擬似時間上のデータが通常の時系列と違うことを押さえる必要がある。実験で 0 分、10 分、20 分、30 分と測定すれば、隣り合う観測点の時間間隔は一定である。一方、擬似時間は細胞集団から推定された進行軸なので、細胞が 0.0、0.1、0.2、0.3 のように等間隔で配置されるとは限らない。状態のよく観測された領域には多くの細胞が集まり、遷移が速い領域や取得しにくい状態では細胞が疎になることもある。その結果、各遺伝子について得られるのは、不等間隔な位置に発現量が並んだ系列になる。
Lomb–Scargle periodogram は、このような不等間隔の観測値から周波数成分を調べるために発展してきた。Scargle は 1982 年、不等間隔に観測されたデータについてスペクトル解析を行う際の統計的性質を整理した[7]。通常の離散フーリエ変換では、観測点が等間隔に並ぶことを前提に扱うのが自然だが、Lomb–Scargle 法は観測時点が不規則でも、その位置情報を保ったまま各周波数に対応する変動の強さを評価できる。擬似時間上の細胞配置が不等間隔であることは、この方法を単一細胞解析へ持ち込む技術的な理由の一つになる。
ただし、本稿でいう「周波数」を、遺伝子が 12 時間周期や 24 時間周期で発現しているかを測る量として読むと、scLS の狙いを取り違える。原論文では、Lomb–Scargle periodogram を生物学的な周期そのものの推定器としてではなく、擬似時間に沿った発現パターンを別の表現へ写すために使っている[6]。元の曲線では、横軸が擬似時間、縦軸が発現量であり、「どの位置でどれだけ発現したか」を直接読む。ピリオドグラムへ変換すると、その曲線に含まれるゆっくりした変化や、より細かな上下動が、それぞれ異なる周波数成分の強さとして表現される。
| 表現 | 主に保持する情報 | 読み取り方 |
|---|---|---|
| 擬似時間上の発現曲線 | どの擬似時間で発現量が高いか、低いかという位置と値の関係を直接表す。 | 中盤で上がる、後半で下がるといった変化を曲線の形として読む。 |
| Lomb–Scargle periodogram | その発現曲線に、どの程度の時間スケールの変動が含まれるかを周波数ごとの強さとして表す。 | 元の曲線形状を別の特徴空間へ変換し、変化の有無や条件間の違いを比較する。 |
たとえば、擬似時間の初期から終盤までゆっくり増え続ける遺伝子と、中盤だけ急に上がって再び下がる遺伝子では、元の時間領域でも曲線形状が異なる。scLS は、この違いを「直線」「一つ山の曲線」といった特定の関数形へ分類してから検定するのではなく、それぞれをピリオドグラムへ変換する。これにより、発現曲線をあらかじめ一種類の形へ固定しなくても、擬似時間に依存した構造的な変化を検出対象にできる。
dynamic expression test が利用するのは、この性質である。発現量が擬似時間と無関係に不規則に並んでいるなら、擬似時間上に一貫した変動構造は生じにくい。一方、初期から後期へ増える、中盤にピークを作る、いったん低下して再び上昇するといった発現パターンには、擬似時間上の位置と対応した構造がある。scLS はその構造をピリオドグラム上で評価するため、検出対象を「単調増加する遺伝子」や「始点と終点が違う遺伝子」に限定する必要がない[6]。
この違いは、第 1 章で示した遺伝子 A と遺伝子 B にそのまま対応する。遺伝子 A の「低い → 中程度 → 高い」という持続的な上昇も、遺伝子 B の「低い → 高い → 低い」という一過性のピークも、擬似時間に沿った構造的変化である。始点と終点の差だけなら両者への感度は大きく異なるが、scLS の dynamic expression test では、どちらも擬似時間に依存した発現パターンとして検定対象に入る。ここで得られる利点は、複雑な曲線をどれも同じものとして扱うことではなく、終点差という一つの特徴だけに情報を縮約しないことである。
shifted expression test では、同じ表現変換を二条件の比較に使う。たとえば野生型とノックアウトについて、ある遺伝子がどちらでも最終的には同じ発現量へ到達したとしても、野生型では早い段階にピークがあり、ノックアウトでは後半まで上昇しないという違いがあり得る。特定時点の発現量だけを比較すれば、選んだ時点によって差が大きく見えたり消えたりする。scLS は両条件の発現パターンをそれぞれパワースペクトルへ変換し、その違いを使って、擬似時間に沿った「変わり方そのもの」が条件間でずれているかを評価する[6]。
この設計は、分岐した軌跡を条件間で比較するときにも意味を持つ。複雑な木構造を持つ軌跡では、野生型の枝 A とノックアウトの枝 A’、枝 B と枝 B’ のように、対応する系譜を先に決めて個別比較する方法も考えられる。しかし、条件によって軌跡の形や枝数が変われば、一対一の対応そのものが曖昧になる。scLS の shifted expression test は、枝を逐一対応付けることを前提にせず、軌跡全体について遺伝子単位の差を先にスクリーニングできる[6]。多数の遺伝子から「条件によって発現動態が変わっていそうなもの」を絞り込む段階では、この性質が実用上の利点になる。
ただし、枝の対応付けを不要にした代わりに、局所的な位置情報まで同時に得られるわけではない。scLS がある遺伝子について条件差を検出しても、それだけでは「左側の枝でだけ変化した」「分岐直後から差が生じた」といった位置までは決まらない。軌跡全体をまとめた遺伝子単位の検定として感度を持たせることと、差が生じた枝や区間を特定することは別の解析課題だからである。後者には、系譜ごとの発現曲線を明示的に扱う解析を追加する必要がある。
scLS の設計を整理すると、変更点は「より複雑な回帰式を用意したこと」ではなく、検定へ渡す表現を変えたことにある。擬似時間上の発現量を直接比較する代わりに、不等間隔な位置情報を保ったまま Lomb–Scargle periodogram へ変換し、その周波数表現から動的発現と条件間差を調べる。この表現変換に実際の利点があるかは、数理的な柔軟性だけでは判断できない。次に確認すべきなのは、実際の単一細胞データで、既存の検定とは異なるどのような遺伝子発現パターンが検出されたかである。
4. 途中で上がって元に戻る遺伝子を、実データで違って拾った
scLS の特徴が最も具体的に見えるのが、ショウジョウバエ精巣の単一細胞 RNA-seq データを使った解析である。元データは Raz らが成体ショウジョウバエ精巣を単一細胞 RNA-seq と単一核 RNA-seq で解析したもので、生殖系列と体細胞系列について、発生段階に沿う細胞状態と分岐した系譜を含んでいる[8]。Iuchi と Hamada はこの公開データを用い、Monocle3 で軌跡と擬似時間を再推定したうえで、scLS の dynamic expression test と tradeSeq の startVsEndTest がどの遺伝子を強く検出するかを比較した[6]。
この比較では、単に二つの方法の P 値を並べて勝敗を決めているのではない。startVsEndTest は、擬似時間上の始点と終点で発現量が異なるかを問う検定である。たとえば発現量が「低い → 中程度 → 高い」と推移し、後半まで高い状態が残るなら、始点と終点の差は大きい。その差は検定が直接評価する対象なので、強い信号として現れやすい。一方、「低い → 高い → 低い」と途中だけ上昇する場合、変化の振幅が大きくても、終点が初期値へ近づけば始点と終点の差は小さくなる。発現曲線のどの部分を統計量へ残すかによって、同じ大きさの変化でも検出のされ方が変わる。
実データでも、この違いが遺伝子ごとの発現曲線として現れた。FBgn0040827、FBgn0031277、FBgn0030158、FBgn0262009 は、擬似時間の後半に発現量が大きく上昇し、その高い状態が終点まで残るパターンとして示された[6]。この形では、途中の詳細を使わなくても始点と終点の差が明瞭になる。startVsEndTest が強く反応するのは、検定の目的と実際の発現形状が一致しているためである。
これに対して、FBgn0034659、FBgn0035569、FBgn0051406、FBgn0034739、FBgn0036488、FBgn0264344、FBgn0038978、FBgn0035585 は、擬似時間の中盤付近で発現量が最大となり、その後に低下するピーク型のパターンを示した[6]。これらの遺伝子では、途中で明瞭な発現変化が存在していても、後半に発現量が低下するため始点と終点の差だけでは変化の全体を表せない。原論文では、この一群が startVsEndTest では相対的に弱い一方、scLS の dynamic expression test では強く検出される例として示されている[6]。
| 比較対象 | 検定が直接評価するもの | 強く拾いやすい発現パターン | 実データで示された例 |
|---|---|---|---|
| tradeSeq startVsEndTest | 擬似時間の始点と終点における発現量の差を評価する。 | 後半まで発現量の変化が残り、始点と終点の差が大きくなるパターンを強く捉えやすい。 | FBgn0040827、FBgn0031277、FBgn0030158、FBgn0262009 が示されている。 |
| scLS dynamic expression test | 擬似時間に沿った発現変化を Lomb–Scargle periodogram へ変換し、その構造的な変動を評価する。 | 単調増減に加えて、中盤でピークを形成した後に低下する一過性・非単調なパターンも検出対象に入る。 | FBgn0034659、FBgn0035569、FBgn0051406、FBgn0034739、FBgn0036488、FBgn0264344、FBgn0038978、FBgn0035585 が示されている。 |
この結果を「scLS は tradeSeq より優れている」と読むのは適切ではない。比較対象になっている startVsEndTest は、tradeSeq が備える複数の検定の一つにすぎない。tradeSeq には、擬似時間に沿った発現変化全体を調べる associationTest、複数の系譜で発現曲線全体を比較する patternTest、終盤の系譜間差を調べる diffEndTest なども用意されている[4]。startVsEndTest と scLS の違いは、同じ問いに対する単純な精度差ではなく、発現曲線のどの特徴を検定対象へ残すかの違いとして読む必要がある。
たとえば、FBgn0034659 のような中盤ピーク型の遺伝子について、startVsEndTest の P 値が相対的に大きかったとしても、それだけで「変化がない」と結論することはできない。その検定が答えているのは、始点と終点の発現量に十分な差があるかという問いだからである。scLS が強く検出した場合も、「この遺伝子が細胞状態遷移の原因である」と分かったわけではない。確認できるのは、擬似時間上に無秩序なばらつきとして説明しにくい発現変化があり、その変化を scLS の統計量が強く捉えたというところまでである。
この区別は、生物学的な候補遺伝子を絞り込む場面で効いてくる。終点で高発現になった遺伝子だけを候補にするなら、途中で一時的に働いて元へ戻る遺伝子は検索範囲から落ちる。一方、状態遷移の特定段階で短時間だけ誘導される遺伝子も候補に含めたいなら、中間部分を保持する検定が必要になる。解析法の選択によって、後続の実験で何を検証対象にするかまで変わり得るため、検定の違いは単なる統計処理上の選択では終わらない。
第 1 章で置いた「低い → 高い → 低い」という単純な例は、この Drosophila データで具体的な形を持つ。始点と終点が近い値なら、その二点からは途中のピークを再構成できない。scLS が実データで拾った中盤ピーク型の遺伝子群は、状態遷移を二点の差へ縮約すると失われる情報が実際に存在することを示している[6]。一過性の発現を解析対象に含めるなら、最終状態へ何が残ったかだけでなく、そこへ至る途中でどの発現パターンを通過したかを見る必要がある。
この実データ解析によって、scLS の周波数領域への変換が単なる数理上の言い換えではなく、実際に強く検出される遺伝子の違いへつながることが確認できる。ただし、ここまで扱ったのは一つの軌跡の中で発現が変化する dynamic expression である。scLS にはもう一つ、野生型とノックアウトのような二条件について、発現量の高低ではなく擬似時間に沿った変化の形そのものを比較する shifted expression test がある。次に見るのは、条件が変わったときに「どれだけ発現したか」ではなく「どう変わったか」を比較する場合である。
5. 野生型とノックアウトでは、「発現量の差」だけでなく「変わり方の差」も比較できる
一つの軌跡の中で発現が変化する遺伝子を探す dynamic expression test に対して、scLS の shifted expression test は二つの条件で発現パターンそのものが変わったかを調べる[6]。たとえば野生型とノックアウトを比較するとき、「ノックアウトでは発現量が低い」という差だけが条件差ではない。野生型では早い段階に発現が上がるのにノックアウトでは遅れる、野生型だけ中盤にピークを作る、両者とも最終的には同じ発現量へ戻る、といった違いもあり得る。ある時点の値だけを比較すれば、どの時点を選んだかによって差が現れたり消えたりする。shifted expression test が対象にするのは、この擬似時間に沿った「変わり方」の違いである。
| 条件差の形 | 一点の発現量比較で見えること | 軌跡全体を比較すると見えること |
|---|---|---|
| 常に一方が高い | 比較する時点を大きく選ばなくても差を捉えやすい。 | 軌跡全体を通じて発現水準がずれていることを確認できる。 |
| ピーク時期がずれる | ピークから外れた時点では差が小さく見えることがある。 | 発現上昇と低下の時期を含むパターンの違いとして捉えられる。 |
| 一方だけ中盤で上がる | 初期や終盤だけを比較すると差が消えることがある。 | 状態遷移の途中だけに存在する条件差を検出対象に含められる。 |
| 終盤だけ差が生じる | 終点付近を比較する検定と相性がよい。 | 軌跡全体の中で終盤に差が集中しているパターンとして扱える。 |
原論文は、この shifted expression test をアフリカトリパノソーマ Trypanosoma brucei の単一細胞 RNA-seq データに適用している。T. brucei は哺乳類宿主内で増殖する slender form から、ツェツェバエへの伝播に適した stumpy form へ状態を移す。Briggs らは 8,599 個の寄生虫について単一細胞の転写状態を測定し、細胞周期と slender form から stumpy form への分化に沿った状態遷移を再構成した[9]。このデータには、状態遷移を制御する RNA 結合タンパク質 ZC3H20 をノックアウトした条件も含まれている。
ZC3H20 はこの比較を考えるうえで都合のよい対象である。ZC3H20 はノックアウト対象そのものであるため、野生型との条件差が検出されることを期待できる。未知の候補遺伝子について「本当に差があるか」を推測する場合とは異なり、少なくとも ZC3H20 自体については、既知の実験操作に対して各検定がどのように反応するかを確認できる。
Iuchi と Hamada は、この ZC3H20 について scLS の shifted expression test と、tradeSeq の patternTest、earlyDETest、diffEndTest を比較した[6]。patternTest は軌跡に沿った発現曲線全体の違いを評価し、earlyDETest は分岐付近の早い段階に現れる系譜間差を調べ、diffEndTest は軌跡の終端付近における発現差を評価する。それぞれ同じ ZC3H20 を入力していても、検定が見ている軌跡上の範囲と曲線の特徴は異なる。
結果も一様ではなかった。ZC3H20 は scLS、patternTest、earlyDETest では非常に小さい P 値を示した一方、終点付近の差に重点を置く diffEndTest では有意性が相対的に弱かった[6]。ノックアウトされた遺伝子なのだから条件間の違いそのものが存在しないわけではない。それでも検定結果に差が生じるのは、ZC3H20 の条件差が軌跡のどこに、どのような形で現れるかを各検定が異なる方法で要約しているからである。
この結果は、P 値を「遺伝子に付随した固定的な重要度」のように読むと誤解しやすいことも示している。同じデータに同じ遺伝子が含まれていても、patternTest が問う「曲線全体は違うか」と、diffEndTest が問う「終盤の発現量は違うか」は別の仮説である。前者で極めて小さい P 値が得られ、後者でそれほど小さくならないことに矛盾はない。P 値は遺伝子そのものの属性ではなく、定義された帰無仮説と観測データの組み合わせから得られる量だからである。
shifted expression test の特徴は、この条件差を時間領域の対応点ごとの差として直接計算するのではなく、各条件の発現パターンを Lomb–Scargle periodogram へ変換して比較するところにある[6]。野生型とノックアウトで発現曲線の形が変われば、それぞれのパワースペクトルにも違いが生じる。その差を用いることで、「ある擬似時間で何カウント違ったか」ではなく、状態遷移全体を通じて発現動態が変化した遺伝子を候補として抽出できる。
この設計は、二条件の軌跡が単純な一本線ではない場合にも意味を持つ。野生型とノックアウトがともに分岐構造を持つとき、系譜ごとの比較では「野生型の枝 A はノックアウトのどの枝に対応するのか」を先に決める必要がある。条件によって分岐の形、枝の長さ、細胞密度が変われば、この対応は必ずしも自明ではない。対応付けを誤れば、本来異なる状態を同じ系譜として比較することになり、その後の遺伝子差も対応関係の仮定に引きずられる。
scLS はこの枝同士の一対一対応を前提とせず、各条件について得られた軌跡全体から遺伝子単位の shifted expression を検定できる[6]。数千から数万の遺伝子の中から、条件を変えたことで発現動態全体が変化した候補を先に絞り込みたい場合には、対応する枝を遺伝子ごとに調べる前段階として使える。分岐構造の違いを吸収したうえで、まず候補遺伝子の集合を得るという位置づけである。
その利便性には、別の解析課題が残る。scLS がある遺伝子について有意な shifted expression を返しても、「どの枝のどの区間で違いが生じたか」はその統計量だけからは分からない。野生型の左枝だけで発現が上がった場合も、複数の枝で小さな変化が積み重なった場合も、軌跡全体として差があれば候補になり得る。軌跡全体の差を枝対応なしで検出することと、差の発生位置を局所化することは別の解析課題であり、後者を知るには系譜ごとの発現曲線や分岐位置を明示的に調べる追加解析が必要になる[6]。
ZC3H20 の例から得られるのは、「scLS ならノックアウトを正しく見つけられる」というだけの確認ではない。同じ既知の条件差に対しても、軌跡全体を見るか、分岐初期を見るか、終点を見るかによって統計的な強さが変わるという具体例になっている。条件差を一つの時点に固定すれば、その時点に残らなかった変化は弱く見える。過程全体を比較すれば、その途中に現れる差も含められる一方、差が生じた位置については別の解析が必要になる。どちらの情報を必要とするかによって、選ぶ検定も変わる。
ここまでの Drosophila と T. brucei の解析を並べると、scLS が狙う範囲がはっきりする。dynamic expression test は「一つの状態遷移の途中で発現した遺伝子」を探し、shifted expression test は「条件を変えると、その状態遷移に沿った発現の仕方まで変わる遺伝子」を探す。いずれも、始点と終点の発現量だけでは消えてしまう情報を過程から拾うための検定である。ただし、拾える形が増えたことは、既存法より常に高い精度で生物学的に正しい遺伝子を特定できることを意味しない。次に確認する必要があるのは、シミュレーション上での性能比較と、その評価をどこまで実データへ一般化できるかである。
6. scLS は既存法を置き換える万能法ではない
新しい解析法を評価するとき、「既存法より高い AUROC を出したか」という比較は分かりやすい。しかし、scLS の位置付けをその一軸だけで判断すると、この手法が何を改善し、何を改善していないのかを取り違える。原論文は scLS を既存法の全面的な置き換えとは位置付けておらず、とくに一過性・非単調な発現パターンを含む候補遺伝子を高速に絞り込み、その後に枝ごとの詳細解析や別の統計手法へ渡す一次スクリーニングとしての利用を想定している[6]。前章までで見たように、scLS の強みは、すべての発現差を一つの基準でより正確に測ることではなく、始点と終点の差だけでは弱く見える変化を別の表現から拾えることにある。
原論文のシミュレーションでは、既知の発現パターンを持つ人工データを作り、scLS、tradeSeq、Lamian、PseudotimeDE などが真の発現差をどの程度識別できるかを AUROC で比較している[6]。ここで AUROC は、真に変化する遺伝子と変化しない遺伝子を、検定統計量や P 値によってどの程度正しく順位付けできるかを表す。1 に近いほど識別性能が高く、0.5 なら無作為な順位付けに近い。
比較結果は、scLS がすべての条件で最高だったという形にはなっていない。比較的単純なシミュレーション条件では複数の手法が高い AUROC を示し、条件によっては scLS も高性能だった。一方、PseudotimeDE の評価用に設計された、より難しいベンチマークでは PseudotimeDE が全条件で最も高い性能を示した[6]。この結果から言えるのは、scLS が既存法を一律に上回ったことではなく、複雑な発現パターンを拾う能力を持ちながら、少なくとも複数のシミュレーション条件で既存法と競争可能な識別性能を維持したという範囲までである。
性能が条件によって変わる理由の一つは、下流検定の前に使う擬似時間そのものが推定結果だからである。原論文では、真の擬似時間を与えた場合だけでなく、Monocle3 や Slingshot で推定した擬似時間を入力した場合も比較している[6]。真の順序を知っているシミュレーションでは発現曲線を正しい軸上で評価できるが、推定擬似時間では細胞の順序や位置に誤差が入る。前段で細胞の並びがずれれば、中盤のピークが広がったり、上昇と低下の位置が移動したりし、その後に計算するピリオドグラムや回帰曲線も変わる。下流検定の性能は、検定法だけではなく、軌跡推定との組み合わせによって決まる。
| 評価軸 | scLS で確認されたこと | 評価時に残る条件 |
|---|---|---|
| 遺伝子識別性能 | 複数のシミュレーション条件で既存法と競争可能な AUROC を示し、一過性・非単調なパターンも検出対象にできた。 | PseudotimeDE の難しいベンチマークでは PseudotimeDE が全条件で最高となり、scLS が常に最良だったわけではない。 |
| 擬似時間への依存 | Monocle3 や Slingshot から得た擬似時間でも解析できた。 | 軌跡推定法が変わると細胞の並びも変わるため、同じ発現データでも下流性能が変化する。 |
| 発現パターンの柔軟性 | 単調増減だけでなく、中盤ピークなどの一過性・非単調な変化を強く検出する例が示された。 | 柔軟に拾えることは、検出した遺伝子が生物学的な原因遺伝子であることを保証しない。 |
| 分岐軌跡への適用 | 枝を条件間で一対一に対応付けず、軌跡全体について遺伝子単位のスクリーニングを行える。 | 有意差がどの枝のどの区間から生じたかは、scLS の統計量だけでは局所化できない。 |
もう一つ具体的な利点として、計算時間がある。原論文では 1000 細胞、2000 遺伝子からなる線形軌跡のシミュレーションを用い、1 コア環境で各手法の実行時間を比較している[6]。dynamic expression test では scLS が 1 秒未満で完了し、比較対象となった tradeSeq の検定や Lamian より短時間だった。shifted expression test でも、tradeSeq の各検定と同程度か、それより短い実行時間を示している。
この速度差は、一次スクリーニングという用途と結び付けて評価する必要がある。単一細胞 RNA-seq では、一つの遺伝子だけを調べるのではなく、数千から数万の遺伝子を同じ軌跡に沿って検定する。さらに条件比較や複数の軌跡推定法を試せば、同じデータに対する検定回数も増える。1 遺伝子ごとの計算が複雑でなくても、遺伝子数と解析条件が積み上がれば総計算量は大きくなる。その前段階で候補遺伝子を短時間に絞れるなら、その後に計算負荷の高い枝別解析や可視化を適用する対象を減らせる。
ただし、1000 細胞・2000 遺伝子・1 コアという論文中のベンチマーク条件から、任意の実データ規模で同じ速度差が維持されるとは言えない。細胞数が増えればピリオドグラムを計算する入力点も増え、遺伝子数が増えれば検定回数も増える。分岐構造の複雑さ、メモリ使用量、並列化の有無、前処理の実装も実運用時間に影響する。論文で確認できるのは、指定されたシミュレーション条件では scLS の計算時間が短かったという事実であり、それ以上の規模での速度は別途測定する必要がある。
精度比較にシミュレーションが使われるのには理由がある。Drosophila や T. brucei の実データでは、各遺伝子について「これは真に擬似時間依存である」「これは真に変化しない」という完全な正解表を持てない。発現曲線を見れば候補は分かっても、全遺伝子について真偽を確定した正解ラベルは存在しない。その状態では、ある手法が拾った遺伝子の数を数えても、それが高感度なのか誤検出が多いのかを厳密には区別できない。
シミュレーションでは、この制約を逆向きに処理する。最初に「どの遺伝子を変化させるか」「軌跡をどの形にするか」「細胞をどこへ配置するか」を決め、その条件から人工的なカウントデータを生成する。生成側が正解を知っているため、解析法がどの程度その正解を回収できたかを AUROC などで定量比較できる。
Splatter は、実際の単一細胞 RNA-seq データから平均発現量、分散、ゼロ値の頻度などの特徴を推定し、それらを反映した人工データを生成するシミュレーション基盤である[10]。PROSSTT はさらに、既知の分岐構造と擬似時間を持つ複雑な分化過程を設定し、その軌跡に沿って遺伝子発現が変わる単一細胞データを生成できる[11]。この種のシミュレーションを使えば、「本当の軌跡」「本当に変化する遺伝子」「真の擬似時間」を知った状態で、軌跡推定法と下流検定の両方を評価できる。
| シミュレーションで固定できるもの | そこから評価できること | 実データへ移したときに残る不確実性 |
|---|---|---|
| 真に変化する遺伝子 | 検出順位から AUROC を計算し、手法間の識別性能を比較できる。 | 実生物では全遺伝子の真の発現差を完全には確定できない。 |
| 真の軌跡構造 | 推定された分岐や細胞順序が正解にどれだけ近いかを測定できる。 | 実データでは推定軌跡そのものの正解を直接観測できない場合が多い。 |
| 真の擬似時間 | 真の時間軸を与えた場合と、推定擬似時間を使った場合の性能差を分離できる。 | 実データの擬似時間は潜在的な進行軸であり、物理時間との一致は保証されない。 |
| ノイズ生成モデル | ノイズ量や発現パターンを制御し、条件ごとの性能低下を再現できる。 | 実細胞にはシミュレーションで明示していない技術変動、生物学的異質性、標本差が含まれ得る。 |
この評価方法には明確な限界もある。シミュレーションで AUROC が高いということは、そのデータ生成モデルが作った発現パターンをうまく回収できたことを意味する。実際の細胞では、転写バースト、細胞周期、個体差、バッチ効果、測定深度、未知の分岐、複数の生物学的過程の重なりなどが同時に存在する。シミュレーションに含まれていない変動は、どれだけ AUROC を精密に比較しても評価対象にはならない。原論文自身も、人工データが実データより単純でノイズが少ない可能性を限界として挙げている[6]。
実データとシミュレーションは、ここで別の役割を持つ。Drosophila の解析は、中盤でピークを作る具体的な遺伝子群を scLS が強く拾うという、方法の挙動を現実のデータ上で示す。T. brucei の ZC3H20 ノックアウトは、既知の条件差に対して複数の検定が異なる反応を示す例になる。一方、AUROC の比較には真の正解が必要なので、定量的な識別性能の評価はシミュレーションが担う。実データだけでも、人工データだけでも、scLS の性質を十分には評価できない。
この区別を置くと、scLS の評価軸は一つではなくなる。第一に、一過性・非単調な発現パターンを検出対象へ含められる。第二に、条件間で枝を明示的に対応付けずに全体スクリーニングを行える。第三に、論文中の条件では比較的短時間で多数の遺伝子を処理できる。第四に、識別性能は複数の既存法と競争可能だが、ベンチマークによって最良法は変わる。これらを合わせると、scLS は「最も高精度な万能法」というより、複雑な発現動態を高速に候補化するための補完的な選択肢として位置付ける方が原論文の結果に合う。
残る課題は、シミュレーションでは真の値として与えられる擬似時間や生物学的反復が、実データでは推定や標本構造の一部になることである。scLS が複雑な発現曲線を拾えることと、その P 値があらゆる実験設計で十分に較正されることは別の問題になる。次章では、擬似時間を同じデータから推定して再び検定へ使うこと、複数の生物学的反復をどう扱うか、分岐上の差をどこまで局所化できるかという、原論文自身が残した統計的な制約を確認する。
7. 擬似時間そのものが推定値である以上、検定にも限界が残る
scLS の限界は、原論文の考察で比較的明確に整理されている。前章までで見た利点は、擬似時間に沿った複雑な発現パターンを拾いやすいこと、条件間で枝を一対一に対応付けず全体差を検索できること、論文中の条件では計算時間が短いことである。一方、これらの利点は、擬似時間そのものが推定値であること、単一細胞データが複数の個体や標本から構成されること、同じ発現行列を軌跡推定と仮説検定の両方に使うことによって生じる統計的な不確実性を消すものではない。scLS がどの形の発現変化を検出できるかと、その P 値をどこまで確かな推論として読めるかは分けて考える必要がある。
7.1 細胞が多くても、生物学的反復が多いとは限らない
第一の制約は生物学的反復（biological replicate）、すなわち異なる個体、培養、処置単位などから独立に取得した生物学的反復を、現在の scLS が明示的にはモデル化しないことである[6]。単一細胞 RNA-seq では、一つのサンプルから数千、数万の細胞が得られるため、観測数だけを見ると非常に大きなデータに見える。しかし、同じ個体から得た 5000 細胞は、5000 個の独立した個体を測ったことにはならない。細胞は同じ個体、同じ処置、同じ実験環境を共有しており、サンプル単位の変動を持つからである。
たとえば野生型 1 個体から 5000 細胞、ノックアウト 1 個体から 5000 細胞を取得した場合、細胞数だけなら各条件に十分な観測があるように見える。しかし両条件の差には、遺伝子ノックアウトの効果だけでなく、個体差、採取時の状態、ライブラリ作製、バッチなど、その 2 個体に固有の差も混ざり得る。細胞をすべて独立な反復として扱えば、同じサンプルに属する細胞の類似性を無視し、実際より情報量が多いように評価する危険がある。
Hou らは、複数の単一細胞 RNA-seq サンプルを対象とする擬似時間差解析の統計枠組みを提案し、Lamian によってサンプル間変動を含めた擬似時間依存の差を検定している[12]。ここで扱われているのは、「細胞数を増やせば解決する」問題ではない。細胞はサンプルの内側に入れ子になっているため、細胞レベルの変動とサンプルレベルの変動を分離してモデル化する必要がある。scLS 原論文も生物学的反復を考慮した推論を今後の課題として挙げている[6]。
この制約は、scLS の使用対象を考えるときに実務的な意味を持つ。探索的に一つのデータセットから一過性発現の候補を順位付けする用途と、複数個体に共通する条件差として統計的に一般化する用途では、要求される推論の強さが違う。前者では遺伝子候補の一次スクリーニングとして利用できても、後者では生物学的反復を考慮した解析や独立データでの再現確認を追加する必要がある。
7.2 軌跡全体の差を拾えても、差が生じた枝までは分からない
第二の制約は、分岐した軌跡に対する軌跡全体を対象とする検定と枝別の解釈の違いである。scLS は、野生型とノックアウトで枝を一対一に対応付けなくても、軌跡全体について遺伝子単位の dynamic expression や shifted expression を計算できる[6]。条件によって枝数や分岐形状が変わる場合に、対応関係を事前に固定しなくてよいことはスクリーニング上の利点になる。
ただし、軌跡全体を一つの統計量へまとめると、差が生じた場所の情報はその統計量だけから復元できない。たとえば、ある遺伝子が左側の系譜だけで大きく上昇した場合と、複数の系譜で小さな変化を示した場合のどちらでも、軌跡全体として十分な変化があれば scLS は候補として検出し得る。検定結果が有意であることから、「第 2 分岐の右側でのみ変化した」「分岐直後に発現がずれた」といった局所的な結論までは導けない。
軌跡全体を対象とする検定と、差が生じた枝や位置を特定する解析は役割が異なる。最初に数千遺伝子から条件依存の発現動態を持つ候補を絞り、その後で候補遺伝子について系譜ごとの発現曲線を描く、分岐ごとの検定を行う、細胞型マーカーと照合するといった二段階の解析が必要になる。scLS の出力は「どの遺伝子を詳しく見るか」を決める材料にはなるが、「軌跡上のどこで何が起きたか」を単独で説明する結果ではない。
7.3 擬似時間を作ったデータで、再び擬似時間との関連を検定する
第三の制約は、擬似時間推定と遺伝子検定に同じ発現データを使うことによるデータの二重利用（double dipping）である。単一細胞 RNA-seq では、遺伝子発現行列から細胞間の類似性を計算し、その構造を使って軌跡や擬似時間を推定する。その後、同じ発現行列の各遺伝子について「推定された擬似時間と関連するか」を検定する。この二段階では、検定の横軸となる擬似時間が、検定対象である遺伝子発現と独立に与えられたものではない。
単純化すると、まずデータを見て「この順序がもっともらしい」と時間軸を作り、その同じデータを使って「この時間軸に沿って変化しているか」を調べている。軌跡推定の段階ですでに発現構造を利用しているため、その構造を固定された外部情報として後段の検定へ渡すと、通常の検定が想定する独立性が崩れる場合がある。その結果、帰無仮説が正しい遺伝子についても P 値が必要以上に小さくなり、第 1 種過誤、すなわち本当は関連がない遺伝子を有意と判定する確率を十分に制御できない可能性が生じる。
Neufeld らは、この問題を単一細胞 RNA-seq における潜在変数推定後の推論として定式化している[13]。細胞の潜在状態、低次元表現、クラスター、擬似時間のような変数を発現データから推定した後、その推定結果と同じ発現データの関連を検定すると、推定工程を無視した通常の P 値では誤差制御が崩れ得る。Neufeld らは Poisson モデルの下でカウントを分割し、一方を潜在変数の推定、もう一方を推論に使うカウント分割（count splitting）を提案している[13]。
scLS 原論文も、この論点を未解決の統計的課題として認めている。scLS は軌跡推定の方法そのものではなく、すでに得られた擬似時間に沿って発現変化を検定する下流手法である。しかし実際の利用では、同じ単一細胞データから Monocle3 や Slingshot などで軌跡を推定し、その結果を scLS へ入力することが多い。そのため、scLS の検定式だけを見て P 値の妥当性を判断するのでは不十分であり、擬似時間がどのデータから、どの手順で推定されたかまで含めて推論系全体を見る必要がある[6]。
| 解析上の制約 | 直接起きる問題 | 結果を読むときの帰結 |
|---|---|---|
| 生物学的反復を明示的に扱わない | 同じサンプルに属する多数の細胞を、独立した生物学的反復と同じようには扱えない。 | 個体や独立サンプルに一般化した条件差を主張するには、反復構造を扱う解析や再現実験が必要になる。 |
| 軌跡全体の統計量を返す | 有意差がどの系譜のどの区間から生じたかを直接局所化できない。 | 候補抽出後に系譜別の解析や発現曲線の確認を追加する必要がある。 |
| 擬似時間が同じ発現データから推定される | 軌跡推定と発現差検定でデータを二重利用し、通常の P 値の誤差制御が崩れる可能性がある。 | P 値を推定工程から独立した確定的な証拠として読むことはできず、推定不確実性を扱う方法が必要になる。 |
| 発現変化を統計的に検出する | 検出された遺伝子が状態遷移の原因なのか、結果なのかを区別できない。 | 因果機構を主張するには遺伝子操作、機能実験、独立データなど別種の証拠が必要になる。 |
7.4 発現が変わったことと、状態遷移を起こしたことは別である
第四の制約は統計手法固有というより、scLS の結果を生物学的に解釈するときの射程に関わる。dynamic expression test が有意なら、その遺伝子の発現には擬似時間に沿った構造的な変化がある。shifted expression test が有意なら、二条件の間で擬似時間に沿った発現動態が異なる。どちらも「発現パターンに統計的な構造がある」ことを検出しているのであって、その遺伝子が状態遷移を開始させた原因であることを直接検定しているわけではない。
たとえば、ある転写因子 X が細胞分化を開始し、その下流で遺伝子 Y の発現が中盤だけ上昇したとする。scLS は Y の「低い → 高い → 低い」という一過性パターンを強く検出できる可能性がある。しかし、この結果だけでは Y が分化を駆動したのか、X によって引き起こされた応答なのかを区別できない。さらに、細胞状態そのものが変化した結果として Y の発現が上がっただけという可能性も残る。擬似時間との関連は、因果の向きを決めない。
ZC3H20 の例で条件差を比較できたのは、ZC3H20 を実際にノックアウトするという介入が元データに含まれていたからである。その介入から ZC3H20 条件と野生型の違いを調べることはできるが、scLS が同じデータから新たに拾ったすべての候補遺伝子について、同じ因果関係が証明されたことにはならない。候補遺伝子が状態遷移を駆動するかを確かめるには、ノックアウトやノックダウン、過剰発現、時系列介入、機能評価など、発現パターンの観察とは別の証拠が必要になる。
ここまでの制約を合わせると、scLS の役割を「状態遷移の鍵となる遺伝子を自動的に確定する方法」と表現するのは強すぎる。より正確なのは、擬似時間に沿った一過性・非単調な発現変化や、条件によって変化の仕方が異なる遺伝子を、軌跡全体から効率よく候補として順位付けする方法と位置付けることである。候補抽出の後には、どの系譜で差が生じたか、複数の生物学的反復でも再現するか、推定擬似時間の不確実性を考慮しても有意か、遺伝子操作によって状態遷移が変わるかという別の検証が続く。
この位置付けにすると、scLS の強みと限界は矛盾しない。周波数領域へ表現を変えることで、終点差だけでは弱く見える一過性発現を候補として浮かび上がらせられる。一方、その候補がどの枝で働くか、生物学的反復を越えて再現するか、状態遷移の原因かどうかは別の問いとして残る。scLS が広げたのは「何を候補として見つけられるか」の範囲であり、「見つかった候補について何まで確定できるか」の範囲まで同時に広げたわけではない。
8. 状態だけを見るか、過程全体を見るかで、見える遺伝子は変わる
ここまでの Drosophila と Trypanosoma brucei の例をつなぐと、scLS が変えたのは細胞の測定方法ではない。入力されるのは、すでに単一細胞 RNA-seq で得られた遺伝子発現データと、そこから推定された擬似時間である。変更されたのは、その擬似時間上の発現変化を何として比較するかである。発現曲線を時間領域のまま始点と終点で比較する方法もあれば、曲線全体を回帰モデルで表す方法もある。scLS はそこへ、Lomb–Scargle periodogram によって周波数領域へ写し、変動構造として検定する選択肢を加えた[6]。
Drosophila の解析では、この違いが実際に検出される遺伝子の違いとして現れた。FBgn0040827、FBgn0031277、FBgn0030158、FBgn0262009 のように擬似時間の後半で発現が上昇し、そのまま高い状態で終わる遺伝子では、始点と終点の差が大きく残るため tradeSeq の startVsEndTest が強く反応しやすい。一方、FBgn0034659、FBgn0035569、FBgn0051406、FBgn0034739、FBgn0036488、FBgn0264344、FBgn0038978、FBgn0035585 のように中盤でピークを作った後に低下する遺伝子では、始点と終点の差が小さくても scLS が強く検出する例が示された[6]。同じ細胞群を測っていても、途中の発現変化を統計量へ残すかどうかで、候補遺伝子の順位が変わる。
T. brucei の ZC3H20 ノックアウトでも、同じ構造が別の形で現れた。ZC3H20 は実験的にノックアウトされた遺伝子なので、野生型との条件差が存在すること自体は既知である。それでも、軌跡全体を比較する scLS や patternTest、分岐初期の差を見る earlyDETest では非常に小さい P 値を示す一方、終点付近の差を見る diffEndTest では有意性が相対的に弱かった[6]。遺伝子もデータも同じままで、検定が見る区間と発現パターンの特徴を変えるだけで、統計的な強さが変わる。
この違いを「scLS は途中まで見ているから正しく、終点を比べる方法は情報を捨てている」と単純化すると、検定の役割を取り違える。終末状態で発現量が異なる遺伝子を知りたいなら、始点や終点に焦点を当てた検定は目的に合っている。分化の途中だけ一時的に発現する遺伝子を探したいなら、途中の形を保持する検定が必要になる。原論文が強調しているのも、すべての問いを一つの検定へ統合することではなく、どの種類の遺伝子を抽出したいかを先に定め、それに適した方法を選ぶことである[6]。
ここで第 2 章の tradeSeq に戻ると、associationTest、startVsEndTest、patternTest、diffEndTest が別々に存在する理由も同じところにある。それぞれは一つの「真の発現差」を異なる精度で推定しているのではなく、擬似時間上の異なる特徴について別の帰無仮説を検定している。scLS もその系列に一つの表現と検定を追加したのであり、従来の問いを無効にしたわけではない。解析法の選択は、計算アルゴリズムの選択であると同時に、どの発現変化を候補として前景化するかを決める工程になる。
既稿「判断が成立する世界はどのように設計されるか」では、現実をそのまま計算へ渡すことはできず、何を変数として取り出し、どの差を比較可能にし、何をモデルの外へ置くかによって、後段で成立する判断が変わる構造を論じた[14]。単一細胞の擬似時間解析では、この構造をかなり具体的に確認できる。RNA-seq で測定された発現値が同一でも、「初期と終盤の発現量」という二点へ縮約すれば中間のピークは検定量から消える。擬似時間上の曲線全体を残せば中間変化を扱えるが、曲線推定という別の仮定が入る。ピリオドグラムへ変換すれば不等間隔な擬似時間上の変動を別の特徴として扱える一方、差が生じた枝や位置を直接読む情報は失われる。どの表現にも、見えるものと見えにくくなるものがある。
既稿「意味は差異の読み取りから生まれる」では、物理的・生物学的な差異が存在することと、その差異をどの仕組みで読み取るかを区別した[15]。今回の例では、各細胞の RNA カウントという観測値がまず存在し、その後に擬似時間、始点と終点の差、回帰曲線、パワースペクトルといった異なる表現を通して差を読む。ただし、この読み取りは任意の解釈ではない。startVsEndTest には始点と終点の差という帰無仮説があり、scLS には Lomb–Scargle periodogram を使った検定統計量がある。それぞれシミュレーションや実データで挙動を検証でき、どの問いには強く、どこでは情報を失うかという技術的制約を持つ。
その制約があるからこそ、解析結果から生物学的な結論へ進むときには一段階を挟む必要がある。scLS が Drosophila で中盤ピーク型の遺伝子を強く検出したことから、「途中だけ発現する遺伝子も候補として拾える」とは言える。しかし、その遺伝子が分化を開始した原因かどうかは分からない。ZC3H20 のように実際のノックアウト介入がある場合と異なり、多くの候補遺伝子では発現変化が状態遷移の原因なのか、その結果なのか、さらに下流の応答なのかを区別する追加実験が必要になる。
冒頭の「低い → 高い → 低い」という例に戻る。初期と終盤だけなら、発現量はほぼ同じである。それでも中盤には大きな発現イベントが存在する。この曲線を二点へ縮約すればイベントは弱く見え、過程全体を残せば検出対象になる。scLS は、その過程を Lomb–Scargle periodogram によって周波数領域へ移し、一過性・非単調な発現を含む候補を拾う具体的な方法を示した[6]。
細胞状態の解析で「何が変わったか」を答えるには、まず何を変化として問うのかを決めなければならない。終末状態の差を知りたいのか、状態遷移の途中に現れる一過性発現を知りたいのか、条件によって発現の時期や形がずれる遺伝子を知りたいのかで、適した検定は変わる。scLS の意義は、同じ単一細胞 RNA-seq データから新しい事実を無条件に引き出すことではなく、これまで別の検定では弱くなりやすかった「過程の形」を検索対象として前に出せるようにしたことにある。
そして、その結果を原因遺伝子の確定と混同しないことが、この手法を正しく使う最後の条件になる。周波数領域への変換によって一過性の発現パターンを見つけ、候補を絞り、次に系譜ごとの解析、独立サンプルでの再現確認、ノックアウトやノックダウンなどの介入実験へ進む。状態遷移の途中で何が変わったかを広く拾うことと、その変化がなぜ起きたかを確定することを分ける。この二段階こそが、scLS の新規性を生かしながら、その統計的・生物学的な射程を越えない読み方になる。
参考文献
- Trapnell C, Cacchiarelli D, Grimsby J, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nature Biotechnology. 2014;32:381–386. https://doi.org/10.1038/nbt.2859
- Saelens W, Cannoodt R, Todorov H, Saeys Y. A comparison of single-cell trajectory inference methods. Nature Biotechnology. 2019;37:547–554. https://doi.org/10.1038/s41587-019-0071-9
- Street K, Risso D, Fletcher RB, et al. Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics. 2018;19:477. https://doi.org/10.1186/s12864-018-4772-0
- Van den Berge K, Roux de Bézieux H, Street K, et al. Trajectory-based differential expression analysis for single-cell sequencing data. Nature Communications. 2020;11:1201. https://doi.org/10.1038/s41467-020-14766-3
- Song D, Li JJ. PseudotimeDE: inference of differential gene expression along cell pseudotime with well-calibrated P-values from single-cell RNA sequencing data. Genome Biology. 2021;22:124. https://doi.org/10.1186/s13059-021-02341-y
- Iuchi H, Hamada M. The Lomb–Scargle periodogram-based differentially expressed gene detection along pseudotime. Nucleic Acids Research. 2026;54(13):gkag682. https://doi.org/10.1093/nar/gkag682
- Scargle JD. Studies in astronomical time series analysis. II. Statistical aspects of spectral analysis of unevenly spaced data. The Astrophysical Journal. 1982;263:835–853. https://doi.org/10.1086/160554
- Raz AA, Vida GS, Stern SR, et al. Emergent dynamics of adult stem cell lineages from single nucleus and single cell RNA-Seq of Drosophila testes. eLife. 2023;12:e82201. https://doi.org/10.7554/eLife.82201
- Briggs EM, Rojas F, McCulloch R, Matthews KR, Otto TD. Single-cell transcriptomic analysis of bloodstream Trypanosoma brucei reconstructs cell cycle progression and developmental quorum sensing. Nature Communications. 2021;12:5268. https://doi.org/10.1038/s41467-021-25607-2
- Zappia L, Phipson B, Oshlack A. Splatter: simulation of single-cell RNA sequencing data. Genome Biology. 2017;18:174. https://doi.org/10.1186/s13059-017-1305-0
- Papadopoulos N, Gonzalo PR, Söding J. PROSSTT: probabilistic simulation of single-cell RNA-seq data for complex differentiation processes. Bioinformatics. 2019;35:3517–3519. https://doi.org/10.1093/bioinformatics/btz078
- Hou W, Ji Z, Chen Z, Wherry EJ, Hicks SC, Ji H. A statistical framework for differential pseudotime analysis with multiple single-cell RNA-seq samples. Nature Communications. 2023;14:7286. https://doi.org/10.1038/s41467-023-42841-y
- Neufeld A, Gao LL, Popp J, Battle A, Witten D. Inference after latent variable estimation for single-cell RNA sequencing data. Biostatistics. 2024;25:270–287. https://doi.org/10.1093/biostatistics/kxac047
- id774, 判断が成立する世界はどのように設計されるか（2026-02-04）. https://blog.id774.net/entry/2026/02/04/3468/
- id774, 意味は差異の読み取りから生まれる（2026-05-09）. https://blog.id774.net/entry/2026/05/09/4740/