シングルセルRNA-seq解析

PCA(主成分分析)により発現プロファイルの次元削減を行い、クラスタリング結果をt-SNEにより可視化しました。T細胞(total 4,721 cells)は11のクラスタに分類され、CD8⁺T細胞は5クラスタ(C3,C4,C5,C7,C10、total 1,668 cells)に分類されました。

CD8Aを発現する細胞障害性T細胞(C3,C4,C5,C7,C10)は大部分がGZMA (Granzyme A)を発現していました。

各クラスタに特徴的な遺伝子の発現をヒートマップで一覧表示しました。C5では、疲弊した(exhausted)T細胞に特徴的に発現するHAVCR2(TIM-3), TNFRSF9, PDCD1(PD-1)などが発現していました。この様にして、細胞集団の性質を知らずにクラスタリング(unsupervised clustering)した後で、その細胞の性質を解析出来ます。さらに、由来組織を色分け表示(図上部、末梢血：緑、正常組織：青、腫瘍組織：赤)すると、C5のほとんどは腫瘍組織由来でした。

疲弊したT細胞に特徴的な遺伝子群(Zheng et al.)の発現量と、ランダムに選択した遺伝子群の発現量の差分をスコアリングして可視化しました。スコアが高い(赤色で表示)細胞ほど、疲弊したT細胞の特徴を有していると考えられます。C5が最も高いスコアを示しました。

由来組織ごとに各クラスタの構成比率を見ると、C5は末梢血(PB)や正常組織(Normal)にはほとんど見られず、腫瘍組織(Tumor)浸潤T細胞に多く含まれていました。この様に、細胞の分類基準や組成が分かっていない状態で細胞を分類し、解析対象(病態モデルやノックアウトなど)に特徴的なクラスタを見出すことが出来ます。

11クラスタに分類されたT細胞 (total 4,721 cells)のうち、CD4⁺T細胞は6クラスタ(C0, C1, C2, C6, C8, C9、total 3,053 cells)に分類されました。

CD4を発現するT細胞(C0,C1,C2,C6,C8,C9)のうち、C2には制御性T細胞(Treg)で特徴的に発現するFOXP3の発現が見られました。

FOXP3を発現しているC2とそれ以外のCD4⁺T細胞の遺伝子発現を比較し、発現変動比と有意確率をボルカノプロットで図示しました。C2は制御性T細胞に特徴的なFOXP3, TNFRSF18, TNFRSF9などの発現量が有意に高い(グラフ右上)ことが示されました。

由来組織ごとにCD4⁺ T細胞の構成比率を見ると、末梢血由来に多いC1とC8は組織中には少なく、また、腫瘍組織由来では正常組織由来と比べて、免疫応答を抑制する制御性T細胞(Treg, C2)の割合が顕著に増加していました。

光刺激0時間後および4時間後の視覚野由来細胞のSingle cell RNA-seqデータをbatch effect correction (実験間誤差補正)後に統合し、クラスタリング結果をUMAPによりプロットしました。細胞タイプ名に続く( )内の数字はクラスタ番号を示します。

各クラスタに特徴的な遺伝子の発現量をヒートマップを用いて可視化しました。各クラスタの細胞数を上限100としてダウンサンプリングした結果を示します。上部のカラーバーは、細胞タイプを示します。

神経細胞と推定されたクラスタについて、各遺伝子の発現レベルをバイオリンプロットを用いて可視化しました。Excitatory neuronにはSlc17a7、Inter neuronにはGad1の特異的な発現が見られました。Inter neuronと推定されたクラスタでは、クラスタ毎に異なる神経ペプチドなど(Sst, Vip, Pvalb)の発現が見られました。

クラスタリング結果を、光刺激0時間後および4時間後の細胞に分けて表示しました。各クラスタの細胞構成比や分布に大きな変化は見られませんでした。図中の数字はクラスタ番号を示します。

細胞タイプごとの平均遺伝子発現量を、光刺激0時間後と4時間後で比較し、Volcano plotにより各遺伝子の発現変動比と統計的有意確率を表示しました。細胞タイプごとに固有の発現変動遺伝子が検出されました。

培養神経細胞を用いた実験などにより、神経活動により発現が変化することが知られている初期応答転写因子の細胞タイプごとの発現変化（log₂FC）をヒートマップ表示しました。

各クラスタに特徴的な遺伝子の発現プロファイルをもとに、分化経路を擬時間に沿って再構成し(pseudo temporal ordering)、二次元に可視化しました。左図はクラスタごと、右図は擬時間に従って細胞を色分け表示しました。

擬時間変化に伴って発現量が変動した遺伝子を網羅的に抽出し、各遺伝子の発現量をヒートマップ表示しました。ヒートマップ中央を分化(および擬時間)の始点(monocyte)として、左側がmacrophage、右側がdendritic cellに分化する方向に擬時間経過を表示しました。

時系列解析1の遺伝子から、既知のdendritic cellおよびmonocyte / macrophageマーカー遺伝子を抽出し、擬時間経過に伴う発現量変化をプロットしました。縦軸は発現量、横軸は擬時間経過、着色された各点は個々の細胞を示しています。擬時間は、始点を0、それぞれの経路の終点を100として表示しました。実線がdendritic cell側、点線がmacrophage側に分岐した細胞での発現量変化を示しています。

