JSBi Bioinformatics Review
Online ISSN : 2435-7022
Review
Trends and outlook of single-cell analysis
Natsu Nakajima
Author information
JOURNAL OPEN ACCESS FULL-TEXT HTML

2022 Volume 3 Issue 2 Pages 61-74

Details
Abstract

生体組織のゲノム情報を個々の細胞レベルで網羅的に計測し解析するシングルセル解析が可能となり、組織を構成する各細胞の挙動や多様な細胞種の性質について分析することができる。シングルセル解析においては、各細胞での遺伝子発現量(Single-cell RNA sequencing)と数理的手法とに基づいて、細胞種分類や軌道推定、遺伝子間相互作用推定などが行われている。

本稿では、scRNA-seqデータに関する基礎知識や、シングルセル解析の動向と注目のトピックスに関して、各解析の概要と適用される数理的手法について概説する。

1.シングルセル遺伝子発現解析

人体の細胞数はおよそ数十兆個であるといわれており、生体組織や腫瘍組織も多様な種類の細胞から構成されている。各細胞では、同じ数の遺伝子が含まれているが、個々の細胞によって発現している遺伝子は異なるため、細胞の種類は、発現している遺伝子パターンにより特徴づけられる。Single-cell RNA sequencing (scRNA-seq)では、数千から数万の細胞について、各細胞での遺伝子発現量を網羅的に計測することができる。scRNA-seqの実験プロトコルとしては、マイクロウェルベース(well-based)であるSMART-Seq、CEL-Seq、MARS-Seq、SCRB-Seq (SMART-Seq2、CEL-Seq2、mcSCRB-Seq)、Quartz-Seq (Quartz-Seq2)、RamDA-Seq、マイクロ流体デバイスによるFluidigm C1、またドロップレットベース(droplet-based)でハイスループットであるDrop-seq、inDrop、Chromium 10Xなどが開発されている[1]。ウェルベースでは、数百細胞程度、ドロップレットベースでは、数千から数万細胞で、発現量はリード数やノイズを低減するUMI1カウント数で計測される。ドロップレットベースでは、細胞数は多く計測できるが、各細胞でのsequence depth2は低くなる。このようにscRNA-seqデータは、対象とする組織の細胞組成や酵素処理[2]などによる影響で、リード数やミトコンドリアRNA量などのデータのクオリティに差が出るため、シングルセル解析においては、各データの特徴によって正規化などの補正を行い、高次元の発現量データの可視化(次元削減)や細胞種分類(クラスタリング)、細胞分化の遷移推定(軌道推定)や遺伝子ネットワーク推定を行うことが重要である(図1)。また発現量の検出感度を向上させる補正(データ補完)については、任意ステップであるため、解析目的に応じた適用が有効である。解析に用いるデータとしては、個々の細胞での各遺伝子のmRNA発現量を表す行列である。本稿では、シングルセル解析について、各解析の概要と適用される数理あるいは情報科学的手法について概説する。2節では、scRNA-seqデータのクオリティを調整するデータの前処理について解説し、3節ではscRNA-seqデータからの特徴抽出と次元削減、4節で細胞種分類であるクラスタリングについて説明する。5節で細胞を擬似時系列に沿って並べる軌道推定と、6節では遺伝子間の相互作用を表す遺伝子ネットワーク推定について述べる。

図1:シングルセル解析の概要

シングルセル解析で用いられる、各細胞での遺伝子発現量行列(gene-cell expression matrix)は、細胞や遺伝子のフィルタリングや正規化を行い、任意のデータ補完やデータ統合により補正する。補正したデータをもとに、次元削減により二次元上に写して、クラスタリングや軌道解析、また遺伝子ネットワーク推定を行う。

2.前処理

2.1 フィルタリング

