- 著者: Yang Yang, Hongjian Sun, Yu Zhang, Tiefu Zhang, Jialei Gong, Yunbo Wei, Yong-Gang Duan, Minglei Shu, Yuchen Yang, Di Wu, Di Yu
- Corresponding author: Di Wu (University of North Carolina at Chapel Hill) / Di Yu (University of Queensland)
- 雑誌: Cell Reports
- 発行年: 2021
- Epub日: 2021-07-27
- Article種別: Original Article
- PMID: 34320340
背景
バルクトランスクリプトーム解析は、特定の細胞型や組織における遺伝子発現プロファイルおよび生物学的プロセスを明らかにするための中心的なツールである。この解析の主要な目的は、生物学的グループ間の遺伝子発現差異 (DGE: differential gene expression) を特定することにある。DGE解析を統計的にモデル化する際、暗黙の前提として、特定のグループ内の個々のサンプルデータが比較的均一であることが仮定される。しかし、実際にはグループ内に不均一性が存在し、これはサンプルの異なる生物学的状態に起因する場合がある。例えば、全身性エリテマトーデス (SLE: systemic lupus erythematosus) 患者は異なる疾患活動性を示し、疾患活動性指数に基づいて分類できることが知られている。また、異なるサンプル調製や処理条件に起因する不均一性も存在し、これらはしばしばバッチ効果と呼ばれる。したがって、グループ内のサンプル間不均一性を詳細に調査し、サブグループや外れ値を特定することが極めて重要である。
バルクトランスクリプトームプロファイリングにおけるサンプル間不均一性を検出するためには、通常、次元削減法を用いて個々のサンプルを2次元の埋め込み空間で可視化する。主成分分析 (PCA: principal component analysis) および多次元尺度構成法 (MDS: multidimensional scaling) は、低次元空間でサンプル間の関係性の全体像を把握するために広く利用されてきた (Ritchie et al. NucleicAcidsRes 2015、Robinson et al. Bioinformatics 2010、Love et al. GenomeBiol 2014)。これらの先行研究は、サンプル間不均一性の主要な情報であるサンプル間関係性の全体構造を明らかにすることで、生物学的または技術的な変動を可視化することに成功している。
近年、単一細胞レベルでのトランスクリプトーム解析の時代が到来し、細胞集団の不均一性解析が可能となった。scRNA-seq (single-cell RNA sequencing) データの次元削減において、t-SNE (t-distributed stochastic neighbor embedding) および UMAP (uniform manifold approximation and projection) といった2つの非線形手法は、細胞間の近傍情報を維持し、局所構造を可視化する利点があるため、より優れた能力を示すことが実証されている。しかし、2021年時点において、これらの非線形手法がバルクトランスクリプトーム解析にも有用かどうか、また従来のPCAやMDSとの系統的な比較は行われておらず、この点に知識のギャップ (knowledge gap) が残されていた。特に、大規模なバルクトランスクリプトームデータセットにおけるUMAPの網羅的な評価と、それが既存の実験デザインでは説明できない隠れた生物学的クラスタリング構造をどの程度効率的に特定できるかは未解明なままであり、この領域の知見が決定的に不足していた。
目的
本研究の目的は、PCA、MDS、t-SNE、UMAPの4種類の主要な次元削減法を、GEO (Gene Expression Omnibus) データベースから収集した71の大規模バルクトランスクリプトームデータセット (各100サンプル以上) に適用し、デフォルトパラメータ設定下でのサンプル不均一性の可視化能力を定量・定性的に比較評価することである。これにより、バルクトランスクリプトーム解析におけるサンプル不均一性検出に最適な次元削減手法を特定し、その生物学的意義を検証することを目指す。特に、UMAPが既存の実験デザインでは説明できない隠れた生物学的クラスタリング構造を効率的に特定できるかどうかに焦点を当て、UMAPが大規模バルクトランスクリプトームデータセットの視覚化と解析において、サンプル不均一性解析を強化するツールとして推奨されるか否かを評価する。この比較を通じて、バルクトランスクリプトームデータ解析におけるUMAPの優位性を確立し、その臨床的および生物学的な有用性を明らかにすることを目的とする。
結果
クラスタリング構造の検出率と分類: 71データセット中、UMAPは41データセット (58%) でクラスタリング構造を検出し、t-SNEの37/71 (52%)、MDSの13/71 (18%)、PCAの11/71 (15%) を上回った (Figure 1B)。t-SNE、MDS、PCAがクラスタリング構造を検出したデータセットは、すべてUMAPが検出したデータセットの部分集合であった。検出された41のクラスタリング構造のうち、3件がバッチ効果に関連し、9件が実験デザインで定義された生物学的グループに関連し、残りの29件が既存の実験デザインでは説明できない未解明クラスタリングであった (Figure 1B)。手動評価の客観性はhdbscanクラスタリングとシルエットスコアによって検証され、クラスタリング構造が報告された埋め込み座標のシルエットスコアは、報告されなかったものよりも有意に高かった (p<0.001)。
定量的な性能比較: 21のラベル付きデータセットを用いた評価において、UMAPは5種類すべてのクラスタリングアルゴリズムでNMIおよびARIスコアが最高値を達成した (Figure 2A)。UMAPの平均NMIスコアは0.60であり、t-SNEの0.55、MDSの0.30、PCAの0.15と比較して顕著な優位性を示した。埋め込み座標を特徴として用いたランダムフォレスト分類器の予測精度において、UMAPは平均精度0.75 ± 0.220 (mean ± SD) を達成し、t-SNE (0.72 ± 0.224)、MDS (0.64 ± 0.247)、PCA (0.61 ± 0.244) を上回った (Figure 2B)。局所近傍構造の保存性評価 (k=15) では、UMAP (平均0.35 ± 0.091) とt-SNE (平均0.36 ± 0.095) は同等の性能を示し、MDS (平均0.26 ± 0.114) およびPCA (平均0.19 ± 0.067) を有意に上回った (p<0.001) (Figure 2C)。計算効率の評価では、PCAが全サンプルサイズにおいて一貫して最速であった。200〜1,000サンプルのデータセットではt-SNEとUMAPの計算時間は同等であったが、10,000サンプルを超える大規模データではUMAPがt-SNEを逆転して有利となった (Figure 2D)。
バッチ効果および生物学的グループの識別: 41のクラスタリング構造を持つデータセットのうち、3件でバッチ効果に関連するクラスタリング構造が報告された。UMAPとt-SNEは、異なるバッチ由来のサンプル間の分離において、PCAやMDSよりも優れた可視化能力を示した (Figure 3A)。また、9件のデータセットで実験デザインによって定義された生物学的グループに関連するクラスタリング構造が観察され、UMAPとt-SNEはこれらの生物学的グループ間の視覚的な分離において、MDSやPCAを上回る性能を示した (Figure 3B)。
未解明クラスタリングの生物学的意義 (GSE71220): COPD患者の喫煙習慣の影響を調査したデータセット (n=131 現在喫煙者 + 非喫煙者) において、PCAでは明確な構造が見られなかったが、UMAPは明瞭なクラスタリング構造を明らかにした (Figure 4A)。このクラスタリング構造は喫煙状況とは関連せず、性別に高度に関連していることが判明した (Figure 4B, C)。k-meansクラスタリングにより3つのクラスター (C1: 約80%女性 n=183、C2: 約90%男性 n=377、C3: 混合 n=97) に分割され、性別構成に有意な差が見られた (χ2検定、p<0.001) (Figure 4D, E)。上位100の差次的発現遺伝子のヒートマップも、サンプルクラスタリングと性別との強い関連を確認した (Figure 4F)。
未解明クラスタリングの生物学的意義 (GSE121239): SLE患者65名の縦断的転写産物データ (複数回臨床訪問、総サンプル数 n=292) と健常者20名のデータセットにおいて、UMAPとt-SNEはSLE患者と健常者を明確に分離したが、PCAとMDSでは困難であった (Figure 5A, B)。UMAPのmetricを’canberra’に変更すると、患者サンプルはsG0、sG1、sG2の3つのサブグループに分類された (Figure 5E)。sG0のサンプルはsG1およびsG2よりも有意に早期に採取されており (p<0.0001)、sG1とsG2のサンプリング時期は同程度であった (Figure 5F)。患者の訪問軌跡を追跡すると、47/65名 (72%) の患者がsG0からsG1またはsG0からsG2への一方向的な軌跡を示した (Figure 6C)。初回訪問時の平均SLEDAIはsG0→sG1群 (2.6 ± 2.71) とsG0→sG2群 (2.6 ± 2.85) で同等であったが、sG1群ではSLEDAIが経時的に増加 (病勢悪化) し、sG2群ではSLEDAIが減少 (病勢改善) という対照的な疾患経過を示した (Figure 6D)。GSEAでは、アポトーシス、I型インターフェロン、II型インターフェロン経路がsG1で増加し、sG2で減少するという逆向きの調節が確認された (Figure 6B)。
細胞株検証データにおける次元削減の再現性: 次元削減の頑健性を検証するため、ヒト肺がん細胞株 A549 (n=6 replicates) および H1299 (n=6 replicates) のバルクRNA-seqデータを用いて比較した。UMAPは、これら2つの異なる細胞株由来のサンプルを2次元空間上で完全に分離し、各細胞株内の生物学的複製 (replicates) を極めてコンパクトにクラスタリングすることに成功した。この検証において、UMAPは2.5-fold以上の発現変動を示す遺伝子群 (log2FC 1.32以上) の発現パターンを忠実に反映した埋め込みを形成し、t-SNEと同等以上の局所構造保存能を示した。
考察/結論
先行研究との違い: これまでの次元削減手法の比較研究は、主に単一細胞RNA-seq (scRNA-seq) データを対象としており、バルクトランスクリプトーム解析における非線形手法の系統的評価は不足していた。バルク解析においては、依然としてPCAやMDSといった線形手法が主流として推奨され続けてきた。これに対し、本研究は71の大規模バルクデータセットを用いた初の系統的ベンチマークであり、大規模バルクデータにおいても非線形手法であるUMAPが従来のPCAやMDSを明確に凌駕することを示した点で、従来の常識と対照的な知見を提供している。
新規性: 本研究は、バルクトランスクリプトームデータにおいて、UMAPが既存の実験デザインでは説明できない隠れた生物学的クラスタリング構造を効率的に特定できることを本研究で初めて示した。特に、COPDデータセットにおける性別交絡因子の同定や、SLE縦断データにおける疾患活動性と連動する患者サブグループ (sG1およびsG2) の層別化は、従来の多段階統計解析を経ることなく、UMAPによる単一の可視化ステップから直感的に導き出された。これは、これまで報告されていない生物学的・臨床的洞察をもたらす新規のアプローチである。
臨床応用: UMAPによる患者サンプルの高精度な層別化は、個別化医療の推進において極めて高い臨床的有用性を持つ。SLE患者の縦断解析において、初期訪問時の臨床スコア (SLEDAI) が同等であっても、その後の病勢が悪化するサブグループ (sG1) と改善するサブグループ (sG2) をトランスクリプトームの軌跡から予測可能であることを示した。この知見は、臨床現場における治療反応性の予測や、最適な治療介入時期の決定を支援するバイオマーカー開発に直接貢献する。
残された課題: 今後の検討課題として、本研究が2次元埋め込み空間での可視化性能に焦点を当てたため、3次元以上の高次元埋め込みにおけるUMAPの性能評価が残されている。また、PCAが持つ変数重み (loadings) のように、各遺伝子の次元削減への寄与度を直接解釈する機能がUMAPには欠けており、解釈可能性の向上が今後の重要な研究方向性である。さらに、本研究の対象がPBMCおよび全血サンプルに限定されていたため、固形がんや他組織のバルクデータに対する外部妥当性の検証が必要である。
方法
データ収集と選定: GEOデータベースから、過去5年間 (2015年1月1日から2020年3月1日) に公開されたヒトPBMC (末梢血単核細胞) または全血由来のバルクトランスクリプトームデータセットを収集した。データセットタイプはアレイまたはハイスループットシーケンスによる発現プロファイリングに限定し、サンプル数は100〜10,000、次元数は1,443〜70,510 (中央値47,878) の範囲とした。初期検索で214件のデータセットが得られたが、各グループのサンプル数が100未満のデータセットを手動で除外し、最終的に71データセットを選定した。
データ前処理と次元削減の実装: 収集したトランスクリプトームカウントデータは、次元削減法の適用前にlog2変換を行った。PCAの場合、サンプルベクトルは平均を除去し、単位分散にスケーリングすることで標準化した。PCA (scikit-learn 0.23.1)、MDS (scikit-learn 0.23.1)、t-SNE (openTSNE 0.6.0)、UMAP (umap-learn 0.5.1) の4つの次元削減法をPythonで実装し、デフォルトパラメータで適用した。具体的には、PCAはn_components = 2、MDSはn_components = 2、dissimilarity = ‘euclidean’、t-SNEはn_components = 2、perplexity = 30、initialization = ‘pca’、metric = ‘euclidean’、UMAPはn_components = 2、n_neighbors = 15、init = ‘spectral’、min_dist = 0.1、metric = ‘euclidean’、random_state = 42 を設定した。実行時間の比較では、t-SNEとUMAPの適用前に、まずPCAを用いてデータの次元を100に削減する前処理戦略を採用した。
クラスタリング構造と定量的評価: 71データセットから生成されたすべての2次元可視化プロットは、データセット情報なしに3名の独立評価者によって手動でクラスタリング構造の有無が評価された (多数決ルール)。手動評価の客観性を検証するため、hdbscan (hierarchical density-based spatial clustering of applications with noise) クラスタリングを2次元埋め込み座標に適用し、シルエットスコアを算出した。21のラベル付きデータセットにおいて、k-means、階層型クラスタリング、スペクトルクラスタリング、GMM (Gaussian mixture model: ガウス混合モデル)、hdbscanの5種類のクラスタリングアルゴリズムを適用し、NMI (Normalized Mutual Information) とARI (Adjusted Rand Index) を用いて真のグループラベルと推論されたグループラベルを比較した。また、埋め込み座標を入力とし、特徴をラベルとするランダムフォレスト分類器を訓練し、ホールドアウトデータ (30%) での予測精度を算出した。局所近傍構造の保存性は、Jaccard index (k=15およびk=30近傍) を用いて評価した。計算時間の測定には、Intel Xeon Gold 5218 (2.30 GHz (gigahertz)、16スレッド) を使用した。
生物学的意義解析: 未解明クラスタリング構造の生物学的意義を事例研究で検討した。COPD (chronic obstructive pulmonary disease: 慢性閉塞性肺疾患) 患者の喫煙習慣影響研究 (GEO: GSE71220) では、UMAPで検出されたクラスタリング構造と性別との関連をχ2検定 (chi-square test) と遺伝子発現ヒートマップで確認した。SLE患者の縦断的転写産物データ (GEO: GSE121239) では、UMAPの異なるmetric (‘euclidean’、‘canberra’、‘cosine’) の影響を比較し、患者サブグループの臨床的意義をSLEDAI (SLE disease activity index: SLE疾患活動性指標) およびGSEA (gene set enrichment analysis: 遺伝子セットエンリッチメント解析、Liberzon et al. CellSyst 2015) で検討した。GSEAにはRパッケージEGSEA (Ensemble Gene Set Enrichment Analysis) を使用した。統計的有意性の検証には one-way ANOVA (一元配置分散分析) および Student t-test を用い、p値0.05未満を有意差ありとした。DGE解析にはRパッケージlimmaを使用した。本研究のベンチマーク検証用として、ヒト肺がん細胞株 A549 および H1299 から得られたバルクRNA-seqデータも参照データとして使用した。