R、Seurat、CellChatを使用してマウスの皮膚創傷治癒の単一細胞トランスクリプトミクスデータセットを解析

Summary

ここでは、Rを使用してマウスの皮膚創傷治癒の単一細胞経時トランスクリプトミクスデータセットを分析するための段階的な視覚的なワークフローを紹介します。このプロトコルには、Seuratを使用したデータセットのダウンロード、品質管理、視覚化、細胞型注釈、およびCellChatを使用した細胞間相互作用分析のための標準パイプラインが含まれています。

Abstract

創傷治癒のプロセスは、空間と時間を超えた異なる細胞型間の複雑な相互作用によって調節されます。単一細胞トランスクリプトミクス法は、複雑な環境内の個々の細胞のプロファイリングを通じて、細胞の不均一性、細胞コミュニケーションネットワーク、創傷治癒プロセスに関与する細胞間相互作用の調査を可能にします。しかし、多くの単一細胞解析ツールはコンピューターコーディング環境内で実行されており、創傷治癒科学者によるより広範な使用は、バイオインフォマティクスの専門知識が明らかに欠如しているため妨げられています。したがって、RStudioと呼ばれるグラフィカルコーディング環境を使用して、時間的マウス切除皮膚創傷治癒データセットの基本的な単一細胞分析を実行する方法を示すステップバイステップのワークフローが提示されます。この視覚的でガイド付きのプロトコルにより、バイオインフォマティクスのバックグラウンドを持たない科学者は、以前に公開された創傷治癒データセットをダウンロードし、重要な品質管理ステップを実行し、Seuratを使用したデータセットの視覚化と細胞型注釈を含む標準的な単一細胞分析ワークフローを実行し、細胞サブタイプ分析を実行し、モジュールスコアリング分析を実行し、CellChatを使用して細胞間相互作用分析を実行し、Seuratを使用して複数のデータセットの統合分析を実行できます。プロトコルの各ステップについて説明が提供され、すべてのコード行からグラフィカルな結果が表示され、ワークフローを通じてユーザーを安全にガイドします。この単一細胞解析パイプラインへの視覚的な導入の目的は、より多くの創傷治癒科学者が自分の研究室でバイオインフォマティクスツールを直接使用できるようにして、独自の単一細胞データセットのより深い分析と、以前に公開された単一細胞データセットのより広範な再分析を促進することです。

Introduction

創傷治癒は哺乳類生物学で最も複雑なプロセスの1つであり、炎症性、増殖性、^{および解決1,2}の3つの治癒段階のスペクトルが含まれます。これらの治癒段階は、創傷修復の空間と時間を超えて、数十の細胞型とその数百の分子産物の協調作用を広く分類します³。治癒の経時過程にわたる創傷組織サンプリングに基づく数十年の組織学的および分子的研究^{により、組織}修復の包括的な細胞パターンが解明されました3、特に切除された皮膚創傷治癒の再現性のあるマウスモデルにおいて^4,5,6。バルク組織^7、8、9および細胞^{10、11、12、13、14}のスケールでの創傷のハイスループットトランスクリプトーム分析の出現から始まる、創傷治癒の複雑さをより十分に理解することが可能になったのは、過去20年のことでした.ごく最近では、いくつかの研究で皮膚の創傷を単一細胞レベルで転写プロファイリングし、新しい創傷細胞のサブタイプを特定し、治癒中にそれらが互いにどのように相互作用するかを示しています^{15、16、17、18、19、20}。Huらは、革新的な空間単一細胞RNA配列決定アプローチを使用して、創傷中心から数放射状の距離で治癒の時間経過を通じて皮膚創傷をプロファイリングし、空間と時間を超えた新しい細胞間および分子の「動き」を明らかにしました²⁰。このような研究は、創傷治癒の複雑さを前例のないほど詳細に解明しており、細胞と分子の途方もない不均一性の絵を描き始めています。

バイオインフォマティクス解析手法の最近の大きな進歩により、創傷治癒研究の分野で生成されている複雑なマルチオミクスデータセットを生物学的に理解することが可能になりました。Seuratのような単一細胞解析パッケージは、創傷などの複雑な組織における細胞型の分類を含む、データセットの堅牢な分析と統合のためのツールを提供します²¹。単一細胞データの下流の解釈のために、CellChatのようなツールを使用して、細胞が創傷を修復するためにどのように調整するかを説明する可能性のある推定細胞間相互作用プログラムを同定します²²。これらのツールは十分に文書化され、よく引用されていますが、ゲノミクスやトランスクリプトミクスのバイオインフォマティクス分野で最も一般的に使用されている統計的およびグラフィカルなプログラミング言語であるRなどのコンピューターコーディング環境内で実行する必要があります。創傷治癒分野の生物学者や臨床医は、組織修復を研究するために単一細胞アプローチを使用することが増えていますが、SeuratやCellChatなどのツールを自分の研究室で直接使用するために必要なバイオインフォマティクスのトレーニングを受けている人はほとんどいません。これらのバイオインフォマティクスツールの使用におけるこのような障壁は、科学者がバイオインフォマティクスの助けを借りずに独自のデータセットをより深く分析することを妨げるだけでなく、科学者が他のグループによってすでに公開されている豊富な単一細胞データを確実に再分析することも妨げます。

したがって、バイオインフォマティクスのバックグラウンドを持たない科学者が、以前に公開され、公開されている単一細胞創傷治癒データセット²⁰を分析できるようにするための段階的なワークフローをここに提示する。このプロトコルは、RStudioと呼ばれる一般的な無料のグラフィカルRコーディング環境を使用し、この環境をナビゲートして、SeuratとCellChatを使用して複雑な単一セルデータセットの基本的な分析を可能にする所定のコード行を実行する方法を示しています。このプロトコルでは、創傷治癒研究に関連する7つの主要な方法が提示されており、これには、1)コーディング環境のインストール、2)データセットのダウンロードと重要な品質管理ステップ、3)視覚化と細胞型注釈を含む単一細胞分析ワークフロー、4)細胞サブタイプ分析、5)モジュールスコアリング分析、6)細胞間相互作用分析、および7)複数のデータセットの統合分析が含まれます。各メソッド内では、ユーザーがプロトコルと並行して実行するための実際のコードが提供され、すべてのコード行からの実際のグラフィカルな結果が表示され、ワークフローをガイドします。RStudioと基本的な単一細胞解析ワークフローのガイド付き視覚的な導入の主な目的は、より多くの創傷治癒科学者がこれらの強力なツールを直接使用できるようにして、研究分野のより迅速な進歩を可能にすることです。

Protocol

手記： 7つのバイオインフォマティクス手法を詳述する以下のワークフローでは、プロトコルのすべてのステップには、リストされている順序でユーザー自身のRStudioインターフェイスで直接実行する必要があるそれぞれのコードブロックが付属しています。このプロトコルをできるだけ使いやすいものにするために、R スクリプト ファイル (補足ファイル 1: JoVE_Rscript.R) が含まれており、ユーザーの RStudio セッションに直接読み込むことができるため、コードの各行を簡単に実行できます。これにより、ユーザーはプロトコル ドキュメントからコードを入力またはコピーして貼り付ける必要があり、エラーが発生する可能性があります。すべてのプロトコル命令は、各コメント行の先頭にハッシュタグ「#」記号で示されるコメントの形式で R スクリプト ファイルにも含まれています。

1. R、RStudio、およびシングルセル解析ワークフローに必要なRパッケージのインストール

1. R (バージョン 4.4.1) をコンピュータにダウンロードしてインストールします。コンピュータのオペレーティングシステムに対応するリンクを使用します。
   1. Microsoft Windows を実行しているコンピューターを使用している場合は、次のリンクを使用します https://cran.rstudio.com/bin/windows/base/
   2. MacOSを実行しているコンピューターを使用している場合は、次のリンクを使用します https://cran.rstudio.com/bin/macosx/
2. 最新バージョンの RStudio をコンピュータにインストールします。次のリンクをクリックして、指示に従ってください。
   https://www.rstudio.com/products/rstudio/download/#download
3. Rtools (バージョン 4.4) をインストールすると、R が特定のパッケージをコンパイルできるようになります。次のリンクをクリックして、指示に従ってください。
   1. Windows を使用している場合は、次のリンクを使用してください: https://cran.rstudio.com/bin/windows/Rtools/
   2. MacOS を使用している場合は、次のリンクを使用します。
      https://mac.r-project.org/tools/
4. ローカル作業ディレクトリを設定します。これは、すべてのファイルがロードされ、保存されるコンピューター上のフォルダーです。RStudio のメニューバーで [セッション] を選択し、[ 作業ディレクトリの設定] > [ディレクトリの選択] をクリックして 、目的のフォルダを選択して、作業ディレクトリを設定します。
   1. Windows コンピューターを使用している場合は、次のコマンドを使用して作業ディレクトリを設定します。次のコード行の [Directory] を実際のディレクトリ構造に変更します。R のディレクトリ区切り文字は文字 "/" であることに注意してください
      setwd("C:/[Directory]")
   2. MacOSコンピュータを使用している場合、次のコマンドは作業ディレクトリも設定します。次のコード行の [Directory] を実際のディレクトリ構造に変更します。R のディレクトリ区切り文字は文字 "/" であることに注意してください
      setwd("~/[Directory]")
   3. R セッション中の任意の時点で、次のコード行を使用して作業ディレクトリを確認します。
      getwd()
   4. RStudio では、右側のウィンドウの [ファイル] タブで、そこに含まれるすべてのファイルとフォルダーを含む作業ディレクトリ構造を視覚的に探索します。RStudio のファイル エクスプローラーを作業ディレクトリに移動するには、 歯車 アイコン -> [ 作業ディレクトリに移動] をクリックします。
5. プロトコルに必要な依存関係である R パッケージ リポジトリ CRAN から次のパッケージをインストールします。これらのパッケージをインストールするには、次のコマンドを実行します。
   install.packages("devtools")
install.packages("readxl")
install.packages("openxlsx")
install.packages("tidyverse")
install.packages("scCustomize")
   手記： Rパッケージのインストール中に、さまざまなウィンドウが表示されたり消えたりするのは正常です。パッケージのコンパイルを求めるウィンドウが表示されたら、[ はい] をクリックします。パッケージをインストールする前に R を再起動するように求めるウィンドウが表示されたら、[ NO] をクリックします。