Droplet-based scRNA-seqでは、ハイスループットにより、数千から数万の細胞での遺伝子発現量を計測できるが、ミトコンドリア由来の発現量が高く転写産物の種類が少ない質の低い細胞や、単一のバーコードに複数の細胞で構成されるダブレットなども含まれ、下流解析に影響を及ぼす可能性がある[3, 4]。ダブレットでは、1つの液滴(droplet)に単一の細胞バーコードと2つ以上の細胞がカプセル化されてラベルされるため、ハイブリッドトランスクリプトームとなり、その割合は、スループットやプロトコルによるが、dropletの10%以下程度がダブレットとなっている[5]。ダブレットは、homotypicとheterotypicとに分類され、heterotypicは、異なる細胞や系統、状態により形成され遺伝子発現パターンも異なるために、特定しやすい[6]。ダブレットを特定するために、RNA発現量やマーカー遺伝子を用いた実験的、数理的手法が開発されているが、必ずしもダブレットでは総リード数が高い値になるとは限らないため、機械学習に基づく手法も開発されている。遺伝子発現パターンが異なる細胞のダブレットを特定する手法としては、DoubletFinder[7]、Scrublet[8]、cxds[9]、bcds[9]、DoubletDetection[10]があり、DoubletFinderとScrubletはartificial doublet (実際のdoubletを推定するために、同一または異なる細胞種の発現量を平均化して作成する擬似的なdoubletのこと)を用いてkNNグラフを構築して、ダブレットの割合として定義されるスコアを計算する手法である。cxdsは、2遺伝子が共発現している細胞のp値によりダブレットスコアを定義し、bcdsは、gradient boosted treesに基づいて確率を計算して分類する。一方で、同種の細胞により形成されるhomotypicダブレットの特定や、ダブレットスコアの閾値、人工的なダブレットを用いた2値分類問題の定義などの課題も含まれており、実験的、数理的手法の発展が期待される。

2.2 正規化

scRNA-seqデータは、シーケンスでの逆転写効率や各細胞でのリード数が低くなるために、多くの低発現遺伝子の発現量が0となる疎なデータとなる[11]。それらは、各細胞での遺伝子発現量のばらつきを減少させる影響があるため、発現量値を補正するために正規化を行う。一般的な正規化法としては、バルクのRNA-seqデータに対する手法も適用されており、総リード数に対する各サンプルに単一のスケールファクター(各サンプルでのリード数の比を調整する係数)を計算するが、scRNA-seqデータは検出感度が低いデータとなるために、一細胞の全遺伝子に対して単一のスケールファクターを用いた正規化では、遺伝子発現における特徴を捉えにくい[12]。

正規化には、グローバルスケーリングと非線形モデルに基づく手法があり、バルクデータに対しては、比較サンプル間での各遺伝子のリードカウントの比の幾何平均を揃えるTMM[13]やリード数分布の75パーセンタイルに比例するように計算するupper quartile normalization[14]、SCDE[15]、リード数の比の中央値からサイズファクターを計算するDESeq[16]などがある[17]。また0の発現量が多いscRNA-seqデータに対しても、頑健に各細胞に対するスケールファクターを計算するscran[18]、BASiCS[19, 20]、GRM[21]、各遺伝子に対してsequence depthと発現量の独立性により推定するSCnorm[12]、バルクデータと同様に、遺伝子のUMI値がべき乗分布に近似できることに基づき、大規模なデータにも拡張可能な正規化法であるPsiNorm[17]などが開発されている。正規化を行った発現量値は、分散安定化として対数変換を行い、各細胞での遺伝子発現量データとして、次元削減やクラスタリングに用いられる。このように、scRNA-seqにおける課題や特徴をふまえて、下流解析に影響するバイアスを補正した正規化法を適用することが有用である。

2.3 データ補完

scRNA-seqデータにおけるdropoutという現象は、シーケンスのプロトコルによるものであり、Fludigm C1では、細胞数は数百程度で、sequence depthは1-2万リードとなるが、dropletマイクロ流体デバイスでは、数万細胞でsequence depthは100-200kリード程度で計測される[22]。Droplet-basedの中でも、inDrop[23]や10X Genomics Chromium[24]では、転写産物の推定率が向上しているものの、比較的低感度で検出される。そこで、dropoutによる低検出感度を向上させるために、様々な補完手法が開発されている。基本的には、バルクやマイクロアレイデータに対する遺伝子発現量にzero-inflated negative binomial (ZINB)を用いてモデル化する方法が適用されるが、0が多いデータとなるため、細胞間、遺伝子間の類似度や色々なモデルに基づく補完手法が開発されている。

MAGIC[25]は、熱拡散の原理を応用し、類似した細胞間での発現量を補完する手法で、発現量データからマルコフ連鎖における遷移確率行列を計算して補完する。混合分布モデルに基づいた手法としては、SAVER[26]とscImpute[22]があり、SAVERは、Bayesの定理をもとに、Poisson-LASSOによる2遺伝子間の回帰により確率分布を計算し、同じ細胞の2遺伝子間の相互作用により補完する。scImputeでは、dropoutにより0となる遺伝子を推定するために、dropout確率を計算し、確率が高い遺伝子については、類似細胞の同じ遺伝子の発現量を用いて、回帰係数が正値である最小二乗法(非負値最小二乗法)により補完される。