分化経路推定を行った二次元プロット上に、dendritic cellおよびmonocyte / macrophageマーカー遺伝子の発現量変化をプロットしました。各マーカー遺伝子の擬時間経過に伴う発現増加が確認できます。

エナメル芽細胞と前駆細胞を含むマウス切歯基部の上皮細胞をクラスタリングしました。PCAにより発現プロファイルの次元削減を行い、クラスタリング結果をUMAPにより二次元上に可視化しました。
3次元プロットはこちらからご覧頂けます。※プロットはドラックして回転、マウスのホイールで拡大、縮小が可能です。読み込みに時間が掛かる場合があります。

切歯基部の前駆細胞はゆっくりと細胞周期が進行していること、エナメル芽細胞は切歯の近位から遠位に向けて分化が進むことが知られています。S期およびG2/M期に特徴的な遺伝子と、エナメル芽細胞に特徴的な遺伝子の発現分布を示しました。

RNA velocity解析により細胞動態を解析しました。矢印は状態変化の方向と大きさを示します。前駆細胞の細胞周期進行と、エナメル芽細胞の分化に伴う状態変化が特徴的に検出されました。

マウス切歯基部由来の上皮細胞をクラスタリングし、経路推定(軌道解析)を行いました。推定された経路を実線(黒)で表示しました。
3次元プロットはこちらからご覧頂けます。※プロットはドラックして回転、マウスのホイールで拡大、縮小が可能です。読み込みに時間が掛かる場合があります。

クラスタ4はエナメル芽細胞の前駆細胞、クラスタ22の先端部が分化が進んだエナメル芽細胞と推定されました。クラスタ4から22に至る、4、9、22の分化経路を解析対象としました。クラスタ4を始点として擬時間を色分け表示しました(右図)。

各クラスタに特徴的な遺伝子の発現量をヒートマップで可視化しました。各遺伝子は、発現プロファイルの類似度により、モジュール(遺伝子群)化しました。縦軸は遺伝子群、横軸はクラスタに対応しています。分化経路のクラスタ(4、9、22)を赤枠、最終分化したエナメル芽細胞で発現が増加するモジュールを青枠で示しました。

分化の最終段階で発現が上昇した遺伝子モジュールには、分化したエナメル芽細胞で発現することが知られているAmelogenin (Amelx)、Enamelin (Enam)などが含まれていました。グラフの各点は細胞を示しており、解析対象クラスタに含まれる細胞のみを表しています。

エナメル芽細胞のマーカーであるAmelogenin と、今回の解析で抽出したClusterin (Clu)遺伝子の発現を二次元プロット上で可視化しました。

MEF細胞にOct4, Sox2, Klf4を導入して0, 3, 6, 8日後、および、化学的誘導処理を行って0, 3, 6, 9日後に取得されたデータを統合してクラスタリングした結果をUMAPにより二次元上にプロットしました。各検体の細胞数を3,000細胞に統一して解析を行いました。

多能性幹細胞のマーカー遺伝子として知られるNanogとZfp42の発現分布を図示しました。一部の共通する細胞集団でNanogとZfp42の発現が見られました。

各遺伝子の発現レベルをバイオリンプロットで図示しました。縦軸はそれぞれの遺伝子発現量、横軸は各クラスタを示します。クラスタ13にはNanogおよびZfp42の特徴的な発現が見られ、多能性幹細胞様の細胞と推定されました。

パネル右上に示したクラスタリング結果の赤い円で囲んだクラスタを選択し、diffusion mapにより3次元上にプロットしました。 各クラスタ番号と細胞をプロットした色はパネル右上と対応しています。 3次元プロット上の細胞について、発現プロファイルの類似性からcell lineage（細胞系譜）を推定し、多能性幹細胞を含む細胞系譜を黒い線で表示しました。

3次元上にプロットした細胞分布を検体別に表示しました。OSK群では時間経過により、まず多能性幹細胞の前駆細胞（青丸内）が誘導され、さらに多能性幹細胞（赤丸内）が誘導されたのに対し、CIP群ではそれらの細胞の誘導は確認されませんでした。

多能性幹細胞に至る系譜として推定された系譜上の各クラスタ（3＞5＞8＞2＞9＞0＞13）の細胞構成比を検体別に表示しました。縦軸は各検体の全細胞数を1とした時の各細胞タイプの割合を示し、横軸はクラスタ番号を示します。OSK群では誘導後の時間経過に伴ってクラスタ0、13に含まれる細胞が増加するのに対し、CIP群では増加が全く見られませんでした。この実験条件下では化学的な刺激による多能性幹細胞の誘導は起こらなかった可能性が考えられました。

細胞を二次元上にプロットし、多能性誘導経路の推定を行いました。多能性幹細胞に向かう経路を赤色の矢印で示しました。

推定した経路に従って擬時間（pseudotime）を推定し、擬時間経過による遺伝子発現変化を群別にプロットしました。縦軸は遺伝子発現量、横軸は擬時間経過を示し、図中の各点は個々の細胞を示します。OSK群では擬時間経過に従って、中胚葉に特徴的に発現するMest (Mesoderm Specific Transcript)の発現が減少し、多能性幹細胞のマーカー遺伝子として知られるNanogおよびZfp42の発現上昇が見られました。一方、CIP群では擬時間が経過してもMestの発現が維持されており、 NanogおよびZfp42の発現上昇は見られませんでした。この条件下では、化学的な誘導による多能性幹細胞の誘導は見られず、中胚葉様の性質を維持していたと考えられます。