6. キュレーションされた R パッケージ リポジトリ Bioconductor から次のパッケージをインストールします
   (https://bioconductor.org/) は、プロトコルに必要な依存関係です。これらのパッケージをインストールするには、次のコマンドを実行します。
   if (! requireNamespace("BiocManager", quietly = TRUE) )
  install.packages("BiocManager")
BiocManager::install("NMF", update=F)
BiocManager::install("ComplexHeatmap", update=F)
BiocManager::install("BiocNeighbors", update=F)
BiocManager::install("SingleCellExperiment", update=F)
BiocManager::install("circlize", update=F)
BiocManager::install("edgeR", update=F)
BiocManager::install("scDblFinder", update=F)
7. この原稿で説明されているワークフローに必要な次のパッケージをインストールします。
   install.packages("Seurat")
devtools::install_github("jinworks/CellChat")
8. 各パッケージにロードして、インストールが成功したことを確認します。いずれかのパッケージで「パッケージが見つかりません」というエラーが発生した場合は、上記の適切なコードを使用して再インストールしてください。
   library(readxl)
library(openxlsx)
library(tidyverse)
library(scCustomize)
library(edgeR)
library(scDblFinder)
library(Seurat)
library(CellChat)

2. 単一細胞創傷治癒データセットへのロードと品質管理ステップの実行

手記： このバイオインフォマティクスワークフローでは、以前に発表された時空間単一細胞皮膚創傷治癒実験の再分析が実施されます²⁰。データセットファイルは、キュレーションされたNCBI遺伝子発現オムニバス(GEO)リポジトリ(https://www.ncbi.nlm.nih.gov/geo/)に保存されます。

1. アクセッション番号GSE204777を使用して、GEOからデータセットファイルに移動します。次のリンクを使用して、研究の実験計画を確認します。 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE204777
2. GEOリポジトリページ内には、5つのシーケンスレーンから取得されたデータの5つの個別のバッチがあります。GSM6190913 というタイトルの最初のデータセットをクリックします。サンプルへの直接リンクは次のとおりです https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM6190913
3. ページの一番下までスクロールし、 ftp または html リンクを使用して次の 3 つのファイルをダウンロードします。コンピュータのファイルエクスプローラ内で、これら3つのファイルを b1というディレクトリに移動します。 b1 フォルダが、手順1.4で設定した作業ディレクトリ内にあることを確認してください。
   ファイル名:GSM6190913_b1_barcodes.tsv.gz / ファイルサイズ:18.5 Mb
   ファイル名: GSM6190913_b1_features.tsv.gz / ファイルサイズ: 254.1 Kb
   ファイル名:GSM6190913_b1_matrix.mtx.gz / ファイルサイズ:151.2 Mb
4. 手順2.3でダウンロードしたシングルセルシーケンシングファイルのディレクトリ情報を取得します。
   b1_data_dir <- file.path(getwd(), "b1")
5. シングルセルシーケンシングファイルをロードします。gene.column パラメーターは、使用する遺伝子/特徴の命名を指定します。この場合、遺伝子シンボルにはgene.column = 2を使用します(gene.column = 1はアンサンブル遺伝子名の場合)。
   b1 <- Read10X_GEO(data_dir = b1_data_dir, gene.column = 2)
   手記： ほとんどの単一セル データセットには追加の多重化データがないため、このステップを使用して生成された 10x ファイルには複数のレイヤーはありません。現在の多重化データセットについては、手順 2.6 に進みます。データを多重化していないデータセットの場合は、手順 2.7 に進みます。
6. 時空間バーコードを使用して、単一セルデータセットを逆多重化します。
   1. 作業データセットの場合は、遺伝子発現と HTO (マルチプレックス) データを分離します。
      dataset_rna_counts <- b1$GSM6190913_b1_$`Gene Expression`
​dataset_barcodes <- b1$GSM6190913_b1_$`Antibody Capture`
   2. 遺伝子発現データを使用してSeuratオブジェクトを作成し、5個未満の細胞で発現する遺伝子と検出された遺伝子が200個未満の細胞を即座に除外します。
      dataset <- CreateSeuratObject(counts = dataset_rna_counts, min.cells = 5, min.features = 200)
   3. 遺伝子発現データセットをレイヤーとして作成し、両方のアッセイに共通する細胞とバーコードのリストを生成します。
      dataset_rna_counts2 <- LayerData(dataset, data = "RNA")
​dataset_joint.barcodes <- intersect(colnames(dataset_rna_counts2), colnames(dataset_barcodes))
   4. 関節細胞バーコードによるサブセット遺伝子発現とHTO数
      dataset_rna_counts2 <- dataset_rna_counts2[, dataset_joint.barcodes]
​dataset_barcodes2 <- as.matrix(dataset_barcodes[, dataset_joint.barcodes])
   5. HTO に予期されるバーコード名があることを確認します。
      ​rownames(dataset_barcodes2)
   6. バーコード情報を保存する新しいアッセイを作成し、このアッセイを以前に作成したSeuratオブジェクトに追加します。
      dataset_barcode_assay <- CreateAssayObject(counts = dataset_barcodes2)
​dataset[["barcodes"]] <- dataset_barcode_assay
   7. オブジェクトに複数のアッセイが含まれていることを確認します。
      DefaultAssay(dataset)
   8. バーコードデータを正規化し、HTODemux関数を介して逆多重化を実行します。この方法は、次の Seurat ビネットで詳しく説明されています https://satijalab.org/seurat/articles/hashing_vignette
      ​dataset <- HTODemux(dataset, assay = "barcodes", positive.quantile = 0.99)
   9. グローバル分類の結果に基づいてセルをグループ化し、バーコード分類のないセルを削除します。
      Idents(dataset) <- "barcodes_classification.global"
​dataset <- subset(dataset, idents = "Negative", invert = TRUE)
   10. 最大HTOシグナルに基づいてセルをグループ化します。
       Idents(dataset) <- "barcodes_maxID"
   11. 多重化されたバーコードごとに細胞内の検出された遺伝子の分布を視覚化します(補足図1)。
       VlnPlot(dataset, features = "nFeature_RNA", pt.size = 0.1, log = TRUE)
   12. バーコードの名前を実際の巻き時間(巻き取り後の日数)とスペース(2〜8 mm)の割り当て(元の原稿から取得)に変更し、time_spaceという新しいメタデータ変数に割り当てます。
       Idents(dataset) <- "barcodes_maxID"
levels(dataset)
dataset <- RenameIdents(dataset,
"Barcode-1" ="D01_2mm",
"Barcode-2" ="D01_4mm",
"Barcode-3" ="D01_6mm",
"Barcode-4" ="D01_8mm",
"Barcode-5" ="D03_2mm",
"Barcode-6" ="D03_4mm",
"Barcode-7" ="D03_6mm",
"Barcode-8" ="D03_8mm",
"Barcode-9" ="D07_2mm",
"Barcode-10" ="D07_4mm",
"Barcode-11" ="D07_6mm",
"Barcode-12" ="D07_8mm",
"Barcode-13" ="D14_2mm",
"Barcode-14" ="D14_4mm",
"Barcode-15" ="D14_6mm",
"Barcode-16" ="D14_8mm",
"Barcode-17" ="UW")
levels(dataset)
dataset[["time_space"]] <- Idents(dataset)
7. データを多重化しないデータセットの場合:Seuratオブジェクトを作成し、5個未満の細胞で発現する遺伝子と検出された遺伝子が200個未満の細胞をすぐに除外します。
   dataset <- CreateSeuratObject(counts = b1, min.cells = 5, min.features = 200)
8. データセットの遺伝子発現アッセイの解析に切り替えます。
   DefaultAssay(dataset) <- "RNA"
Idents(dataset) <- "data"
9. 重要な品質管理ステップとして、各細胞内のミトコンドリア遺伝子の割合を計算し、それをメタデータ変数として割り当てます。この方法は、次の Seurat ビネットで詳しく説明されています https://satijalab.org/seurat/articles/pbmc3k_tutorial#qc-and-selecting-cells-for-further-analysis
   dataset[["percent.mt"]] <- PercentageFeatureSet(dataset, pattern = "^mt-")
10. 検出された遺伝子の分布、RNAの数、およびすべての細胞におけるミトコンドリアの割合を視覚化します(補足図2)。
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
plot1 + plot2
11. ミトコンドリア含有量が多い細胞が多数あり、これはRNA数の低下と相関しています。これらは死んだ細胞または死にかけている細胞です。これらの低品質のセルは、公正なカットオフを使用してデータセットから削除します。この場合、25%を超えるミトコンドリア遺伝子を持つ細胞が除去される、以前に発表された研究²⁰で説明されている元のデータセットの値を使用します。
    dataset <- subset(dataset, subset = nFeature_RNA > 200 & percent.mt < 25)
12. 低品質の細胞を除去した後、検出された遺伝子の分布、RNAの数、およびすべての細胞におけるミトコンドリアの割合を視覚化します(補足図3)。
    plot1 <- FeatureScatter(dataset, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(dataset, feature1 = "nFeature_RNA", feature2 = " percent.mt ")
plot1 + plot2
13. 追加の重要な品質管理ステップとして、データセット内のダブレットの可能性を検出します。これらは液滴配列決定中に結合された細胞であるため、単一細胞レベルではない遺伝子発現が得られます。この問題に対処するために、いくつかのツールが開発されています。各ツールは単一セル データについて特定の仮定を立てているため、ユーザーは 1 つのツールを使用する前に、関連するすべてのドキュメントを読むことが重要です。このワークフローでは、scDblFinder²³ というメソッドを実装します。この方法では、固定の予想ダブレットレートを課すアルゴリズムを使用することに注意してください。次のコマンドを使用して、scDblFinderパイプラインを実行します。
    DefaultAssay(dataset) <- "RNA"
sce <- scDblFinder(GetAssayData(dataset, assay="RNA", slot="counts"), samples=Idents(dataset))
    手記： このような多重化された単一セルデータセットは、複数のバーコードを発現するセルを削除することにより、ダブレットについてスクリーニングすることもできます。ほとんどの単一セル データセットにはこの固有の機能がないため、このワークフローではこの方法が実証されませんでした。代わりに、ダブレットをスクリーニングするためのより一般化可能なパイプラインが、scDblFinderメソッドを使用して示されています。
14. ダブレットスコアを新しいメタデータ変数に割り当てます。
    dataset$scDblFinder.score <- sce$scDblFinder.score
15. すべてのセルにおけるダブレットスコアの分布を視覚化します(補足図4)。
    VlnPlot(dataset, features = "scDblFinder.score", raster=FALSE, pt.size=0.5)
16. 0.25 ダブル スコアのしきい値を超えるセルを削除します。このしきい値は、上記で生成されたバイオリンプロットに基づいて選択され、データセット内のほとんどのセルが非常に高いまたは非常に低いダブルスコアに割り当てられることを示しており、0.25は、このデータセットの妥当なカットオフであり、ありそうもないダブレットを多く削除することなく、ダブレットの可能性が高い大部分を削除します。
    dataset <- subset(dataset, scDblFinder.score < 0.25)
17. データセット Seurat オブジェクトを RDS ファイルとして作業ディレクトリに保存します。
    saveRDS(dataset, "dataset_post_Method2.rds")

3. Seuratを使用した単一細胞創傷治癒データセットの分析

注: (オプションの手順)ここでワークフローを開始する場合は、保存したRDSファイルをSeuratオブジェクトとしてロードします。

dataset <- readRDS("dataset_post_Method2.rds")

1. 単一セルデータセットの正規化、スケーリング、主成分分析(PCA)のための標準的なSeuratワークフローを実行します。この標準ワークフローは、次の Seurat ビネットで説明されています。
   Seurat - ガイド付きクラスタリングチュートリアル: https://satijalab.org/seurat/articles/pbmc3k_tutorial
   Seurat コマンド リスト: https://satijalab.org/seurat/articles/essential_commands
   DefaultAssay(dataset) <- "RNA"
   dataset <- NormalizeData(object = dataset)
dataset <- FindVariableFeatures(object = dataset)
dataset <- ScaleData(object = dataset)
dataset <- RunPCA(object = dataset)
2. 最初の50個のPCAディメンションに関するデータセットの変動量を視覚化します(補足図5)。
   ElbowPlot(dataset, reduction = "pca", ndims = 50)
   主要な変動の多くは、最初の13次元内で発生します。
3. 1〜13の設定PCA次元範囲と0.1の設定解像度を使用して、データセットのセルクラスタリングを実行します。
   dataset <- FindNeighbors(dataset, verbose = TRUE, dims = 1:13)
dataset <- FindClusters(dataset, verbose = TRUE, resolution = 0.1)
   手記： このワークフローは、細胞タイプの大規模な違いに焦点を当てています。したがって、PCAディメンション範囲にはかなり保守的なパラメータが使用されており、最初の13個のディメンションはデータセット内の変動の大部分を表すことが示されています。細胞をより小さく希少なサブタイプに区別するために、これらのまれな細胞サブタイプはデータセットの変動レベルが低い可能性が高いため、ユーザーは下流の分析により多くのディメンションを使用できます。resolution パラメーターの範囲は 0 から 1 で、データセットに課されるカテゴリ分離の大きさを決定します。このパラメータの設定は、ユーザーの研究課題によって異なります。細胞を多数の小規模で希少なサブタイプにクラスタリングするには、下流の分析に高解像度を使用します。このワークフローは、主要な細胞タイプ間のより広範な違いを調査することを目的としているため、0.1というかなり小さな解像度値を使用し、細胞をより少ない大きなグループにクラスター化することが期待されます。手順3.3で説明した設定を使用すると、8つの一意のセルクラスターが識別されます。これらのクラスターは、"seurat_clusters" と呼ばれるメタデータ変数に自動的に割り当てられます。
4. 最初の 13 個の PCA 次元を使用して、UMAP 次元削減と近傍検索解析を実行します。シード番号 123 を追加して、結果のデータ投影の再現性を確保します。
   dataset <- RunUMAP(dataset, verbose = TRUE, dims = 1:13, seed.use = 123)
   手記： UMAP アルゴリズムは確率的であり、次元削減にランダム性を導入します (https://cran.r-project.org/web/packages/umap/vignettes/umap.html の「安定性と再現性」を参照)。一貫したシードを使用すると、アルゴリズムに「最小限の再現性」が得られますが、結果のプロットは、代表的な図や下流の結果に示されているものとはわずかに異なる場合があります。テストの結果、Windows を実行しているコンピューターと MacOS を実行しているコンピューターでは、おそらくこれらのオペレーティング システムでのランダム性の実装が異なるため、結果が特に異なることがわかりました。
5. UMAPプロットで細胞のクラスタリングを視覚化します(図1)。
   DimPlot(dataset, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
   手記： UMAP アルゴリズムのランダム性により、Windows および MacOS を実行しているコンピューターで上記と同じコードを使用して生成された代替 UMAP プロットを示す図に示すように、わずかに異なるプロットが生成される場合があります。クラスターの形状にわずかな違いがあることに注意してください。したがって、ユーザーは、生成されたすべてのデータとプロットを保存してタイムスタンプを付け、クラスターに関するすべての下流分析は、細胞型アノテーションについて以下で説明するように、生物学的理解を念頭に置いて慎重に行うことが不可欠です。
6. 創傷治癒中に細胞がいつどこから来たかを示す元の実験ラベルが含まれているため、UMAPプロットで細胞の創傷時間/空間アノテーションを視覚化します(図2)。
   DimPlot(dataset, group.by = "time_space", raster = FALSE, label = FALSE)
7. 巻線時間/空間注釈を使用したセルクラスターの関連付けのテーブルを生成します。
   table(dataset$time_space, dataset$seurat_clusters)
8. データセット内の主要な細胞型の ID を決定します。これを行うには、すべてのクラスター間で差次発現遺伝子(DEG)を計算します。クラスターの DEG リストを取得し、それらを変数に割り当て、出力を区切り文字で区切られたテキスト ファイルとして作業ディレクトリ内に保存します。
   手記： このステップは CPU を集中的に使用し、ユーザーのハードウェアによっては時間がかかる場合があります。
   Idents(dataset) <- "seurat_clusters"
Cell_markers <- FindAllMarkers(dataset, max.cells.per.ident = 500, only.pos = TRUE, min.pct = 0.10, logfc.threshold = 0.25)
write.csv(Cell_markers, file = file.path(getwd(), "dataset_cluster_markers.txt"))
9. テキストファイルの内容をコピーし、テキストインポートウィザードを使用してカンマ区切り文字と遺伝子名列のIDをテキストとして指定することで、付属のdataset_cluster_markers.txtファイルをスプレッドシート(Excelなど)にダウンロードして開きます。遺伝子名が「テキスト」であることを示すことが重要ですが、そうしないと、Excelは特定の遺伝子名を自動的に日付に変換します(例:Sept7が9月7日に変更されます)。
10. スプレッドシートで、次の推奨パラメータに従って結果をフィルタリングします。
    1. 対数₂ 倍の変化 (log2FC) の減少に従ってすべての行を並べ替えるために、avg_log2FC列を最大から最小にランク付けします。
    2. クラスター列を最小から最大にランク付けして、増加する Seurat クラスター数に従ってすべての行を配置します。
    3. avg_log2FC列で2.5以上の数値をフィルタリングして、示されたクラスターと他のクラスターで最も差次的に発現する遺伝子(DEG)のみを表示します。
    4. pct.1 列で 0.4 以上の数値をフィルタリングします。この列は、示されたクラスター内の細胞が示された遺伝子を発現する割合(クラスター%)を参照し、しきい値を0.4に設定すると、示されたクラスター内の細胞の少なくとも40%で発現する遺伝子のみが表示されることを意味します。
    5. pct.2 列をフィルタリングして、0.2 以下の数値を探します。この列は、示されたクラスターに示された遺伝子を発現するNOT細胞の割合(Non-cluster %)を指し、しきい値を0.2に設定すると、示されたクラスターにない細胞の最大20%で発現する遺伝子のみが表示されることを意味します。
    6. p_val_adj列をフィルタリングして、0.01 以下の数値を探します。この列は、調整されたP値または偽発見率(FDR)を参照し、示されたDEGの統計的強度を示し、しきい値を0.01に設定すると、FDRが0.01<遺伝子のみが表示されることを意味します。
       注:補足表1 (JoVE_DEGs_cellMarkers.xlsx)には、プロトコルのステップ3.10で使用されたランク付けされた差次的に発現する遺伝子の完全な出力が含まれています。 補足表2は 、この分析における各クラスターの上位5つの遺伝子を示しており、太字の遺伝子はその後の視覚化に使用されています。
11. クラスターの偏りのないセル型アノテーションには、EnrichR Webベースのエンリッチメント分析ツールを使用します。
    リンクを使用する: https://maayanlab.cloud/Enrichr/
12. 各クラスターの DEG のリストを個別の EnrichR ウィンドウにコピーし、[分析] をクリックします。EnrichRツールは、何百ものキュレーションされたデータベースを通じて遺伝子リストを実行し、各カテゴリの濃縮されたすべての用語をランク付けします。
13. セルタイプアノテーションの目的で、上の[セルタイプ]タブをクリックし、左側にある3つのセルマーカーキュレーションデータベースの上位5つのエンリッチメントに焦点を当てます(図3)。
    セルマーカー 2024 (http://bio-bigdata.hrbmu.edu.cn/CellMarker/)
    タブラ・サピエンス (https://tabula-sapiens-portal.ds.czbiohub.org/)
    PanglaoDB拡張(https://panglaodb.se/)
14. これらのデータベース内の DEG のエンリッチメントに基づいて、8 つのクラスターの考えられる同一性を確認します。線維芽細胞として濃縮する2つのクラスター(2、6)があることに注意してください。したがって、これらのクラスターを単一セルタイプの注釈に結合します。セルタイプの ID をラベルとして cell_types という新しいメタデータ変数に割り当てます。
    Idents(dataset) <- "seurat_clusters"
dataset[["cell_types"]] <- Idents(dataset)
Idents(dataset) <- "cell_types"
dataset <- RenameIdents(dataset,
"0" = "Macrophage",
"1" = "Neutrophil",
"2" = "Fibroblast",
"3" = "Epithelial cell",
"4" = "Endothelial cell",
"5" = "T cell",
"6" = "Fibroblast",
"7" = "Smooth muscle cell"
)
levels(dataset)
dataset[["cell_types"]] <- Idents(dataset)
15. 名前を変更した細胞クラスターを UMAP プロット上の注釈として視覚化します (図 4)。
    DimPlot(dataset, group.by = "cell_types", raster = FALSE, label = FALSE)
16. 表1の上位(太字)クラスターマーカー遺伝子の局在を一連のUMAPプロットで視覚化します(図5)。
    DefaultAssay(dataset) <- "RNA"
FeaturePlot(dataset, features = c("Arg1", "Retnlg", "Fgf7", "Dsp", "Tie1", "Cd3g", "Cdh4", "Rgs5"), raster=FALSE, ncol = 4)
17. 元のクラスター番号でグループ化されたドットプロット上の上部クラスターマーカーDEGを視覚化します(補足図6)。
    DotPlot(dataset, group.by = "seurat_clusters", features = c("Arg1","Retnlg","Fgf7","Dsp","Tie1", "Cd3g","Cdh4","Rgs5")) + RotatedAxis() + scale_colour_gradient2(low = "blue", mid = "grey", high = "red")
18. 注釈付きの細胞タイプでグループ化されたドットプロット上の上部クラスターマーカーDEGを視覚化します(図6)。
    DotPlot(dataset, group.by = "cell_types", features = c("Arg1","Retnlg","Fgf7","Dsp","Tie1", "Cd3g","Cdh4","Rgs5")) + RotatedAxis() + scale_colour_gradient2(low = "blue", mid = "grey", high = "red")
19. 時系列解析を実行するには、まずデータセットを簡略化して空間コンポーネントを削除します。経時分析では、「DPW」と呼ばれる新しいメタデータ変数を使用して、創傷の時間/空間アノテーションを創傷後全体の日数 (DPW) にグループ化します。
    Idents(dataset) <- "time_space"
dataset[["DPW"]] <- Idents(dataset)
Idents(dataset) <- "DPW"
new.cluster.ids <- c("D1", "D1", "D1", "D1", "D3", "D3", "D3", "D3", "D7", "D7", "D7", "D7", "D14", "D14", "D14", "D14", "UW")
names(new.cluster.ids) <- levels(dataset)
dataset <- RenameIdents(dataset, new.cluster.ids)
dataset[["DPW"]] <- Idents(dataset)
Idents(dataset) <- "DPW"
20. UMAPプロットで新しい創傷の時間経過グループを視覚化します(補足図7)。
    DimPlot(dataset, group.by = "DPW", raster = FALSE, label = FALSE)
21. 各DPWで発生する各タイプのセルの数を示すテーブルを生成します。
    ​table(dataset$DPW, dataset$cell_types)
    1. オプションの手順: 巻線の時間経過グループの DEG リストも取得するには、それらを変数に割り当て、出力を区切り文字で区切られたテキスト ファイルとして保存します。
       Idents(dataset) <- "DPW"
Cell_DPW_markers <- FindAllMarkers(dataset, max.cells.per.ident = 500, only.pos = TRUE, min.pct = 0.10, logfc.threshold = 0.25)
write.csv(Cell_DPW_markers, file = file.path(getwd(), "dataset_DPW_markers.txt"))
22. 細胞数をカテゴリごとの比率に変換して、治癒の経過全体にわたる細胞型組成の相対的な変化をよりよく理解します。各細胞タイプにおけるDPWの割合を視覚化します(補足図8):
    pt1 <- table(dataset$DPW, dataset$cell_types)
pt1 <- as.data.frame(pt1)
pt1$Var1 <- as.character(pt1$Var1)
ggplot(pt1, aes(x = Var2, y = Freq, fill = Var1)) +
theme_bw(base_size = 15) +
geom_col(position = "fill", width = 0.5) +
xlab("Sample") + RotatedAxis() +
ylab("Proportion") +
theme(legend.title = element_blank())
23. 各DPWの細胞型の割合を視覚化します( 図7)。
    pt2 <- table(dataset$cell_types, dataset$DPW)
pt2 <- as.data.frame(pt2)
pt2$Var1 <- as.character(pt2$Var1)
ggplot(pt2, aes(x = Var2, y = Freq, fill = Var1)) +
theme_bw(base_size = 15) +
geom_col(position = "fill", width = 0.5) +
xlab("Sample") +
ylab("Proportion") +
theme(legend.title = element_blank())
24. データセット Seurat オブジェクトを RDS ファイルとして作業ディレクトリに保存します。
    saveRDS(dataset, "dataset_post_Method3.rds")

4. Seuratを使用した細胞サブタイプの分析

手記： 単一細胞解析の力により、上記で解析した主要な細胞型内の希少サブタイプの発見と解析が可能になります。この例は、最初に2つのスーラクラスターにクラスター化されてから1つのカテゴリに結合された線維芽細胞に焦点を当てています。プロトコルのこの部分は、他のすべての細胞型を除外して、特に線維芽細胞に焦点を当て、創傷治癒中の線維芽細胞の同一性と時間的特性を調査します。オプションの手順として、必要に応じて、保存したRDSファイルをSeuratオブジェクトとしてロードします。

dataset <- readRDS("dataset_post_Method3.rds")

1. 線維芽細胞の同一性に従って元のデータセットをサブセット化します。
   アイデント(データセット) <- "cell_types"
   dataset_fibroblast <- subset(dataset, idents = "Fibroblast")
2. この小さなデータセットでPCAを実行し、PCA次元に対するデータセットの変動量を視覚化します(補足図9)。
   dataset_fibroblast <- RunPCA(dataset_fibroblast, verbose = TRUE)
ElbowPlot(dataset_fibroblast, reduction = "pca", ndims = 50)
   手記： 主要な変動の多くは、最初の9次元内で発生します。
3. 1〜9の設定PCA次元範囲と0.1の設定解像度を使用して、データセットのセルクラスタリングを実行します。
   dataset_fibroblast <- FindNeighbors(dataset_fibroblast, verbose = TRUE, dims = 1:9)
dataset_fibroblast <- FindClusters(dataset_fibroblast, verbose = TRUE, resolution = 0.1)
   手記： これらの設定を使用して、アルゴリズムは3つの固有の線維芽細胞クラスターを識別します。これらのクラスターは、seurat_clusters というメタデータ変数に自動的に割り当てられます。
4. 最初の 9 つの PCA 次元を使用して、UMAP 次元削減と近傍検索解析を実行します。シード番号 123 を追加して、結果のデータ投影の再現性を確保します。
   dataset_fibroblast <- RunUMAP(dataset_fibroblast, verbose = TRUE, dims = 1:9, seed.use = 123)
5. UMAPプロットで細胞のクラスタリングを視覚化します(図8)。
   DimPlot(dataset_fibroblast, group.by = "seurat_clusters", raster = FALSE, label = TRUE)
6. UMAPプロットで細胞の創傷経時アノテーションを視覚化します(補足図10)。
   DimPlot(dataset_fibroblast, group.by = "DPW", raster = FALSE, label = FALSE)
7. 3つの線維芽細胞サブタイプのDEGリストを取得し、作業ディレクトリのテキストファイルに保存します。
   Idents(dataset_fibroblast) <- "seurat_clusters"
fibroblast_markers <- FindAllMarkers(dataset_fibroblast, max.cells.per.ident = 500, only.pos = TRUE, min.pct = 0.10, logfc.threshold = 0.25)
write.csv(fibroblast_markers, file = file.path(getwd(), "fibroblast_cluster_markers.txt"))
8. スプレッドシートで上記と同様の手順(ステップ3.9〜3.10)に従って、DEGを上位の線維芽細胞サブタイプマーカーにフィルタリングします。
9. DEGテキストファイルから3つのクラスターのそれぞれの上位5つの線維芽細胞サブタイプマーカーをコピーすることにより、カスタム遺伝子リストを変数として定義します。
   FB_type_marker <- c("Plac8", "Cthrc1", "Adam12", "Ifi211", "Plat", "Cdh4", "Aldh3a1", "Ppp1r14a", "Oxtr", "Serpina3c", "Sostdc1", "Scube3", "Ptprz1", "Corin", "Alx4")
10. 線維芽細胞のみのデータセットのリスト内の遺伝子を視覚化するには、ドットプロットのfeaturesパラメータで変数を呼び出します(図9)。
    DotPlot(dataset_fibroblast, group.by="seurat_clusters", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
DotPlot(dataset_fibroblast, group.by="DPW", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
11. ドットプロットのfeaturesパラメータで変数を呼び出して、元の単一細胞データセットのリスト内の遺伝子を視覚化します(補足図11)。
    DotPlot(dataset, group.by="cell_types", features = FB_type_marker) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
12. 各DPWカテゴリーにおける線維芽細胞サブタイプの割合を視覚化します(補足図12)。
    pt3 <- table(dataset_fibroblast$seurat_clusters, dataset_fibroblast$DPW)
pt3 <- as.data.frame(pt3)
pt3$Var1 <- as.character(pt3$Var1)
ggplot(pt3, aes(x = Var2, y = Freq, fill = Var1)) +
theme_bw(base_size = 15) +
geom_col(position = "fill", width = 0.5) +
xlab("Sample") +
ylab("Proportion") +
theme(legend.title = element_blank())
13. 各線維芽細胞サブタイプカテゴリにおけるDPW線維芽細胞の割合を視覚化します(補足図13)。
    pt4 <- table(dataset_fibroblast$DPW, dataset_fibroblast$seurat_clusters)
pt4 <- as.data.frame(pt4)
pt4$Var1 <- as.character(pt4$Var1)
ggplot(pt4, aes(x = Var2, y = Freq, fill = Var1)) +
theme_bw(base_size = 15) +
geom_col(position = "fill", width = 0.5) +
xlab("Sample") +
ylab("Proportion") +
theme(legend.title = element_blank())
14. データセット Seurat オブジェクトを RDS ファイルとして作業ディレクトリに保存します。
    saveRDS(dataset_fibroblast, "dataset_fibroblast_post_Method4.rds")

5. モジュールスコアリングによるフォローアップ分析の例

手記： 単一セル データセットを分析するための便利な方法の 1 つは、モジュール スコアリングと呼ばれます。このワークフローでは、以前の知識に従って遺伝子リストを定義し、モジュールスコアを計算して、各細胞内の遺伝子リストの潜在的な濃縮を特定できます。これらのスコアは、細胞アノテーション全体で平均化して、濃縮の潜在的なパターンを明らかにすることができます。

ここでは、治癒連続体全体からのバルクRNA配列決定サンプルを使用して創傷治癒期特異的遺伝子が同定された、以前に発表された研究²の遺伝子リストを使用します。遺伝子リストはタブ区切りのテキストファイル(補足ファイル2:JoVE_PhaseSpecificGenes.txt)に保存され、作業ディレクトリにダウンロードして、3つの主要な治癒段階を特定する遺伝子リストを生成するために使用できるようになりました。

1. TEXTファイルを読み込んで、遺伝子リストを変数に読み込みます。
   PhaseSpecificGenes <- read_delim("JoVE_PhaseSpecificGenes.txt", delim="\t", col_names = T)
2. 列を個々の遺伝子リスト変数に分離し、遺伝子を最初の文字が大文字のマウス名に変更します。
   PS_Inflammatory <- PhaseSpecificGenes[1]
PS_Inflammatory_ms <- lapply(PS_Inflammatory, str_to_sentence)
PS_Proliferative <- PhaseSpecificGenes[2]
PS_Proliferative_ms <- lapply(PS_Proliferative, str_to_sentence)
PS_Resolution <- PhaseSpecificGenes[3]
PS_Resolution_ms <- lapply(PS_Resolution, str_to_sentence)
   注: (オプション) 必要に応じて、保存した RDS ファイルを Seurat オブジェクトとしてロードします。
   dataset <- readRDS("dataset_post_Method3.rds")
3. 遺伝子リストをモジュールとして使用し、治癒の3つの段階に従ってデータセット内の各細胞をスコアリングします。
   dataset <- AddModuleScore(
object = dataset,
features = PS_Inflammatory_ms,
ctrl = 100,
name = 'Inflammatory'
)
dataset <- AddModuleScore(
object = dataset,
features = PS_Proliferative_ms,
ctrl = 100,
name = 'Proliferative'
dataset <- AddModuleScore(
object = dataset,
features = PS_Resolution_ms,
ctrl = 100,
name = 'Resolution'
)
4. DPWや主要な細胞タイプを含む、細胞カテゴリーごとの集計モジュールスコアを視覚化します(図10)。
   DotPlot(dataset, group.by="DPW", features = c("Inflammatory1","Proliferative1", "Resolution1")) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")
DotPlot(dataset, group.by="cell_types", features = c("Inflammatory1","Proliferative1", "Resolution1")) + RotatedAxis() + scale_colour_gradient2(midpoint = 0, low = "blue", mid = "grey", high = "red")

6. CellChatによるフォローアップ分析の例

手記： 単一細胞データセットの分析に有用でよく引用されているもう 1 つの方法は、細胞間相互作用の推論です。このワークフローでは、細胞群間の異なるリガンド-受容体相互作用を分析することにより、細胞間コミュニケーションを推測するパッケージCellChatを使用します²²。最近、CellChat の開発者は、その一般化された使用のための詳細なステップバイステップのプロトコルを公開しました²⁴、これは、ユーザーが次のワークフローを実行し、それをデータセットに適用する際の優れたリソースです。例として、次のワークフローは、創傷後 1 日と 14 日(DPW)における創傷内のすべての主要細胞の相互作用を比較します。すべてのステップは、公式の CellChat 出版物²⁴ とチュートリアルですでに説明されているため、各ステップは詳細に説明されていません。

を用いた細胞間コミュニケーションの推論と解析
セルチャット: https://github.com/jinworks/CellChat/blob/master/tutorial/CellChat-vignette.html

CellChatを使用した複数のデータセットの比較分析:
https://github.com/jinworks/CellChat/blob/master/tutorial/Comparison_analysis_of_multiple_datasets.html

オプションの手順: 必要に応じて、保存した RDS ファイルを Seurat オブジェクトとしてロードします。

dataset <- readRDS("dataset_post_Method3.rds")

1. 元のデータセットを DPW アノテーションに従って 2 つのデータセットにサブセット化します。
   Idents(dataset) <- "DPW"
dataset_D1 <- subset(dataset, ident = "D1")
dataset_D14 <- subset(dataset, ident = "D14")
2. CellChat を実行する注釈を定義します---、この場合は主要なセルタイプを使用します。
   Idents(dataset_D1) <- "cell_types"
Idents(dataset_D14) <- "cell_types"
3. CellChat オブジェクトを作成し、一般的な CellChat ワークフローに従います。各ステップの詳細なリファレンスとして、上記のリンク先のチュートリアルを参照してください。
   cellchat_D1 <- createCellChat(dataset_D1, group.by = "ident", assay = "RNA")
cellchat_D14 <- createCellChat(dataset_D14, group.by = "ident", assay = "RNA")
CellChatDB <- CellChatDB.mouse
CellChatDB.use <- CellChatDB
cellchat_D1@DB <- CellChatDB.use
cellchat_D14@DB <- CellChatDB.use
cellchat_D1 <- subsetData(cellchat_D1)
cellchat_D14 <- subsetData(cellchat_D14)
future::plan("multisession", workers = 4)
cellchat_D1 <- identifyOverExpressedGenes(cellchat_D1, do.fast = F)
cellchat_D14 <- identifyOverExpressedGenes(cellchat_D14, do.fast = F)
cellchat_D1 <- identifyOverExpressedInteractions(cellchat_D1)
cellchat_D14 <- identifyOverExpressedInteractions(cellchat_D14)
cellchat_D1 <- computeCommunProb(cellchat_D1, type = "triMean", population.size = TRUE)
cellchat_D14 <- computeCommunProb(cellchat_D14, type = "triMean", population.size = TRUE)
cellchat_D1 <- filterCommunication(cellchat_D1, min.cells = 10)
cellchat_D14 <- filterCommunication(cellchat_D14, min.cells = 10)
cellchat_D1 <- computeCommunProbPathway(cellchat_D1)
cellchat_D14 <- computeCommunProbPathway(cellchat_D14)
cellchat_D1 <- aggregateNet(cellchat_D1)
cellchat_D14 <- aggregateNet(cellchat_D14)
cellchat_D1 <- netAnalysis_computeCentrality(cellchat_D1, slot.name = "netP")
cellchat_D14 <- netAnalysis_computeCentrality(cellchat_D14, slot.name = "netP")
4. 各創傷治癒時点でのすべての主要な細胞型における入ってくる相互作用と出る相互作用の強さを視覚化します(補足図14)。
   netAnalysis_signalingRole_scatter(cellchat_D1)
netAnalysis_signalingRole_scatter(cellchat_D14)
   線維芽細胞は、D1 DPW と D14 DPW の間の相互作用を劇的に増加させます。
5. 推定されるすべての重要な細胞間コミュニケーション経路のリストを表示します。
   cellchat_D1@netP$pathways
cellchat_D14@netP$pathways
   コラーゲン経路は、D1 DPW と D14 DPW の両方で重要な経路の 1 つです。
6. コラーゲンシグナル伝達経路と線維芽細胞との相互作用に焦点を当てます。
   pathways.show <- c("COLLAGEN")
7. 円図を使用して、細胞型間のコラーゲンシグナル伝達経路の相互作用を視覚化します(補足図15)。
   par(mfrow=c(1,2))
netVisual_aggregate(cellchat_D1, signaling = pathways.show, layout = "circle")
netVisual_aggregate(cellchat_D14, signaling = pathways.show, layout = "circle")
par(mfrow=c(1,1))
8. コード図を使用して、細胞型間のコラーゲンシグナル伝達経路の相互作用を視覚化します(補足図16)。
   par(mfrow=c(1,2))
strwidth <- function(x) {0.5}
netVisual_aggregate(cellchat_D1, signaling = pathways.show, layout = "chord", vertex.label.cex = 0.6)
netVisual_aggregate(cellchat_D14, signaling = pathways.show, layout = "chord", vertex.label.cex = 0.6)
par(mfrow=c(1,1))
9. 線維芽細胞をソース細胞とするCOLLAGENシグナル伝達経路の相互作用を視覚化します(補足図17)。
   手記： cellchatオブジェクトの細胞タイプは、元のSeuratオブジェクトで割り当てられた順序でIDとしてリストされます:1 =マクロファージ、2 =好中球、3 =線維芽細胞、4 =上皮細胞、5 =内皮細胞、6 = T細胞、7 =平滑筋細胞。
   par(mfrow=c(1,2))
strwidth <- function(x) {0.5}
netVisual_aggregate(cellchat_D1, signaling = pathways.show, layout = "chord", vertex.label.cex = 0.6, sources.use = 3)
netVisual_aggregate(cellchat_D14, signaling = pathways.show, layout = "chord", vertex.label.cex = 0.6, sources.use = 3)
​par(mfrow=c(1,1))
10. 線維芽細胞をソース細胞として、COLLAGENシグナル伝達経路における各リガンド-受容体ペアの寄与を視覚化します。
    1. バブルプロットの使用(補足図18):
       gg1 <- netVisual_bubble(cellchat_D1, sources.use = 3, targets.use = NULL, signaling = pathways.show, remove.isolate = FALSE)
gg2 <- netVisual_bubble(cellchat_D14, sources.use = 3, targets.use = NULL, signaling = pathways.show, remove.isolate = FALSE)
​gg1 + gg2
    2. コード図の使用 (補足図 19):
       par(mfrow=c(1,2))
strwidth <- function(x) {0.4}
netVisual_chord_gene(cellchat_D1, sources.use = 3, targets.use = NULL, signaling = pathways.show, lab.cex = 0.6, show.legend= F)
netVisual_chord_gene(cellchat_D14, sources.use = 3, targets.use = NULL, signaling = pathways.show, lab.cex = 0.6, legend.pos.x = 60)
​par(mfrow=c(1,1))
11. COLLAGEN シグナル伝達経路内の Col1a1-Cd44 リガンドと受容体の相互作用に焦点を当てます。
    ​LR.show <- "COL1A1_CD44"
12. コード図を使用して、細胞型間のCol1a1-Cd44リガンド-受容体相互作用を視覚化します(補足図20)。
    strwidth <- function(x) {0.5}
netVisual_individual(cellchat_D1, signaling = pathways.show, pairLR.use = LR.show, layout = "chord")
netVisual_individual(cellchat_D14, signaling = pathways.show, pairLR.use = LR.show, layout = "chord")
13. 結合された CellChat オブジェクトを生成することにより、差分 CellChat 分析を実行します。
    object.list_D14_v_D1 <- list(D1 = cellchat_D1, D14 = cellchat_D14)
cellchat_D14_v_D1 <- mergeCellChat(object.list_D14_v_D1, add.names = names(object.list_D14_v_D1))
14. 創傷治癒時点における細胞間相互作用の総数と相対的な強さを視覚化します(補足図21)。
    gg1 <- compareInteractions(cellchat_D14_v_D1, show.legend = F)
gg2 <- compareInteractions(cellchat_D14_v_D1, show.legend = F, measure = "weight")
gg1 + gg2
15. 円を使用して視覚化し、創傷が1日目から14日目に移行する際の各細胞型間の細胞間相互作用の強さの違いをプロットします(補足図22)。
    netVisual_diffInteraction(cellchat_D14_v_D1, weight.scale = T, measure = "weight")
16. ヒートマップを使用して、創傷が1日目から14日目に移行する際の各細胞型間の細胞間相互作用の強さの違いを視覚化します(補足図23)。
    netVisual_heatmap(cellchat_D14_v_D1, measure = "weight")
17. ランクプロットを使用して、14日目と1日目の線維芽細胞との細胞間相互作用に対する個々の経路の相対的な寄与を視覚化します(補足 図24)。
    rankNet(cellchat_D14_v_D1, mode = "comparison", measure = "weight", sources.use = 3, targets.use = NULL, stacked = T, do.stat = TRUE)
18. バブルプロットを使用して、1日目と比較して、線維芽細胞をソース細胞として14日目のコラーゲンシグナル伝達経路における個々のリガンド-受容体ペアの相対的な寄与を視覚化します(補足図25)。
    gg1 <- netVisual_bubble(cellchat_D14_v_D1, sources.use = 3, targets.use = NULL, signaling = pathways.show, comparison = c(1, 2), max.dataset = 2, title.name = "Increased signaling in D14", angle.x = 45, remove.isolate = F)
gg2 <- netVisual_bubble(cellchat_D14_v_D1, sources.use = 3, targets.use = NULL, signaling = pathways.show, comparison = c(1, 2), max.dataset = 1, title.name = "Decreased signaling in D14", angle.x = 45, remove.isolate = F)
gg1 + gg2
19. Seurat オブジェクトと同様に、CellChat オブジェクトは RDS ファイルとして保存して開くことができます。
    saveRDS(cellchat_D1, "cellchat_D1.rds")
saveRDS(cellchat_D14, "cellchat_D14.rds")
saveRDS(cellchat_D14_v_D1, "cellchat_D14_v_D1.rds")
20. オプションの手順: CellChat オブジェクトは RDS ファイルから開くこともできます。
    cellchat_D1 <- readRDS("cellchat_D1.rds")
cellchat_D14 <- readRDS("cellchat_D14.rds")
cellchat_D14_v_D1 <- readRDS("cellchat_D14_v_D1.rds")

7. 複数の単一セルデータセットを組み合わせた統合解析の例

手記： 単一細胞データセットは、バッチまたはグループで配列決定されるため、多くの場合、複数のファイルに分割されます。このワークフローは、創傷治癒データセット ²⁰の5つのバッチのうち2つを統合する方法を示しています。データセット統合の現在の方法は、次の Seurat ビネットで説明されています。

scRNA-seq統合の概要:
https://satijalab.org/seurat/articles/integration_introduction

Seurat v5 での統合分析:

https://satijalab.org/seurat/articles/seurat5_integration

手記： 単一セル データセットの統合方法には数多くあり、それぞれに独自の長所と短所があります。詳細については、積分方法の包括的なベンチマーク²⁵ を参照してください。ユーザーは、1 つの統合方法に依存する前に、関連するすべてのドキュメントを読むことが重要です。

1. Hu et al. データセット²⁰ の別のバッチに対して、方法 2 のすべての手順を繰り返します。次のプロトコルでは、バッチ #3 が使用されます。補足 R スクリプト ファイルが含まれており、バッチ #3 (補足ファイル 3: JoVE_Rscript_b3 の処理に使用できます。R)。以下のコードでは、データセットのバッチ #3 に "dataset_b3" を使用して、データセットに新しい変数を作成して使用することを忘れないでください---。
   1. オプションの手順: 必要に応じて、作業ディレクトリに保存されている RDS ファイルから 2 つのデータセットを Seurat オブジェクトとして開きます。
      dataset <- readRDS("dataset_post_Method2.rds")
dataset_b3 <- readRDS("dataset_b3_post_Method2.rds")
2. 「バッチ」という名前の新しい変数を各データセットに割り当てて、後続の解析で原点のデータセットにラベルを付けます。
   dataset@meta.data$batch <- "b1"
dataset_b3@meta.data$batch <- "b3"
3. 方法3で説明されているように、2つのデータセットのSeuratマージを実行し、バッチベースのセルIDアノテーションを追加してから、マージされたデータセットに対して標準のSeuratワークフローを実行します。
   dataset_merged <- merge(x = dataset, y = c(dataset_b3), add.cell.ids = c("b1", "b3"), merge.data = TRUE)
DefaultAssay(dataset_merged) <- "RNA"
dataset_merged <- NormalizeData(dataset_merged)
dataset_merged <- FindVariableFeatures(dataset_merged)
dataset_merged <- ScaleData(dataset_merged)
dataset_merged <- RunPCA(dataset_merged)
ElbowPlot(dataset_merged, reduction = "pca", ndims = 50)
4. データ統合の前に、結合されたデータセットに対してクラスタリングと UMAP 分析を実行します。
   dataset_merged <- FindNeighbors(dataset_merged, dims = 1:10, reduction = "pca")
dataset_merged <- FindClusters(dataset_merged, resolution = .1, cluster.name = "unintegrated_clusters")
dataset_merged <- RunUMAP(dataset_merged, dims = 1:10, seed.use = 123, reduction = "pca", reduction.name = "umap.unintegrated")
5. クラスター番号とバッチ番号に従ってUMAPプロットを視覚化します(補足図26)。
   DimPlot(dataset_merged, reduction = "umap.unintegrated", group.by = c("seurat_clusters", "batch"))
6. バッチ番号に従って、各クラスターのセル番号の分布を表示します。
   table(dataset_merged$batch, dataset_merged$seurat_clusters)
   手記： UMAP プロットと表から、これら 2 つのデータセットに有意なバッチ効果はないようです。バッチ効果の証拠は、2つのデータセット間のクラスター分布の予期しない不一致として現れ、これはデータセット間に実際の生物学的類似性を上書きする潜在的な技術的違いがあることを意味する可能性があります。
7. RPCA方式を使用してSeuratデータ統合を実行します。この方法やその他のデータ統合方法の詳細については、上記のリンク先の Seurat ビネットをお読みください。
   dataset_merged <- IntegrateLayers(
object = dataset_merged, method = RPCAIntegration,
orig.reduction = "pca", new.reduction = "integrated.rpca",
verbose = TRUE
)
8. データ統合後に結合されたデータセットに対してクラスタリングと UMAP 分析を実行します。
   FindNeighbors(dataset_merged, dims = 1:10, reduction = "integrated.rpca")
dataset_merged <- FindClusters(dataset_merged, resolution = .1, cluster.name = "rpca_clusters")
dataset_merged <- RunUMAP(dataset_merged, dims = 1:10, seed.use = 123, reduction = "integrated.rpca", reduction.name = "umap.rpca")
9. 積分後のクラスター番号とバッチ番号に従ってUMAPプロットを視覚化します(補足図27)。
   DimPlot(dataset_merged, reduction = "umap.rpca", group.by = c("seurat_clusters","batch"))
10. 統合後のバッチ番号に従って、各クラスターのセル番号の分布を表示します。
    table(dataset_merged$batch, dataset_merged$seurat_clusters)
    UMAP プロットと統合データの表から、異なるクラスターにわたる 2 つのバッチ間に優れた重複があります。興味深いことに、データの統合により、同じクラスタリングパラメータを使用して追加のクラスタが識別されました。
11. データセットの統合後、ダウンストリーム解析の前に、マージされたデータセットのレイヤーを結合する必要があります。
    dataset_merged <- JoinLayers(dataset_merged)
12. データセット Seurat オブジェクト を RDS ファイルとして作業ディレクトリに保存します。
    saveRDS(dataset_merged, "merged_dataset_post_Method7.rds")

Results

メソッド#2から始めて、プロトコルは、単一細胞創傷治癒データセットにロードして品質管理ステップを実行するための手順を順を追って説明します。Seuratオブジェクトを作成した後(ステップ2.6.2)、一連のステップでデータセット内の2つのアッセイ(RNAとタンパク質、ステップ2.6.3-2.6.7)をマージし、時空間バーコードに従ってタンパク質アッセイの脱複合体化を実行します(ステップ2.6.8-2.6.9)。解複雑関数は、各セルの最も可能性の高い時空間バーコードを識別する「barcodes_maxID」を含む、データセット内の各セルにいくつかのメタデータラベルを割り当てます(ステップ2.6.10)。ステップ2.6.11では、バイオリンプロット機能を実行し、多重化されたバーコードに基づいて細胞内の検出遺伝子の分布を視覚化します。このステップの代表的な結果(補足図1)は、各バーコードで検出された遺伝子がかなり均一に分布していることを示しており、これはデータセットの整合性と創傷治癒時点の下流分析にとって重要です。タンパク質バーコードに適切なラベルを割り当てた後(ステップ2.6.12)、プロトコルは、各細胞内のミトコンドリア遺伝子の割合の計算(ステップ2.9)から始めて、データセットのRNAアッセイで品質管理ステップを実行する方法を示します。ステップ2.10では、特徴散布図機能を実行して、検出された遺伝子の分布、RNAの数、およびすべての細胞におけるミトコンドリアの割合を視覚化します。このステップの代表的な結果(補足図2)は、ミトコンドリア含有量が大きい細胞が多数存在し、それがRNA数の低下と相関し、死んだ細胞または死にかけている細胞を同定することを示しています。RNA数が少なく、ミトコンドリア含有量が多い細胞を除去した後(ステップ2.11)、ステップ2.12でサブセットデータセットに対して別の特徴散布図関数を実行し、このステップの代表的な結果(補足図3)は、検出された遺伝子の分布と細胞あたりのミトコンドリアRNAの割合がより正常になったことを示し、堅牢な下流分析への道を切り開きます。次に、プロトコルは、データセット内のダブレットの可能性を識別するためのscDblFinder関数の使用を説明し、「scDblFinder.score」と呼ばれる新しいメタデータを各セルに割り当てます(ステップ2.13〜2.14)。ステップ2.15では、バイオリンプロット関数を実行してデータセット内のダブレットスコアの分布を可視化し、このステップの代表的な結果(補足図4)は、比較的高いダブレットスコアを持つセルが多数存在し、0.25が自然なカットオフであり、それを超えるとダブレットの集団が存在するように見えることが示されています。したがって、次の手順では、このパラメーターを使用してデータセットをカットオフの下のセルにサブセット化し(ステップ2.16)、それによってこの単一セルデータセットの品質管理ステップを完了します。

メソッド#3から始めて、プロトコルは、Seuratパッケージとワークフローを使用して、品質管理された単一細胞創傷治癒データセットを分析するための手順を順を追って説明します。RNAデータの正規化とスケーリングの後、PCA解析が実行されます(ステップ3.1)。ステップ3.2では、エルボプロット関数を使用して、最初の50個のPCA次元に関するデータセットの変動量を視覚化し、このステップの代表的な結果(補足図5)は、主要な変動の多くが、グラフの曲がりで識別される最初の13次元内で発生することを示しています。次に、プロトコルは、最初の13のPCA次元と比較的低いクラスタリング解像度パラメータ0.1を使用して、ネイバーを見つけ、データセットの細胞クラスタリング(ステップ3.3)とUMAP次元削減(ステップ3.4)を実行する方法を示しています。ステップ3.5では、UMAPプロット上のセルのクラスタリングを視覚化するために次元プロット機能を実行し、このステップの代表的な結果(図 1)は、データセット内のすべてのセルが、Windows(左)とMacOS(右)を実行しているコンピューターから取得されたUMAPプロットがわずかに異なる8つの主要な色分けされたSeuratクラスターグループを中心にクラスター化されていることを示しています。ステップ3.6では、別の次元プロット機能を実行して、細胞の創傷時間/空間アノテーションを視覚化し、このステップの代表的な結果(図 2)は、データセット内のすべてのセルが時間/空間の原点に従って分散されており、時間/空間アノテーションに従って明らかなクラスタリングがないことを示しています。次に、プロトコルでは、差次的に発現した遺伝子のリストを取得してテキストファイルに保存する方法(ステップ3.8)、スプレッドシートでデータテーブルを開き、各細胞クラスターの上位ランクのクラスターマーカーを取得するためにさまざまなフィルタリングステップを実行する方法について説明します(ステップ3.9〜3.10.6)。これらのステップの代表的な結果 (補足表1)は、ランク付けされた差次発現遺伝子の完全な出力を含む最終的なスプレッドシートファイルであり、別の代表的な結果(補足表2)は、各Seuratクラスターの上位5つのアップレギュレーションおよび発現遺伝子を示す簡略化された表です。次に、プロトコルは、EnrichRと呼ばれるWebベースの機能濃縮分析ツールを使用して、上位のクラスターマーカー遺伝子に従って推定細胞型を同定する方法(ステップ3.11〜3.12)、およびこれらのステップの代表的な結果(図 3)は、8つの細胞クラスターのそれぞれについて、上位に濃縮された細胞タイプを示すEnrichR出力の切り取られたスクリーンショットです。次に、プロトコルは、最も濃縮された細胞タイプの注釈に従って、それぞれのSeuratクラスター内のすべての細胞に「cell_types」と呼ばれる新しいメタデータラベルを割り当てます(ステップ3.14)。ステップ3.15では、次元プロット機能を実行して、名前が変更されたセルクラスターをUMAPプロット上のセルタイプアノテーションとして視覚化し、このステップの代表的な結果(図4)は、データセット内のすべてのセルが主要な色分けされたセルタイプの周りにクラスター化していることを示しました。ステップ3.16では、特徴プロット機能を使用して、一連のUMAPプロット上の上位クラスターマーカー遺伝子(補足表2から)の局在化を視覚化し、代表的な結果(図5)は、それぞれの主要な細胞型クラスター位置内での上位細胞マーカー遺伝子の高い発現を示すUMAPプロットのグリッドです。ステップ3.17および3.18では、ドットプロット機能を実行して、細胞内のトップクラスターマーカー遺伝子の相対的な発現レベルを視覚化し、最初に元のSeuratクラスター番号でグループ化し(ステップ3.17)、次に注釈付き細胞型ラベルによってグループ化しました(ステップ3.18)。これらのステップの代表的な結果は、それぞれのSeuratクラスターでのみトップ細胞マーカー遺伝子の高レベルの発現を確認しました(補足図6)とそれぞれの主要な細胞型(図 6).プロトコルの次のステップでは、元の時空間タンパク質ベースのラベルを厳密に時間的注釈に簡素化し、細胞が発生した創傷後日数(DPW)に基づいて細胞を識別します。ステップ3.20では、次元プロット機能を実行して、細胞をUMAPプロット上のDPWアノテーションとして視覚化し、このステップの代表的な結果(補足図7)は、単一細胞創傷治癒データセット全体にわたる創傷経時アノテーションの局在化を示しました。予想通り、1日目(D1)の注釈は好中球とマクロファージのクラスターを支配していましたが、その後の創傷治癒の時点は他の細胞型でより多く表されました。プロトコルの次のステップでは、積み上げ棒グラフを使用して、最初に異なる細胞型にわたるDPWの比率を視覚化し(ステップ3.22)、次に異なる時点にわたる細胞型の比率を視覚化しました(3.23)。これらのステップの代表的な結果は、各主要な細胞タイプカテゴリ(補足図8)および各DPWカテゴリの主要な細胞型の相対数(図 7).これらの結果は、炎症期の初期の時点で免疫細胞(好中球とマクロファージ)が優勢であり、他の細胞型(上皮細胞と内皮細胞)が増殖期に現れ始め、線維芽細胞が創傷解消の後の時点で特に優勢であるという、皮膚創傷治癒の既知の細胞カスケードを裏付けました。

方法#4から始めて、プロトコルは、創傷治癒中の潜在的な細胞サブタイプを特定するために、Seuratを使用して単一細胞データセット内の個々の主要な細胞型に焦点を当てる手順を概説します。このプロトコルは、最初に2つのSeuratクラスターにクラスター化されてから1つのカテゴリに結合された線維芽細胞に焦点を当て、元のデータセットの線維芽細胞のみを含む新しいSeuratオブジェクトを作成する方法を説明しています(ステップ4.1)。Seuratワークフローは、この線維芽細胞固有のデータセット(ステップ4.2〜4.4)で実行され、ステップ4.2ではエルボープロット(補足図9)が得られ、線維芽細胞データセットの主要な変動の多くが最初の9つのPCA次元内で発生することを示しています。ステップ4.5では、UMAPプロット上の細胞のクラスタリングを視覚化するために次元プロット機能が実行され、このステップの代表的な結果(図8)は、データセット内の線維芽細胞が3つの色分けされた細胞サブタイプの周りにクラスター化されていることを示しました。DPWアノテーション(ステップ4.6)に従って線維芽細胞データセットを視覚化すると、DPWアノテーションに従ってデータセット内の線維芽細胞が全体に分布していることを示すUMAPプロット(補足図10)が得られました。次に、プロトコルは、差次的に発現した遺伝子のリストを取得してテキストファイルに保存する方法(ステップ4.7)、Excelでデータテーブルを開き、各細胞クラスターの上位ランクのクラスターマーカーを取得するためにさまざまなフィルタリングステップを実行する方法(ステップ4.8)を説明し、「FB_type_marker」という名前の上位線維芽細胞マーカー遺伝子をリストする新しい変数を割り当てる方法(ステップ4.9)について説明します。ステップ4.10では、ドットプロット関数を使用して、featuresパラメータの「FB_type_marker」変数を呼び出すことで、線維芽細胞のみのデータセットのリスト内の遺伝子を視覚化し、このステップの代表的な結果(図9)は、線維芽細胞サブタイプマーカーの発現が高いことをそれぞれのクラスターカテゴリでのみ確認するドットプロット(上)ですが、DPWカテゴリ全体にかなり分布しています(下)。ステップ4.11では、創傷治癒データセット全体の線維芽細胞マーカー遺伝子を可視化するために同じ特徴変数を呼び出し、代表的な結果(補足図11)は、主に元の線維芽細胞における線維芽細胞サブタイプマーカーの高発現を確認したドットプロットです。最後に、プロトコルの次のステップでは、積み上げ棒グラフを使用して、最初に3つの線維芽細胞サブタイプにわたるDPWの割合を視覚化し(ステップ4.12)、次に異なる時点にわたる線維芽細胞サブタイプの比率を視覚化しました(ステップ4.13)。これらのステップの代表的な結果は、各線維芽細胞サブタイプカテゴリーにおけるDPW細胞の相対数(補足図12)および各DPWカテゴリーにおける線維芽細胞サブタイプの相対数を示す比率プロットである(補足図13)。これらの結果は、治癒の経時的に線維芽細胞サブタイプの比率に有意な変化があることを示しており、最初の線維芽細胞サブタイプ(クラスター0)は初期段階の創傷(D1およびD3)で強く優勢であり、2番目のサブタイプ(クラスター1)は創傷治癒の増殖期(D14)で優勢であり、3番目のサブタイプ(クラスター2)は創傷治癒の増殖期(D7)で最も高かった。

メソッド #5 から始めて、プロトコルは、Seurat のモジュール スコアリング関数を使用して単一細胞創傷治癒データセットを分析する手順を順を追って説明します。このプロトコルは、最初にタブ削除されたテキストファイルを使用して遺伝子セットをRの変数にアップロードする手順(ステップ5.1〜5.2)を説明し、続いて、創傷治癒の3つの主要な段階に関連する3つの遺伝子セットにモジュールスコアリング関数を適用する(ステップ5.3)。ステップ5.4では、ドットプロット機能を使用して、2つの異なるメタデータカテゴリにわたる集計モジュールスコアを視覚化し、このステップの代表的な結果(図10)は、創傷後日数カテゴリ(DPW、右)および主要な細胞タイプカテゴリ(左)の細胞全体の主要な治癒期モジュールの平均発現を示すドットプロットです。これらの結果は、バルクシーケンシングベースの遺伝子発現プロファイルを単一細胞発現データセットに擬似バルク方式で適用することが、創傷治癒の分野で以前に公開されたデータセットを活用した比較バイオインフォマティクスアプローチの強力な方法であることを示しています。

メソッド#6から始めて、プロトコルは、初期と後期の創傷に由来する細胞を比較するという特定の科学的問題に従って、CellChatパッケージとワークフローを使用してSeurat由来の単一細胞創傷治癒データセットを分析するための手順を順を追って説明します。プロトコルはまず、Seuratデータセット全体を損傷後の2つの時点、1つは炎症期(1日目(D1))ともう1つは創傷解消中(14日目(D14))の2つの時点にサブセットします(ステップ6.1)。2つのCellChatオブジェクトが作成され、プロトコルはCellChatプロトコルのすべての典型的な機能を実行して、プロトコルの方法#3(ステップ6.2〜6.3)で特定された細胞型間の推定相互作用をすべて計算します。ステップ6.4では、シグナル伝達散布図機能を実行して、各創傷治癒時点におけるすべての主要な細胞型における入出力相互作用の強さを視覚化します。このステップの代表的な結果 (補足図14)は、D1(左)とD14(右)の時点での主要細胞型の入ってくる(y軸)と出る(x軸)の相互作用の強さを示す散布図です。これらの結果は、好中球やマクロファージなどの免疫細胞が炎症期に最も細胞間相互作用の強さを持っているが、創傷解消期には線維芽細胞が細胞間相互作用を支配していることを示し、数十年にわたる創傷治癒研究を裏付けています。次のステップでは、大幅に濃縮された経路の1つであるコラーゲン経路(ステップ6.5〜6.6)に焦点を当てて分析します。ステップ6.7では、円図機能を実行して、2つの時点での細胞型間のコラーゲンシグナル伝達経路の相互作用を視覚化します。このステップの代表的な結果 (補足図15)は、D1(左)とD14(右)ですべての細胞型間の推定されたコラーゲン経路シグナル伝達相互作用を示す円プロットです。ステップ6.8では、コードダイアグラム関数を使用して同じ相互作用を視覚化し、代表的な結果(補足図16)は、各時点でのすべての細胞型間の推定されたコラーゲン経路シグナル伝達相互作用を示すコード図です。予想通り、これらの結果は、線維芽細胞がコラーゲンシグナル伝達経路の主要なソース細胞であることを示しましたが、情報の流れはD14と比較してD1の免疫細胞に制限されていました。細胞間相互作用のソース細胞としての線維芽細胞に焦点を当てるために、ステップ6.9は、ソースセルパラメータを追加するコードダイアグラム関数を繰り返し、代表的な結果(補足図17)は、各時点で線維芽細胞をソース細胞として行う推定コラーゲン経路シグナル伝達相互作用を示すコード図です。ステップ6.10では、線維芽細胞をソース細胞としてコラーゲンシグナル伝達経路における各リガンド-受容体ペアの寄与を視覚化するために、1つはバブルプロット(ステップ6.10.1)を使用し、もう1つはコードダイアグラム(ステップ6.10.2)を使用して視覚化するために2つの機能を実行します。代表的な結果は、両方のバブルプロット(補足図18)とコードダイアグラム(補足図19).これらの結果は、D1では線維芽細胞由来のコラーゲン経路がCd44およびSdc4受容体が優勢な好中球とマクロファージに限定されていたが、D14では他の細胞がインテグリンを含むさまざまな受容体を介してレシーバーとして機能することを示しました。線維芽細胞相互作用において強い強みを示したCol1a1-Cd44リガンド-受容体相互作用に焦点を当てるために、パラメータを設定し(ステップ6.11)、ステップ6.12でコードダイアグラム関数を使用して、すべての細胞型間のこの特定のリガンド-受容体相互作用を視覚化し、代表的な結果(補足図20)は、D1(左)とD14(右)の時点でのすべての細胞型間の推定されたCol1a1-Cd44リガンド-受容体相互作用を示すコード図です。これらの結果は、D1ではこの相互作用がソース細胞としての線維芽細胞に限定されているのに対し、D14ではマクロファージと平滑筋細胞もソース細胞として機能することを示しました。次に、このプロトコルでは、最初に D1 と D14 の CellChat オブジェクトをマージすることにより、差分 CellChat 分析を実行する方法を記述します (ステップ 6.13)。ステップ6.14では、相互作用の比較機能を実行して、2つの創傷治癒時点間の細胞間相互作用の総数と相対的な強さを視覚化し、代表的な結果(補足図21)は、D1とD14の創傷を含む細胞における推定された相互作用の総数(左)と強度(右)を示す結果の棒グラフであり、D1の相互作用の相対的な強さが高いのとは対照的に、D14の相互作用の数が多い。ステップ6.15および6.16では、2つの関数を使用して、創傷が1日目から14日目に移行する際の各細胞型間の異なる細胞間相互作用の強さを、それぞれの代表的な結果とともに視覚化し、最初の関数は円プロットです(ステップ6.15、 補足図22)で、2つ目はヒートマップ(ステップ6.16、 補足図23)で、D1と比較してD14で増加した相互作用は赤で示され、減少した相互作用は青で示されています。予想通り、好中球とマクロファージを介した相互作用はD1で増加し、線維芽細胞を介した相互作用はD14で増加します。ステップ6.17では、ランク付け関数を使用して、D1と比較してD14で線維芽細胞をソース細胞とする細胞間相互作用に対する個々の経路の相対的な寄与をランク付けするプロットを作成し、代表的な結果(補足図24)は、D1が上部に赤で、D14が下部に青で表され、いくつかの経路がD1またはD14のいずれかで排他的に表され、他の多くの経路が活性化の勾配を示す結果のランクプロットを示しています。最後に、ステップ6.18では、2つのバブルプロット関数を使用して、D1と比較してD14で線維芽細胞をソース細胞とするコラーゲンシグナル伝達経路における個々のリガンド-受容体ペアの相対的な寄与を示し、対応する代表的な結果(補足図25)は、X軸上の多くの細胞間相互作用にわたって、D1と比較してD14のシグナル伝達ペアの増加(左)と減少(右)を示しています。予想通り、線維芽細胞は、炎症期にコミュニケーションが好中球とマクロファージに限定されていたD1創傷と比較して、D14創傷のいくつかのレシーバー細胞間で発信リガンド-受容体ペアの相互作用がはるかに増加しました。

メソッド #7 から始めて、プロトコルは Seurat を使用して 2 つの単一細胞創傷治癒データセットを統合する手順を順を追って説明します。このプロトコルは、最初に、公開された単一細胞データセットの2つのバッチをマージし、マージされたデータセットに標準のSeuratワークフローを適用する手順を記述します(ステップ7.1〜7.4)。ステップ7.5では、次元プロット機能を使用して、マージされたがまだ統合されていない創傷治癒データセットのクラスター番号とバッチ番号に従ってUMAPプロットを視覚化します。このステップの代表的な結果(補足図26)は、Seuratクラスター(左)とバッチ番号(右)の分布を可視化したUMAPプロットであり、データ統合前のこれら2つのデータセットに有意なバッチ効果は見られないことを示しています。次に、プロトコルは、RPCA法と統合データセットのフォローアップSeuratワークフローを使用してデータ統合を実行します(ステップ7.7〜7.8)。ステップ7.9では、次元プロット機能を使用して、統合された創傷治癒データセットのクラスター番号とバッチ番号に従ってUMAPプロットを視覚化します。このステップの代表的な結果(補足図27)は、Seuratクラスター(左)とバッチ番号(右)の分布を視覚化したUMAPプロットであり、異なるクラスター間で2つのバッチ間の重複がさらに大きくなったことを示しています。結果はまた、データの統合後に追加のクラスターが出現していることを示しており、これは、データバッチの技術的効果が制御された後に、潜在的に重要な細胞サブタイプを特定する能力が向上したことを示している可能性があります。

図1:データセット内のすべてのセルが8つの主要な色分けされたクラスターグループを中心にクラスター化したことを示すUMAPプロット。 Windows(左)とMacOS(右)を実行しているコンピューターから取得した結果。この図は、ステップ3.5に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

図2:データセット内のすべてのセルが時間/空間の原点に従って広がっていることを示すUMAPプロットで、時間/空間アノテーションによる明らかなクラスタリングはありません。 この図は、ステップ3.6に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

図3:EnrichR出力のトリミングされたスクリーンショットで、各細胞クラスターの上位に濃縮された細胞タイプを示しています。 この図は、ステップ3.13に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

図4:主要な色分けされた細胞型の周りにクラスター化されたデータセット内のすべての細胞を示すUMAPプロット。 この図はステップ3.15に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

図5:主要な細胞型クラスター内の上位細胞マーカー遺伝子の高い発現を示すUMAPプロットのグリッド。 この図は、ステップ3.16に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

図6:それぞれの主要な細胞型でのみ、上位の細胞マーカー遺伝子の高レベルの発現を確認するドットプロット。 この図は、ステップ3.18に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

図7:各DPWカテゴリーの主要細胞型の相対数を示す比率プロット。 この図はステップ3.23に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

図8:3つの色分けされた細胞サブタイプの周りにクラスター化されたデータセット内の線維芽細胞を示すUMAPプロット。 この図は、ステップ4.5に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

図9:線維芽細胞サブタイプマーカーの高発現をそれぞれのクラスターカテゴリーでのみ確認するドットプロットですが、DPWカテゴリー全体に公平に分布しています。 この図はステップ4.10に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

図10:DPWごとおよび主要な細胞タイプごとの細胞全体の主要な治癒期モジュールの平均発現を示すドットプロット。 この図は、ステップ5.4に対応します。 この図の拡大版を表示するには、ここをクリックしてください。

補足図1:各バーコードで検出された遺伝子がかなり均等に分布していることを示す結果は、データセットの整合性と創傷治癒時点の下流分析にとって重要です。 この図は、手順 2.6.11 に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図2:ミトコンドリア含有量が多い細胞が多数存在することを示す散布図であり、これらは死んだ細胞または死にかけている細胞である---、RNA数の低下と相関しています。 この図は、ステップ2.10に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図3:検出された遺伝子の分布と細胞あたりのミトコンドリアRNAの割合がより正常になったことを示す散布図は、堅牢な下流分析への道を切り開きました。 この図はステップ2.12に対応します。この図をダウンロードするには、ここをクリックしてください。

補足図4:比較的高いダブレットスコアを持つセルが多数あり、0.25が自然なカットオフであるように見え、それを超えるとダブレットの可能性のある集団があることを示すバイオリンプロット。 この図はステップ2.15に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図5:主要な変動の多くが最初の13次元内で発生することを示すエルボプロット。 この図は、ステップ3.2に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図6:それぞれのSeuratクラスターでのみ、上位細胞マーカー遺伝子の高レベルの発現を確認するドットプロット。 この図はステップ3.17に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図7:創傷治癒データセット全体の創傷経時アノテーションの局在化を示すUMAPプロット。 この図はステップ3.20に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図8:各主要な細胞タイプカテゴリにおけるDPW細胞の相対数を示す比率プロット。 この図は、ステップ3.22に対応します。この図をダウンロードするには、ここをクリックしてください。

補足図9:線維芽細胞データセットの主要な変動の多くが最初の9次元内で発生することを示すエルボプロット。 この図は、ステップ4.2に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図10:DPWアノテーションに従って全体に分布したデータセット内の線維芽細胞を示すUMAPプロット。 この図は、ステップ4.6に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図11:主に元の線維芽細胞クラスターにおける線維芽細胞サブタイプマーカーの高い発現を確認するドットプロット。 この図は、ステップ4.11に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図12:各DPWカテゴリにおける線維芽細胞サブタイプの相対数を示す比率プロット。 この図はステップ4.12に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図13:各線維芽細胞サブタイプカテゴリにおけるDPW全体の線維芽細胞の相対数を示す比率プロット。 この図はステップ4.13に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図14:1日目(D1、左)および14日目(D14、右)の時点での主要細胞型の入ってくる(y軸)と出る(x軸)相互作用の強さを示す散布図。 この図は、ステップ6.4に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図15:各DPWカテゴリーのすべての細胞型間の推定されたコラーゲン経路シグナル伝達相互作用を示す円プロット。 この図は、ステップ6.7に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図16:各DPWカテゴリーのすべての細胞型間の推定されたコラーゲン経路シグナル伝達相互作用を示す弦図。 この図は、ステップ 6.8 に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図17:各DPWカテゴリーのソース細胞としての線維芽細胞との推定されたコラーゲン経路シグナル伝達相互作用を示す弦図。 この図は、ステップ6.9に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図18:各DPWカテゴリーのソース細胞として線維芽細胞を使用したコラーゲン経路シグナル伝達における各リガンド-受容体ペアの推定寄与を示すバブルプロット。 この図は、手順6.10.1に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図19:各DPWカテゴリーのソース細胞として線維芽細胞を使用したコラーゲン経路シグナル伝達における各リガンド-受容体ペアの推定寄与を示すコード図。 この図は、手順 6.10.2 に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図20:各DPWカテゴリーのすべての細胞型間の推定されたCol1a1-Cd44リガンド-受容体相互作用を示す弦図。 この図は、ステップ6.12に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図21:1日目と14日目の創傷における推測された相互作用の数(左)と強度(右)を示す棒グラフ。 この図は、ステップ6.14に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図22:創傷が1日目(青)から14日目(赤)のDPWに移行する際の、各細胞型間の細胞間相互作用の強さの違いを示す円プロット。 この図はステップ6.15に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図23:創傷が1日目(青)から14日目(赤)に移行する際の各細胞型間の異なる細胞間相互作用の強さを示すヒートマップ。 この図は、ステップ6.16に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図24:1日目と14日目のDPWにおける線維芽細胞と他の細胞型の間の細胞間相互作用に対する個々の経路の相対的な寄与を示す順位プロット。 この図は、ステップ6.17に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図25:1日目と14日目のDPWで線維芽細胞をソース細胞とするコラーゲンシグナル伝達経路における個々のリガンド-受容体ペアの相対的な寄与を示すバブルプロット。 この図は、ステップ6.18に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図26:データ統合前のSeuratクラスター(左)とバッチ番号(右)の分布を示すUMAPプロット。 この図は、ステップ 7.5 に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足図27:データ統合後のSeuratクラスター(左)とバッチ番号(右)の分布を示すUMAPプロット。 この図は、ステップ7.9に対応します。 この図をダウンロードするには、ここをクリックしてください。

補足ファイル1:JoVE_Rscript.R: メイン R コード スクリプト ファイルには、プロトコルのすべての部分について説明されているすべての手順と説明が含まれています。 このファイルをダウンロードするには、ここをクリックしてください。

補足ファイル2:JoVE_PhaseSpecificGenes.txt。 プロトコルのステップ5.1でロードされた遺伝子のリストを含むタブ区切りのテキストファイル。 このファイルをダウンロードするには、ここをクリックしてください。

補足ファイル3:JoVE_Rscript_b3.R. プロトコルのステップ 7.1 で使用するデータセットのバッチ #3 を分析するために必要なすべての手順と説明を含む、補足的な R コード スクリプト ファイル。 このファイルをダウンロードするには、ここをクリックしてください。

補足表1:JoVE_DEGs_cellMarkers.xlsx。 プロトコルのステップ3.10で使用されたランク付けされた差次的に発現した遺伝子の完全な出力を含むExcelファイル。 この表をダウンロードするには、ここをクリックしてください。

補足表2:各スーラクラスターの上位5つのアップレギュレーションおよび発現遺伝子。この表をダウンロードするには、ここをクリックしてください。

Discussion

このプロトコルでは、RStudioを使用して、Seuratを使用して複雑な単一セルデータセットの基本的な分析を可能にする所定のコード行を実行します。Rコーディング環境のインストール、以前に公開された単一細胞創傷治癒データセットのダウンロード、重要な品質管理ステップの実行、および視覚化、主要な細胞型注釈、細胞サブタイプ分析、およびSeuratを使用した統合分析を含む標準的な単一細胞分析ワークフローの実行、およびCellChatを使用した細胞間相互作用分析の実行など、創傷治癒研究に関連するいくつかの方法が提示されています。

ここで紹介する方法は、Rとその人気のあるオープンソースの科学パッケージであるSeurat²¹ およびCellChat²²を使用した単一細胞分析の典型的なワークフローの簡略化されたビネットです。実際、このワークフローは、複雑な単一細胞創傷治癒データセットで達成できる分析タイプの一例にすぎません。この方法に可能な変更はほぼ無限であり、唯一の制約はユーザーの特定の科学的調査です。たとえば、ユーザーは、このデータセットに尋ねたい研究課題に応じて、細胞の種類や時点などの主要なパラメーターの一部を変更できます。著者らはまた、ユーザーがこのワークフローを関心のある独自の単一セルデータセットに適応させるのに十分快適であることを望んでいます。ただし、このワークフローを使用して他のデータセットを分析する場合は、各実験が技術的な問題やサンプル前処理の問題をデータ自体に伝播する可能性があるため、注意が必要です。したがって、以前に公開および再分析された単一セルデータセットの結果を解釈する前に、ユーザーはすべての実験の詳細を読んで理解することが不可欠です。バイオインフォマティクス ツールは、生物学的プロセスを探索し、仮説を立てるための強力な方法であり、結果の重要な生物学的解釈は、フォローアップ実験で検証する必要があることを覚えておくことが重要です。

プロトコル全体を通して、他のタスクを実行するために、ワークフローの特定の領域に大きな変更が加えられる可能性があることに注意してください。ただし、ワークフローに対する変更の可能なすべての組み合わせの詳細は、この原稿の範囲を超えています。たとえば、細胞クラスタリングに使用される解像度とUMAP分析に使用される次元は必然的に主観的であり、ここで紹介するツールは、大規模な分析(広義に定義された主要な細胞タイプについてここで示されたように)と、より大きなデータセット内のよりまれな亜集団への細胞のサブクラスタリングを伴う可能性のある非常に特殊な分析の両方を可能にします。単一細胞解析法のこの側面の詳細、および単一細胞解析パイプラインで変更される可能性のある他のすべてのパラメーターの詳細については、著者らは、この進化するツールの著者が詳細な説明、ビネット、およびチュートリアルを提供するSeuratの出版物^21,26およびWebサイト(https://satijalab.org/seurat/)を参照します。

この原稿では、単一細胞トランスクリプトミクスの文献で最も引用され、使用されているツールのいくつか、すなわち、それぞれ単一細胞および細胞間相互作用分析のためのSeurat²¹ とCellChat²²を紹介しました。ただし、わずかに異なる方法で同様の機能を実行する他のツールも存在します。単一セルデータセット分析には、Scran²⁷、Scater ²⁸、および Python ベースの ScanPy²⁹ があり、データセット統合にさまざまな方法を使用します²⁵。このプロトコルでは、クラスター細胞マーカーの濃縮を解釈するためのユーザーの判断に依存する細胞タイプの手動アノテーションが実証されましたが、現在では、SingleR³⁰ やscGate³¹など、細胞タイプの自動分類を可能にするさまざまなツールが存在します。細胞間コミュニケーション解析については、CellChatがこのプロトコルで実証されましたが、CellPhoneDB³²、Cytotalk³³、およびLIANA(LIgand-receptor ANalysisフレームワーク)コンセンサスフレームワーク³⁴内に実装されている他のリガンド受容体データベースなど、細胞間コミュニケーションを推定するための他のツールが存在します。すべてのバイオインフォマティクスツールはユニークであり、独自の特性と変更可能なパラメータを備えています。したがって、ユーザーは、使用から生成された出力を解釈する前に、各ツールの関連ドキュメントを注意深く読んで、そのニュアンスを理解することが重要です。最後に、どのようなバイオインフォマティクス ツールを使用する場合でも、そのようなツールは継続的に進化しており、パッケージのバージョンが異なれば出力も異なる可能性があることを覚えておくことが重要です。

R では構文が重要であり、句読点、引用符、括弧、または大文字の区別が間違っているとエラーが発生します。したがって、ユーザーがコードを入力するときに細部に注意を払い、コード行をコピーするときは、新しい科学的な質問やデータセットに適応させるために特に注意を払うことが重要です。発生する可能性のある特定のエラーのトラブルシューティングのために、著者は、最も一般的に発生するエラーは知識のあるパワーユーザーによってすでに回答されている可能性が高いため、エラーメッセージをユーザーのお気に入りのWeb検索エンジンにコピーして貼り付け、GitHubやStack Overflowなどのバイオインフォマティクスフォーラムの結果を参照することを推奨しています。一部のフォーラムでは、最も成功した回答は、問題に最適な解決策を見つけた他のユーザーによって「賛成票」が投じられます。ユーザーは、インターネットで見つけたコード行を自分のコンピューターにコピーして貼り付けないように注意する必要があります (特に、ソリューションで R プログラミング言語の外部でシステム設定を変更する必要がある場合)、そのようなプログラムが悪意のあるものである可能性があるためです。コーディング エラーをトラブルシューティングするエキサイティングな新しい方法は、OpenAI の ChatGPT、Microsoft の Copilot、Google の Gemini などの強力な生成大規模言語 AI モデルを使用することです。これらのモデルは、ソフトウェアエンジニアリング全般、特にトラブルシューティングに特に役立つことが証明されています。このために、ユーザーは、コードに対するユーザーの意図についてチャットボットに簡単なプロンプトを提供した後、コードの全行をコピーして貼り付けることができます。これらのモデルは絶対確実ではなく、問題の解決に適した回答を生成するために、ユーザーは複数のプロンプトを試す必要がある場合があるという通常の注意点が存在します。

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
ラップトップまたはデスクトップコンピューター	該当なし	該当なし	WindowsまたはMacOS を実行している。
R	該当なし	バージョン 4.4.1	https://cran.rstudio.com/ から無料でダウンロード可能
Rstudio スタジオ	Positソフトウェア、PBC	バージョン 2024.09.0	https://posit.co/download/rstudio-desktop/ から無料でダウンロード
オフィスエクセル	マイクロソフト	任意のバージョン	テーブルデータの分析用
インターネットブラウザ	該当なし	該当なし	Web サイトへの移動
R パッケージ	リポジトリ	バージョン
開発ツール	クラン	2.4.5
読み取りXL	クラン	1.4.3
OpenXLSX	クラン	4.2.7.1
ティディバース	クラン	2.0.0
scカスタマイズ	クラン	2.1.2
バイオマネージャー	生体伝導体	1.30.25
NMFの	生体伝導体	0.28
コンプレックスヒートマップ	生体伝導体	2.20.0
バイオックネイバー	生体伝導体	1.22.0
シングルセル実験	生体伝導体	1.26.0
循環	生体伝導体	0.4.16
エッジR	生体伝導体	4.2.1
scDblファインダー	生体伝導体	1.18.0
スーラ	クラン	5.1.0
セルチャット	Github(英語)	2.1.2