さらに、発現量が0の場合に、細胞種特異的に発現していない(true zero)か、dropoutによるもの(false zero)かを推定して、dropoutによるzeroを補完する手法も開発されている。PRIME[27]とkNN-smoothing[28]は、細胞間類似度に基づいた手法であるが、PRIMEでは、細胞間のネットワークを推定して、遺伝子がある細胞種で発現していて、他細胞種で発現していない場合にはdropoutによる0であると仮定している。kNN-smoothingは、k-nearest neighborにより、細胞間類似度から隣接細胞群を推定して、それらの発現量を合算することで補完する。

また、deep learningによる大規模なデータにも拡張可能な補完手法も開発されており、DCA[29]は、ニューラルネットワークのオートエンコーダを用いた手法で、各遺伝子に対してZINB回帰により、誤差関数を最小化するような分布のパラメータである平均、分散、dropout確率を計算して発現量を補完する。DeepImpute[30]は、ディープニューラルネットワークと分割統治法に基づき、遺伝子発現データを分割して、複数の部分ネットワークを構築し、dropout率などを推定して発現量を補完する。このように、クラスタリングや低発現遺伝子の発現を考慮に入れた解析を行う場合には、データ補完手法の適用を検討すべきであるが、任意のステップであり、補完によるバイアスがかかる場合もあるため、データにおけるスパース率や解析目的に応じた適用を考慮することが有効であると考えられる。

2.4 データ統合とバッチエフェクト

シーケンス技術の進歩により、シングルセルレベルでもmRNA発現量の他に、免疫形質[31, 32]、ゲノムシーケンス[33, 34]、DNAメチル化[35, 36]、クロマチン近接可能性[37, 38, 39]、空間的位置[40, 41, 42]などの異なる細胞のモダリティを計測できるようになっている[43]。例えば、scATAC-seqは、エンハンサー領域や発現制御機構に関する情報を有するが、細胞種分類では、mRNA発現量も加味することで情報量が多くなる。またSlide-seqによる空間的なmRNA発現量は、組織構成を捉えることができるが、遺伝子全体の発現量は計測されにくい。従って、scRNA-seqと空間的情報やscATAC-seqなどの異なるモダリティデータを統合して解析することで、細胞種分類のみならず、発現制御機構や空間的構成に関する理解へと繋がることが期待できる[43]。また異生物種や刺激応答データに関しても、統合解析を行うことで、それらに共通する細胞種を分類することができる。

異なるモダリティデータの統合に関しては、正準相関分析(CCA)と複数のデータ間で次元削減を行い、共通する多様体(データが、見かけの次元より十分少ないパラメータで記述できる曲がった空間に収まっていること)を学習する機械学習法(多様体アラインメント)により、異なるデータ間で共通する遺伝子間の相関を推定し、低次元空間に写像する手法[44]や、mutual nearest neighbors (MNNs)を用いて、データ間に共通の部分的な細胞群を推定する手法[45]、確率的最適化とニューラルネットワークに基づいて、類似した細胞群を推定する手法[46]などが開発されている。mRNA、エピゲノム、プロテオミクス、空間的なデータに対しては、訓練されたモデルを他の関連する課題に応用する転移学習において、訓練データとテストデータの確率分布が異なる場合の学習法(ドメイン適応)を用いて、データ間に共通の隣接細胞群から変換行列を計算してバッチエフェクトを補正し統合する[43]。バッチエフェクトとは、各データが実験ごとに異なる細胞状態やシーケンスプラットフォーム、実験プロトコルにより測定されるため、実験条件の違いが要因のデータ間の差異を指す。Human Cell Atlas (HCA)でも、異なるモダリティデータが蓄積されており、数百万細胞からなり異なる実験環境で測定されているため、組織構成の差異や細胞種分類を行うためには、バッチエフェクトを補正する方が良い。バッチエフェクト補正法としては、MNNsを用いて各データに含まれる類似する細胞種や状態の細胞群を推定し、それらの発現量の差異がバッチエェクトであると定義して、ガウシアン平滑化により補正する方法[45]や、発現量の中央値を用いて、各細胞をリファレンスデータでのクラスタや類似細胞に写像するscmap[47]などがある。またbayNorm[48]は、ベイズ統計の手法をもとに、分布パラメータを計算し、正規化やデータ補完、バッチエフェクト補正を統合的に行う。

