- 著者: Cole Trapnell, Davide Cacchiarelli, Jonna Grimsby, Prapti Pokharel, Shuqiang Li, Michael Morse, Niall J Lennon, Kenneth J Livak, Tarjei S Mikkelsen, John L Rinn
- Corresponding author: John L Rinn (john_rinn@harvard.edu, Harvard University)
- 雑誌: Nature Biotechnology
- 発行年: 2014
- Epub日: 2014-03-23
- Article種別: Original Article (Methods)
- PMID: 24658644
背景
細胞分化・増殖・リプログラミングなどの動的プロセスは複雑な遺伝子制御ネットワークによって駆動されるが、個々の細胞はこれらの変化を大きく異なる速度で実行するため、同時点に回収した細胞集団は複数の分化中間状態の混合物となる。バルク細胞のRNA-Seq時系列解析は、この混合効果によりシンプソンのパラドックスが生じ、個々の細胞で正に相関する因子がバルク平均では逆相関しているように見えるなど、転写カスケードの早期・後期相の区別が困難であった。質量分析フローサイトメトリーを用いた SPADE (Spanning Progression Analysis Density-normalized Events) は複合細胞階層の再構成に有効だが、測定できるタンパク質が最大32種に限られ、全トランスクリプトームのダイナミクスを追跡できない (Satija et al. NatBiotechnol 2015)。マイクロアレイのバルク測定値をバイオロジカルプロセスに沿って時間順に並べる計算幾何学的手法 (Magwene et al. 2003) が先行して開発されていたが、単一細胞の高い発現変動と多系統への分岐に対応できず、重要な制御因子のタイミングを正確に特定できなかった。筋芽細胞分化においては筋原性制御因子 MyoD, MyoG, MRF4 (muscle regulatory factor 4; 筋制御因子4), Myf5 (myogenic factor 5; ミオゲニックファクター5) などの主要転写因子は特定されていたものの、広域的な筋原性制御ネットワークの解明には低い時間分解能が障壁となっており、未発見の制御因子が多数存在すると予想されていた。このような方法論的な不足から、事前知識なしに全トランスクリプトームから分化の進行順を教師なし算出する未解明な課題が残されており、新たな計算ツールが切実に求められていた (Bray et al. NatBiotechnol 2016)。
目的
単一細胞RNA-Seq (scRNA-Seq; single-cell RNA sequencing) データを用いて、マーカー遺伝子の事前知識なしに細胞を生物学的プロセスの進行順に並べる教師なしアルゴリズム「Monocle」を開発すること。ヒト骨格筋芽細胞の分化モデルに適用し、バルク解析では捉えられない転写ダイナミクスと新規制御因子の発見を実証すること。
結果
Monocleアルゴリズムの概要 — 擬似時間による分化軌跡の再構成：
Monocleは5ステップからなる教師なし細胞順位付けアルゴリズムである (Fig 1, 2)。(1) 各細胞の発現プロファイルを高次元ユークリッド空間の点として表現する。(2) ICA (independent component analysis; 独立成分分析) で次元削減し、細胞の発現空間における構造を低次元で可視化する。(3) 細胞上にMST (minimum spanning tree; 最小全域木) を構築する。(4) MSTの最長路 (diameter path) を分化主軌跡として同定する。(5) 分岐を含む擬似時間 (pseudotime; 生物学的進行度の定量指標) を算出する。分岐処理にはPQツリーを用い、ノイズによる小枝を含む実データでも頑健な順位付けを数秒で完了した。発現プロファイルはTopHat2 + CufflinksによりFPKM (fragments per kilobase per million mapped reads; 百万マップリード当たりキロ塩基当たりフラグメント数) として定量し、GAM (generalized additive model; 一般化加法モデル) のTobitファミリーで擬似時間依存的発現変動を検定した。
筋芽細胞分化軌跡の再構成 — 2相構造と間葉系汚染細胞の分離：
ヒト骨格筋芽細胞 HSMM (human skeletal muscle myoblast) を高ミトゲン培地 GM (growth medium) から低血清分化培地 DM (differentiation medium) に切り替え、Fluidigm C1マイクロ流体システムで0, 24, 48, 72時間の各時点にそれぞれn=49-77の細胞を捕捉し (合計n=270細胞)、1細胞あたり約400万リードでHiSeq 2500シーケンスした (Fig 1)。バルクRNA-Seq (n=3生物学的複製) vs. 単一細胞解析の比較でFPKM≥1の遺伝子は各細胞で十分検出された。Monocleは分化軌跡を2つの主要フェーズと1つの分岐に分解した: 第1相は活発な増殖マーカーである CDK1 (cyclin-dependent kinase 1; サイクリン依存性キナーゼ1) 陽性の細胞が主体、第2相はDM移行後の分化細胞 (MyoG陽性) が主体であった。主軌跡から分岐した細胞群は PDGFRA (platelet-derived growth factor receptor alpha; 血小板由来増殖因子受容体α) と Sphk1 (sphingosine kinase 1; スフィンゴシンキナーゼ1) を発現し、間質間葉系細胞 (汚染物) と同定された。これらは免疫精製を行わずにin vitroの細胞間クロストークを維持したまま除外でき、Monocleの多系統分岐処理能力を実証した (Fig 2)。
擬似時間による動的遺伝子の発見 — バルク解析では見えない逐次的制御波：
Tobitモデル + 尤度比検定 (FDR < 5%) により、擬似時間に沿って動的に発現変動する遺伝子が1,061個同定された (Fig 2)。バルクRNA-Seq (n=3複製) vs. 単一細胞擬似時間解析の比較では、免疫蛍光法による Mef2C (myocyte class-C enhancer factor 2; 筋細胞エンハンサー因子2) 陽性細胞増加後に Myh2 (myosin heavy chain 2; ミオシン重鎖2) 陽性細胞が増加するという順序が、擬似時間順序付けにより初めて明確に示された (Fig 2c, d)。また重要なスイッチ的イベントとして、ID1 (inhibitor of DNA binding 1; DNA結合阻害因子1) の不活性化がMyoGの活性化に先行する転写的カスケード (ID1/MyoGスイッチ) が時収集順位付けでは完全に隠蔽されていたのに対し、Monocle擬似時間順位付けでは鮮明に現れた (Fig 2e, f)。バルクRNA-Seq解析では識別できなかった一過性の発現クラスターが単一細胞擬似時間解析で検出され、逐次的な転写再構成の「波」が可視化された。
6つの動態クラスターとcis制御解析 — 175転写因子モチーフが上流候補：
1,061個の動的遺伝子をGAM応答曲線のK-メドイドクラスタリング PAM (Partitioning Around Medoids; 中心点分割法、Rパッケージ) で6クラスターに分類した (Fig 3)。早期下降 / 後期上昇クラスターは筋形成 (細胞周期終了・筋肉特異構造タンパク質) の機能的遺伝子群を高度に濃縮した。各クラスターのDNase I過感受性部位 (HSMM + 筋管、ChromHMM機能分類) にconserved転写因子結合モチーフが175個有意に濃縮された (p<0.001 Bonferroni補正ハイパーゲオメトリー検定、Fig 3)。上昇遺伝子群でMyoD, MyoG, Pbx1 (pre-B-cell homeobox 1; プレB細胞ホメオボックス1), Meis1 (meis homeobox protein 1; ホメオドメイン転写因子), Mef2 (myocyte enhancer factor 2; 筋細胞エンハンサー因子2) ファミリーを含む広範なシス制御ネットワークが明確化された (Fig 3)。バルクRNA-Seqの制御エレメント解析ではこれらの一部しか検出できず、単一細胞擬似時間解析による高解像度の優位性が示された (Supplementary Fig. 11)。
shRNAスクリーニングによる新規分化制御転写因子の実験検証：
上記cis制御解析から選定した11候補転写因子に対して、計44種の短鎖ヘアピンRNA (shRNA; short hairpin RNA) を組み込んだレンチウイルスベクターで増殖期HSMMを感染し、5日間DM培養後にMyh2陽性細胞の頻度・面積・核数をCeligo Sサイトメーターで自動定量した (Fig 4, n=4技術反復)。2種以上の独立shRNAによって分化動態が有意に変化した (FDR < 5%) 因子は MZF1 (myeloid zinc finger 1; 骨髄系亜鉛フィンガー転写因子1), Zic1 (zinc-finger cerebellum 1; 亜鉛フィンガー型転写因子), XBP1 (X-box binding protein 1; Xボックス結合タンパク1), USF1 (upstream stimulatory factor 1; 上流刺激因子1) の4種であり、これらのノックダウンで筋管面積が対照比2.3倍以上増加し、より大型の多核筋管形成を促進した (fold increase ≥2.3)。Cux1 (cut-like homeobox 1; カットホメオボックス1), Arid5b (adenine-thymine-rich interaction domain binding protein), Pou2f1 (pituitary octamer unit-domain factor 1; 下垂体型転写因子), AHR (aryl hydrocarbon receptor; アリール炭化水素受容体) のノックダウンも分化効率を上昇させたが有意差は小さかった。核計数は群間で差がなく (Fig 4a)、増殖速度でなく分化効率の変化であることを確認した。MZF1とUSF1のモチーフはMyoDモチーフを持つエンハンサーで共占有率が高く、競合による分化抑制モデルと一致した (Fig 4b)。
Monocleの頑健性評価 — n=50細胞のサブセットでも高精度を維持：
異なる細胞数でのサブセットシミュレーション実験により、わずかn=50細胞のサブセットでも全データセット (n=270細胞) との順位相関 Spearman r=0.82 を維持し (Supplementary Fig. 7)、全サブセットサイズで動的遺伝子検出の精度 precision=97%を達成した。バルクRNA-Seq (n=3生物学的複製) vs. 単一細胞擬似時間の比較では一過性遺伝子クラスターを追加同定し、バルクでは検出できなかった制御カスケードが高精度で再現された。細胞数増加とともにリコールが改善し、スループットと精度のトレードオフが定量化された。
考察/結論
① 先行研究との違い: バルクRNA-Seq時系列は細胞集団の平均的な転写変化しか捉えられず、シンプソンのパラドックスにより重要な制御因子のタイミングや一過性クラスターを見逃すという限界があった。これと異なり、Monocleは各細胞を独立した「時間点」として扱い、細胞変動そのものを解析の資源とすることで、バルク解析では隠蔽されていたID1/MyoG転換などの逐次的スイッチ状の制御事象を初めて可視化した。またcis制御解析でもバルクと比べて大幅に多くの転写因子モチーフ (175種) が同定され、高い解像度の転写情報が上流制御因子発見に直結することを示した。
② 新規性: 単一細胞発現プロファイルをマーカー遺伝子の事前知識なしに生物学的プロセスの進行度に沿って新規な擬似時間軸上に並べる、初の頑健な教師なしアルゴリズム「Monocle」を提案した。特に多系統分岐の取り扱い (PQツリー)、Tobitモデルに基づく擬似時間依存的差次発現解析、GAM応答曲線を用いたK-メドイドクラスタリングは先行研究にはなかった新規な方法論的貢献である。
③ 臨床応用: 分化・増殖・腫瘍転換など多様な動的プロセスへの適用可能性があり、単一細胞ゲノミクスの基盤技術として広く活用される。希少な汚染細胞サブタイプを分離同定できる点は、がん幹細胞や薬剤耐性集団の解析においても臨床的意義を持つ。cis制御解析との統合により、細胞療法や再生医療に向けた新規分化制御因子候補の同定を加速する橋渡し研究 (translational) ツールとして位置付けられる。
④ 残された課題: 本研究はin vitro培養系のHSMMと骨格筋系列に限定されており、他の細胞種・臓器・生物種での普遍性は今後の検討課題である。Monocle v1ではICAの2次元削減が前提であり、より複雑な多分岐・再収束型分化には対応が困難で、後続研究 (Monocle 2/3) でdiffusion map・UMAP対応が必要となった。また単一細胞RNA-Seqの高ドロップアウト率や技術的ノイズへの対応、大規模 (10,000細胞以上) データへのスケーラビリティも今後の研究に委ねられる (Aibar et al. NatMethods 2017)。
方法
ヒト骨格筋芽細胞 (HSMM、Lonza CC-2580、大腿四頭筋生検、17歳女性、BMI 19、5継代以内) を10% FBS SkGM-2で拡大培養後、MEM-α (minimum essential medium alpha; 最小必須培地α改変型) + 2%馬血清へ切り替えて分化誘導した。Fluidigm C1マイクロ流体システムで0, 24, 48, 72時間の各時点にn=49-77細胞を捕捉し (1チップ/時点、合計n=270細胞)、SMARTer (Switching Mechanism At RNA Template) ULR (ultra-low RNA) キット (Clontech) + Nextera XTでcDNAライブラリを作成し、Illumina HiSeq 2500で100 bp paired-end シーケンス (約4×10^6リード/細胞)。バルクRNA-Seq: n=3生物学的複製。マッピング: TopHat 2.0.9 + Bowtie 2.0.6、ヒトゲノムhg19/GENCODE v17。発現定量: Cuffdiff 2.2 (FPKM)。ICA次元削減: fastICAパッケージ (R)。差次発現: Tobitファミリー GAM + VGAM (vector-generalized additive models; ベクター一般化加法モデル) lrtest()、Benjamini-Hochberg補正 (FDR < 5%)。クラスタリング: K-メドイド PAM (Partitioning Around Medoids)、Rパッケージ。shRNAスクリーン: TRC (The RNAi Consortium; ブロード研究所RNAiコンソーシアム) レンチウイルスライブラリ、対照 = 非標的ルシフェラーゼshRNA、n=4技術反復、pooled分散t検定 + Benjamini-Hochberg補正。データ寄託: NCBI Gene Expression Omnibus (GEO), accession GSE52529。