シングルセルデータは、実験環境や個人差による影響も大きく、クラスタリングでは細胞種ごとに分類されないことも多いため、異なるデータの統合解析を行う際には、バッチエフェクトの補正も検討した方が効率的であると考えられる。

3.次元削減

次元とは、データを特徴づける特徴量数のことであり、scRNA-seqデータによる細胞種分類を行う際には、遺伝子が特徴量となる。scRNA-seqデータは、各細胞での数万程度の遺伝子の発現量を表すため、数万次元となる高次元データである。生物学的プロセスは多くの遺伝子により構成されているが、各細胞種を特徴づける遺伝子はそれほど多くはなく、発現量データの2値化を行うことでも、シングルセルデータの本質的な構造を捉えることは可能であると報告されている。例えば、造血幹細胞の分化では、2次元か数次元で、データが持つ特徴は捉えることが可能で、ある次元は、細胞種にどれくらいの経路で分化するかを表し、その他の次元は、各細胞種での細胞周期を表している[49]。このように、少ない次元数にすることでデータの特徴が捉えやすくなるため、データが持つ位相幾何学的な特徴を可視化することが可能である。ゆえに、次元削減とは高次元データの構造をできるだけ維持して低次元に写すことであり、細胞のクラスタ分類や分化軌道推定の可視化を行うためには必須のステップである。

高次元データを低次元に削減する際に、特徴選択と特徴抽出を行うが、特徴選択は、高次元のデータを最も良く説明する特徴量を選択することであり、特徴抽出は、高次元データの特徴量と組み合わせて、データを要約する新たな特徴量を抽出することである。特徴選択法には、Filter法、Wrapper法、Embedded法があり、各変数について分散などを計算し、変数を組み合わせて学習することで選択を行う[50]。発現変動遺伝子も、細胞間での発現量の変動が大きい遺伝子群で、特徴量として抽出される。発現変動遺伝子を計算する手法では、発現量の分散と平均の関連をフィッティングして、発現量の変動が大きい遺伝子について統計テストを行う[51]。シングルセルデータの次元削減における特徴抽出法としては、線形変換のPCA[52]、非線形変換のt-SNE[53]、UMAP[54]、deep learningやカーネル法に基づくDCA[29]やSIMLR[55]が代表的である(図2)。

図2:scRNA-seqデータを用いた次元削減による二次元への写像

膵臓のscRNA-seqデータ[57]を用いて、PCAで20次元まで次元削減し、t-SNEとUMAPにより2次元に写す。(A)t-SNE(B)UMAP

―PCA

PCAはデータの分散を最大化するような変数の線形結合により、分散が最大となる軸を主成分として抽出する。数万次元から数十次元までPCAで次元削減を行ってから、t-SNEやUMAPにより2次元に写す場合が多い。

―t-SNE

t-distributed Stochastic Neighbor Embedding (t-SNE)は、SNE[56]に基づいて、高次元データのローカル構造(近傍点の配置)を捉えるように、非線形変換によって次元削減を行う。SNEでは、高次元と2次元において、ユークリッド距離を用いて2細胞間の距離を条件付き確率とした確率分布で表現し、低次元空間での確率分布とのKLダイバージェンスを最小化するように(高次元で近くの点は、低次元でも近くに配置)2次元に写す。t-SNEでは、距離の対称性とt分布を用いて、低次元空間でのデータ点の混雑を軽減して配置するように写像する。

―UMAP

Uniform Manifold Approximation and Projection (UMAP)は、リーマン多様体学習に基づいて、高次元データのグローバル構造を維持するように、非線形変換を行う。多様体とは、データが、見かけの次元より十分少ないパラメータで記述できる曲がった空間に収まっていることを意味する。UMAPでは、高次元で細胞間類似度を表すkNNグラフを作成し、グラフの辺がある部分についてのみ2細胞間の確率分布を計算して、辺の密度が高い部分を近くの点として配置するように2次元に写す。UMAPは、辺の密度が高い部分グラフに対して変換を行うため計算時間を少なくすることができ、写像空間の次元数が高くても適用することができる。

―SIMLR

SIMLRは、マルチカーネル学習に基づいた手法で、2細胞間の類似度を用いて、対角成分として細胞種数のブロック構造を持った対称行列を計算する。類似度を表す親和性行列をインプットとして、t-SNEを用いて次元削減を行う。

このように、次元削減による細胞種分類や高次元データ構造の可視化については、少ない遺伝子の発現パターンにより特徴づけられ、それらは特徴量として抽出される。一方で、発現量が0の遺伝子や低発現遺伝子などは抽出されない傾向があるため、個々の細胞での遺伝子発現の差異までは捉えにくいが、例えば、細胞ラベルなしのデータを用いた細胞のクラスタ分類やマーカー遺伝子発現差解析などを行う際には、次元削減によってデータのおよその特徴を捉えることが可能である。

4.クラスタリング

多細胞生物は、分化した細胞群である組織により構成され、1つの組織が複数の細胞種から構成されていたり(細胞不均一性)、同じ細胞種でも同一の刺激に対する個々の細胞での遺伝子発現は異なる。ゆえに、scRNA-seqデータを用いたクラスタリングを行うことによって、組織における細胞種の多様性や各細胞種に特異的なマーカー遺伝子の推定、時空間情報に基づいた分化における異なる細胞系統の軌道を辿ることができる[58]。クラスタリングにおいても、遺伝子を発現パターンによって分類することは重要であるが、dropoutや技術的誤差による影響で疎なscRNA-seqデータでは、発現変動遺伝子の推定や、類似した細胞種なども分類されにくい傾向がある。

シングルセル解析でのクラスタリングは、古典的なクラスタリング法やscRNA-seqに対する手法、コミュニティ抽出法を適用して行われる。古典的な手法では、k-means法[59]、spectral clustering[60]、EM法[61]があり、次元削減により写した2次元座標に対して適用する。k-means法ではクラスタ数を定義し、spectral clusteringでは、ラプラシアン行列の固有値分解によりクラスタリングを行う。

scRNA-seqデータに対する手法としては、sc3[62]は、コンセンサスクラスタリングに基づいて、ユークリッド距離、ピアソン、スピアマン相関により、2細胞間類似度を表す行列を作成し、各行列について、PCAまたはラプラシアン固有マップによる次元削減を行って、変換した行列の各固有ベクトルについて、k-meansクラスタリングを行う。各クラスタリングにより、2細胞が同じクラスタに含まれているかを2値の類似度で表すコンセンサス行列を作成して階層的クラスタリングを行う。SINCERA[63]は、主要な細胞種分類や各細胞種でのマーカー遺伝子推定、分化を制御する転写制御ネットワーク推定を行うパイプラインで、ピアソン相関による類似度をもとに、クラスタ分類を行い、各クラスタについて、ウィルコクソン検定により発現変動遺伝子の推定を行って細胞種を推定する。またscRNA-seqでのスパース性や異なる遺伝子が同じ機能を有していたり、1つの遺伝子が複数機能を司っているという性質にテキストマイニングを応用して、LCA[64]では、情報検索におけるトピックモデルの潜在的意味解析(Latent semantic indexing (LSI))に基づいた手法を開発している。特異値分解により次元削減を行い、二重の低次元空間であるLCの次元空間では2細胞間のコサイン類似度を計算して、PCAによる次元空間で、2つのクラスタ間での非類似度をシルエット係数により計算して、細胞種分類のクラスタ数を推定する。

また特定の細胞状態や疾患の細胞種において、同一の転写制御シグナルに対して機能的に関連する、もしくは共発現している状態特異的な機能的遺伝子群(condition-specific functional gene modules (FGMs))を推定することも応用されている。つまり、FGMsの推定を行うことで、マーカー遺伝子の推定、さらには、刺激応答としての個々の細胞間相互作用についての理解に繋がると考えられる[65]。QUBIC2[65]は、グラフ理論に基づいて、2次元のデータマイニングにより、ある細胞状態で共発現する遺伝子群を同時に推定するバイクラスタリング(biclustering)を行う手法である。発現量の2値化を行い、ある状態で同一のシグナルに応答する遺伝子をバイクラスタとし、そうでない遺伝子は異なるシグナルに応答するとして、それらの遺伝子とバイクラスタ間の類似度をKLダイバージェンスにより計算してコアバイクラスタを推定する。また水平、垂直方向にもクラスタリングを行うことでデュアルバイクラスタを計算している。

さらに、ネットワークのコミュニティ抽出法であるLouvain法[66]やLeiden法[67]も、クラスタ数によらずモジュラリティを最大化するように分類を行うため、シングルセル解析においても適用されることが多い。細胞ラベルなしのデータの細胞種分類においては、次元削減後に、コミュニティ抽出法や古典的手法も適用することが有用であるため、マーカー遺伝子推定や分化軌道解析などの下流解析の目的に応じたクラスタリング手法の適用を検討されたい。

5.軌道解析

scRNA-seqは、個々の細胞での遺伝子発現量を表すと同時に、細胞分化における分化の段階や進行度についての情報も与える。mRNA、プロテオミクス、エピゲノミクスデータなどを含むシングルセルオミクスデータは、細胞周期や細胞分化、細胞活性などの細胞状態の遷移を捉えることが可能であり、そのようなダイナミクスは、軌道推定である擬似時系列解析によってモデル化される。軌道推定とは、遺伝子の発現パターンによる細胞類似度をもとに、軌道に沿って細胞を並び替えることで、70以上の推定手法が開発されている[68]。軌道の形状としては、linear、bifurcation、multifurcation、treeなどが多く、cycle、connected graph、disconnected graphなどの複雑な形状にも対応して推定される。推定手法には、形状は固定して細胞を並び替える手法[69, 70, 71, 72]や近年では、枝における細胞の配列順と枝に接続する形状の両方を推定して頑健に軌道を推定する手法[73, 74, 75, 76]も開発されている。Wanderlust[70]では、始点の細胞も入力として、kNNで次元削減を行い細胞間の最短距離によって並び替える。Monocle[77]では、独立成分分析(independent component analysis (ICA))により次元削減を行い、各細胞状態での高発現変動遺伝子を推定して、それらの細胞での最小全域木を計算し、パスの最近接点に各細胞を写像する。SCORPIUS[78]では、全細胞ペア間についてスピアマン相関で距離を計算し、多次元尺度構成法(multi dimensional scaling (MDS))により次元削減を行い、k-means法により細胞のクラスタリングを行って、各クラスタの中心からの最短距離を計算することで初期軌道を推定する。

―擬似時系列解析

擬似時系列解析とは、横断的に計測された高次元データを、擬似時系列という1次元の特徴として捉えることであり、分化や疾患などの時系列の挙動を要約することができる。擬似時系列とは、疾患の進行度や細胞の分化での細胞状態の遷移を理解するため、発現量データをもとに、個々の細胞間の類似度に基づいて、擬似的な時間軸に沿って並び替えることである。推定手法としては、発現パターンが類似する2細胞は、それに応じて類似した軌道上にあるようなグローバル構造を保持して各細胞を並び替えるが、その類似度の定義が各手法により異なっている。PhenoPath[79]では、異なる刺激に対して個々の細胞がどのように応答するかなどの遺伝的、表現型的、環境的なファクター間の関連も含めたベイジアンモデルに基づく確率論的な推定手法を開発している。線形回帰と擬似時系列の進行度を調節する潜在変数モデルとにより、分化や疾患の進行に沿って相対的な位置に配置されるような軸を推定する。Diffusion pseudo times (DPT)[71]では、kNNグラフと遷移確率により細胞間類似度を計算し、Monocle2[76]では、次元削減を行って、reverse graph embedding (RGE)により、枝の分岐点を教師なし学習で推定する。モデルに基づく分岐点の推定としては、DeLorean[80]やSlingshot[81]があり、各クラスタで最小全域木により軌道のグローバルな構造を推定し、複数の分岐系において各系統でのprinciple curve (主曲線)を同時にフィッティングする手法(simultaneous principal curves)により各細胞系統での擬似時系列を推定する。CellTrails[82]では、ファジイ相互情報量によりグラフを作成し、固有値分解により低次元でのデータ構造を捉えて、spectral clusteringで各クラスタでの細胞間類似度を計算する。kNNにより各細胞状態での類似細胞を推定し2細胞間の近接度(ジオメトリ近接)と各状態での類似細胞の割合(カーディナリティ)を最大化するような状態木を推定する。各細胞は距離が最も近くなるフィッティング直線上に直角に写し、直線に沿って並べ分岐軌道を推定する。

―RNA velocity

また分化や組織再生は、時間的な尺度に基づくダイナミクスであるため、遺伝子スプライシングなどの速度を計算して軌道を推定することができる[83]。RNA velocityは、遺伝子発現の(状態)空間における、個々の細胞の状態変化の向きや速度(velocity)を表し、分岐軌道の推定や各細胞のダイナミックな遷移を写して二次元上に可視化される[84]。RNA velocityの推定については、遺伝子のスプライシングの特徴を用いて、遺伝子特異的な定常状態の平衡定数を分位点回帰により推定し、推定されたvelocityは、低次元空間での隣接細胞と他細胞との類似度により写像する[83]。VeTra[84]では、RNA velocityをもとに有向グラフを作成し、TENET[85]により軌道における同じ細胞系に属する各細胞群での主要な制御因子を推定する。有向グラフは、低次元での類似細胞のスプライシングリード数によりvelocityベクトルを計算し、コサイン類似度により同じ方向の類似細胞を選択して作成する。グラフから連結成分を推定し、近接した類似細胞に分類して軌道を推定する。CellPath[86]は、発現量データとRNA velocityにより、k-means法とLeiden法でメタ細胞群を推定し、それらでのkNNグラフを作成する。ダイクストラ法により多くの頂点が含まれるようにパスを推定し、メタ細胞群の中央を通って各細胞を並べ替えて軌道を推定する。

このように、軌道推定手法は、細胞間類似度やRNA velocityをもとにした、様々な分岐の形状に対応した手法が開発されているため、データの系に適した推定手法の適用を検討するのが良い。

6.遺伝子ネットワーク推定

6.1 遺伝子相互作用ネットワーク推定

遺伝子ネットワークは、遺伝子を表す頂点と遺伝子間の相互作用を表す辺により構成されており、生物学的プロセスに影響を与える制御遺伝子の相互作用なども含まれている。ダイナミクスにおける相互作用ネットワークは、各発生段階において、異なる細胞種での、時空間的な遺伝子発現を制御しているため、それらの遺伝子の機能や相互作用を理解することは、生物学的プロセスを解釈する上でも重要となる。

遺伝子ネットワークは、バルクデータに対する推定手法[87, 88, 89, 90, 91, 92, 93]や転写制御の相互作用[94, 95, 96, 97]を推定する手法が数多く開発されている[98]。scRNA-seqデータに対しても、正常組織と疾患組織での細胞種の差異に関連する個々の細胞状態を制御する転写制御因子を推定する手法[99, 100, 101, 102, 103]が開発されている[98]。推定手法は、遺伝子ネットワーク推定の仮定によって異なり、入力は発現量や発現量を2値化したデータ、擬似時系列や遺伝子間の相関によって変換したデータである。遺伝子間の相互作用は、連続値や2値の隣接行列、スコアによる2頂点間の辺またはブール関数として出力され、頂点と辺から構成されるグラフであるネットワークを作成する。遺伝子間の制御規則は、ブーリアンモデル、微分方程式、相関に基づいたモデル化が代表的である。ブーリアンモデルでは、発現量を2値化して、遺伝子間の制御はブール関数によって表され、関数により各頂点の状態を最適化する。微分方程式では、各遺伝子の発現量の時間変化量を、他の遺伝子の発現量の関数で回帰してモデル化を行う。また2遺伝子間の相関は、ベイジアンモデルや相互情報量により計算される[98]。ブーリアンモデルでは、SCNS[104]やBTR[105]があり、SCNSは、各細胞状態の遷移の規則であるブール関数については、状態遷移グラフを作成し、各細胞の状態空間に沿った関数を推定する。BTRでは、発現量データから状態空間でのブーリアンモデルを作成して、距離スコアを計算し、各細胞の状態を推定する。GRISLI[106]とSCODE[107]は、微分方程式に基づいて、GRISLIは、各細胞において、各遺伝子の発現の増加や減少の速度を推定し、転写因子の遺伝子発現量を、各細胞の速度行列と制御因子の発現量行列から計算し、ラッソ(lasso)回帰により相互作用を推定する。速度は各細胞が軌道上のどの位置にいるかの時刻ラベルをもとに、2細胞間の距離を計算し、それらの加重平均として推定する。SCODEでは、各細胞について、軌道上での時刻を計算して、各転写因子と全遺伝子の発現量から線形回帰により発現量の二乗誤差を最小化して相互作用の係数行列を計算し、転写因子間の相互作用を推定する。

このように、scRNA-seqデータを用いたネットワーク推定としては、各細胞種で特異的に発現する転写因子に関する相互作用や、細胞状態の遷移から遺伝子間相互作用の推定を行う手法が多く開発されている。またバルクデータに対する手法も、サンプル数を細胞数として発現量データを入力とすることで適用することができる。ネットワーク推定では、遺伝子間の相互作用の全体像を捉えることができるため、細胞数が多いscRNA-seqデータに適用することで、各細胞種を特徴づける遺伝子やその機能の理解への発展が期待される。

6.2 遺伝子共発現ネットワーク推定

生物学的プロセスにおける遺伝子の機能や疾患に関連する遺伝子を特定するためには、遺伝子共発現ネットワークを構築して解析することが重要である。共発現ネットワークは、サンプル間を通して、同時に活性化したり、発現パターンが類似した遺伝子間に辺を張ることで作成され、さらにネットワークのクラスタリングを行うことでモジュールとして分類される[108]。共発現遺伝子のモジュール分類によって、機能的に関連のある遺伝子群や、疾患に関連するマーカー遺伝子、ハブとしての機能を持つ遺伝子などを推定するため、遺伝子の機能について概観することができる。

サンプル間を通しての発現パターンの類似性は、ピアソン相関係数[109]や相互情報量[110, 111]により計算され、共発現する遺伝子から構成されるネットワークを推定し、発現が類似する遺伝子群をモジュールとして分類する[107]。モジュールはk-means法やLouvain法などで分類され、機能的に関連のある遺伝子群の抽出や、differential co-expressionでは、異なる疾患や組織での共発現遺伝子群を推定する。また、各モジュールの機能を制御する遺伝子であるハブ遺伝子も、媒介中心性により計算され、疾患に関連する遺伝子群などを推定することができる。scRNA-seqデータに対する推定手法として、LEAP[112]は、細胞の擬似時系列を推定して、各時間区間での2遺伝子間の共発現度を表す最大絶対相関(maximum absolute correlation)に基づき、p値を計算する。scLink[113]では、2遺伝子間の共発現行列を計算して、共発現度をもとに正則化項の割合を計算する手法(データアダプティブ尤度法)によりスパースなネットワークを推定する。また、疎なscRNA-seqデータから細胞集団特異的な遺伝子群を特定するため、Co-Dependency Index[114]に基づいて、遺伝子共発現ネットワークを構築し、共発現遺伝子群をモジュールとして分類することで機能的に関連のある遺伝子群を推定する手法も開発されている[115]。さらに、相互排他的発現の遺伝子を推定する指標、Exclusively Expressed Indexで計算される遺伝子群を細胞分類の特徴量として、クラスタリングに応用する手法も開発している。

scRNA-seqデータからの共発現遺伝子ネットワーク推定やモジュール抽出によって、遺伝子の機能的分類を行うことができるため、各細胞種でのマーカー遺伝子や機能的特徴に関する知見をより深めることができる。

7.まとめ

シングルセル解析では、単一細胞レベルでの遺伝子発現量をもとに、個々の細胞の個性や遺伝子発現の多様性を理解することができる。また、計測技術の改良や解析のための数理的手法も数多く開発されているため、シングルセルデータの特徴と解析目的とに応じた数理的手法の適用が重要である。また、疾患に関する臨床データの解析や身体を一細胞という解像度で地図化するHuman Cell Atlas (HCA)、複数のモダリティデータの統合などは、生物学的な知見を深め、数理的手法についても更なる発展が見込まれる。特に、HCAデータは大規模なデータとなるため、ユーザーに使いやすく、スケーラブルといった手法の有用性を高くすることである。またツールとしては、各解析に特化している場合が多いため、ツール間での互換性を高くすることで、解析パイプラインの構築や統合的解析も可能となると考えられる。シングルセルデータの有する情報を、細胞から組織へと統合的な解析を行うことで、様々な生命現象や疾患におけるメカニズムの解明が期待される。

脚注

1 UMI (Unique Molecular Identifiers):各RNA分子に付加された分子バーコード(塩基配列)で、PCR増幅によらない各細胞でのRNAのコピー数を正確に測定することができる。

2 各細胞でのリード数

References
著者略歴

仲嶋 なつ
2015年に京都大学大学院情報学研究科博士後期課程修了し、同年に京都大学化学研究所バイオインフォマティクスセンター研究員を経て、2017年から2021年5月まで東京大学定量生命科学研究所特任研究員としてシングルセル解析研究に従事。遺伝子ネットワーク推定手法やデータ解析手法の開発を行っている。

 
© 2022 Japan Society for Bioinformatics

This article is licensed under a Creative Commons [Attribution-NonCommercial-ShareAlike 4.0 International] license.
https://creativecommons.org/licenses/by-nc-sa/4.0/
feedback
Top