教程背景
這個(gè)例子演示了允許用戶使用Seurat分析和探索多模式數(shù)據(jù)的新特性。雖然這只是一個(gè)初始版本，但我們很高興在未來(lái)為多模態(tài)數(shù)據(jù)集發(fā)布重要的新功能。
在這里，我們分析了由CITE-seq產(chǎn)生的8,617個(gè)臍帶血單核細(xì)胞(CBMCs)的數(shù)據(jù)集，我們同時(shí)測(cè)量了單細(xì)胞轉(zhuǎn)錄組和11個(gè)表面蛋白的表達(dá)，其水平通過(guò)DNA-barcode抗體來(lái)量化。首先，我們加載兩個(gè)計(jì)數(shù)矩陣:一個(gè)用于RNA測(cè)量，另一個(gè)用于抗體衍生標(biāo)簽ADT。你可以在這里下載ADT文件和RNA文件
1.加載RNA UMI矩陣
注意，這個(gè)數(shù)據(jù)集還包含了大約5%的小鼠細(xì)胞，我們可以將其作為蛋白質(zhì)測(cè)量的陰性對(duì)照。因此，基因表達(dá)矩陣在每個(gè)基因的開(kāi)頭都附加了HUMAN_或MOUSE_。
cbmc.rna <- as.sparse(read.csv(file = "Analysis/data/GSE100866_CBMC_8K_13AB_10X-RNA_umi.csv.gz", sep = ",", header = TRUE, row.names = 1))
為了讓生活更容易繼續(xù)下去，我們將拋棄除前100個(gè)高度表達(dá)的老鼠基因外的所有基因，并從CITE-seq前綴中去除“HUMAN_”
# 這個(gè)函數(shù)有幾個(gè)默認(rèn)參數(shù)，ncontrols = 100，prefix = "HUMAN_"，controls = "MOUSE_"
cbmc.rna <- CollapseSpeciesExpressionMatrix(cbmc.rna)
2.加載ADT UMI矩陣
cbmc.adt <- as.sparse(read.csv(file = "Analysis/data/GSE100866_CBMC_8K_13AB_10X-ADT_umi.csv.gz", sep = ",", header = TRUE, row.names = 1))
在向Seurat添加多模態(tài)數(shù)據(jù)時(shí)，可以使用重復(fù)的特性名稱(chēng)。每一組模態(tài)數(shù)據(jù)(例如。RNA、ADT等)儲(chǔ)存在自己的Assay對(duì)象中。其中一個(gè)Assay對(duì)象被稱(chēng)為“default assay”，意思是它用于所有的分析和可視化。若要從非默認(rèn)Assay中提取數(shù)據(jù)，可以指定一個(gè)與Assay鏈接的鍵，用于提取特征。要查看所有對(duì)象的所有鍵，請(qǐng)使用Key function。
最后，我們觀察到CCR5、CCR7和CD10的低富集，因此將它們從矩陣中去除(可選)。(關(guān)鍵是這三個(gè)蛋白得低富集怎么看出來(lái)的？)
cbmc.adt <- cbmc.adt[setdiff(rownames(x = cbmc.adt), c("CCR5", "CCR7", "CD10")), ]
3.設(shè)置一個(gè)Seurat對(duì)象，并基于RNA表達(dá)對(duì)細(xì)胞進(jìn)行聚類(lèi)
下面的步驟表示基于scRNA-seq數(shù)據(jù)對(duì)pbmc進(jìn)行快速聚類(lèi)。有關(guān)單個(gè)步驟或更高級(jí)選項(xiàng)的詳細(xì)信息，請(qǐng)參見(jiàn)這里的PBMC集群指導(dǎo)教程。
cbmc <- CreateSeuratObject(counts = cbmc.rna) # 創(chuàng)建seurat對(duì)象
cbmc <- NormalizeData(cbmc) # 標(biāo)準(zhǔn)化
cbmc <- FindVariableFeatures(cbmc) # 選擇高變基因
cbmc <- ScaleData(cbmc) # 歸一化
cbmc <- RunPCA(cbmc, verbose = FALSE) # 運(yùn)行PCA降維
ElbowPlot(cbmc, ndims = 50) # 確定取多少個(gè)PCs
這里教程說(shuō)選擇13個(gè)PCs，但是看代碼選的是25個(gè)，蜜汁迷惑。
cbmc <- FindNeighbors(cbmc, dims = 1:25)
cbmc <- FindClusters(cbmc, resolution = 0.8)
cbmc <- RunTSNE(cbmc, dims = 1:25, method = "FIt-SNE")
# 找到定義每個(gè)cluster的標(biāo)記，并使用它們對(duì)cluster進(jìn)行注釋?zhuān)覀兪褂胢ax.cells.per加快進(jìn)程
cbmc.rna.markers <- FindAllMarkers(cbmc, max.cells.per.ident = 100, min.diff.pct = 0.3, only.pos = TRUE)
head(cbmc.rna.markers)
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
IL7R 4.280959e-21 1.1368139 0.872 0.274 8.776395e-17 0 IL7R
LTB 3.149773e-20 1.0246049 0.985 0.525 6.457351e-16 0 LTB
LDHB 8.557593e-20 0.9786093 0.964 0.527 1.754392e-15 0 LDHB
TRBC2 9.896619e-18 0.7739586 0.868 0.425 2.028906e-13 0 TRBC2
IL32 9.951232e-18 0.9599525 0.896 0.327 2.040102e-13 0 IL32
ITM2A 2.335897e-17 0.8866571 0.786 0.298 4.788822e-13 0 ITM2A
注意，為了簡(jiǎn)單起見(jiàn)，我們將兩個(gè)CD14+單核細(xì)胞簇(HLA-DR基因表達(dá)不同)和NK細(xì)胞簇(細(xì)胞周期階段不同)合并。
這里直接進(jìn)行了細(xì)胞類(lèi)型注釋?zhuān)⑨尩眠^(guò)程？
new.cluster.ids <- c("Memory CD4 T", "CD14+ Mono", "Naive CD4 T", "NK", "CD14+ Mono", "Mouse", "B",
"CD8 T", "CD16+ Mono", "T/Mono doublets", "NK", "CD34+", "Multiplets", "Mouse", "Eryth", "Mk",
"Mouse", "DC", "pDCs")
names(new.cluster.ids) <- levels(cbmc)
cbmc <- RenameIdents(cbmc, new.cluster.ids)
DimPlot(cbmc, label = TRUE) + NoLegend()
4.將蛋白質(zhì)表達(dá)水平添加到Seurat對(duì)象
Seurat v3.0允許您將來(lái)自多個(gè)分析的信息存儲(chǔ)在同一個(gè)對(duì)象中，只要數(shù)據(jù)是多模態(tài)的(在同一組細(xì)胞中收集)。您可以使用SetAssayData和GetAssayData訪問(wèn)器函數(shù)來(lái)添加和從其他分析中獲取數(shù)據(jù)。
我們將定義一個(gè)ADT試驗(yàn)，并為它存儲(chǔ)原始計(jì)數(shù)。
如果你對(duì)這些數(shù)據(jù)內(nèi)部如何存儲(chǔ)感興趣，您可以查看Assay類(lèi)，它在object . R中定義。注意，所有的單細(xì)胞表達(dá)數(shù)據(jù)，包括RNA數(shù)據(jù)，仍然存儲(chǔ)在Assay對(duì)象中，也可以使用GetAssayData訪問(wèn)。
cbmc[["ADT"]] <- CreateAssayObject(counts = cbmc.adt)
現(xiàn)在我們可以重復(fù)我們通常在RNA上運(yùn)行的預(yù)處理(標(biāo)準(zhǔn)化和縮放)步驟，但是需要修改‘a(chǎn)ssay’的參數(shù)。對(duì)于CITE-seq數(shù)據(jù)，我們不建議使用典型的對(duì)數(shù)標(biāo)準(zhǔn)化。相反，我們使用中心對(duì)數(shù)比(CLR)歸一化，獨(dú)立計(jì)算每個(gè)特性。與原始版本相比，這是一個(gè)略微改進(jìn)的過(guò)程，我們將很快發(fā)布更高級(jí)的CITE-seq規(guī)范化版本。
# 這里的標(biāo)準(zhǔn)化方法是CLR
cbmc <- NormalizeData(cbmc, assay = "ADT", normalization.method = "CLR")
cbmc <- ScaleData(cbmc, assay = "ADT")
5.觀察RNA簇上的蛋白水平
您可以使用任何ADT標(biāo)記的名稱(chēng)?！癮dt_CD4”)，在FetchData、FeaturePlot、RidgePlot、FeatureScatter、DoHeatmap或任何其他可視化特性中。
Q1：這里得ADTfeatures為什么突然多了三個(gè)字母adt_？
Q2：這里得蛋白和RNA對(duì)應(yīng)得表達(dá)水平，從下圖可以看出來(lái)什么？
# in this plot, protein (ADT) levels are on top, and RNA levels are on the bottom
FeaturePlot(cbmc, features = c("adt_CD3", "adt_CD11c", "adt_CD8", "adt_CD16", "CD3E", "ITGAX", "CD8A",
"FCGR3A"), min.cutoff = "q05", max.cutoff = "q95", ncol = 4)
在這個(gè)圖中，蛋白質(zhì)(ADT)水平在上面，RNA水平在下面。
山脊圖繪制
RidgePlot(cbmc, features = c("adt_CD3", "adt_CD11c", "adt_CD8", "adt_CD16"), ncol = 2)
繪制ADT散點(diǎn)圖(如FACS的雙軸圖)。注意，如果需要的話，你甚至可以使用HoverLocator和FeatureLocator來(lái)“gate”細(xì)胞。
FeatureScatter(cbmc, feature1 = "adt_CD19", feature2 = "adt_CD3")
查看蛋白質(zhì)和RNA之間的關(guān)系
FeatureScatter(cbmc, feature1 = "adt_CD3", feature2 = "CD3E")
查看T細(xì)胞中CD4和CD8的水平
tcells <- subset(cbmc, idents = c("Naive CD4 T", "Memory CD4 T", "CD8 T"))
FeatureScatter(tcells, feature1 = "adt_CD4", feature2 = "adt_CD8")
讓我們看看原始的(非標(biāo)準(zhǔn)化的)ADT的count值。你可以看到這些值相當(dāng)高，特別是與RNA值相比。這是由于細(xì)胞中蛋白質(zhì)拷貝數(shù)顯著增加，這大大減少了ADT數(shù)據(jù)中的“drop-out”。
對(duì)于drop-out如何理解？
FeatureScatter(tcells, feature1 = "adt_CD4", feature2 = "adt_CD8", slot = "counts")
如果你仔細(xì)觀察，你會(huì)發(fā)現(xiàn)我們的CD8 T細(xì)胞群富含CD8 T細(xì)胞，但仍然包含許多CD4+ CD8- T細(xì)胞。這是因?yàn)槌跏糃D4和CD8 T細(xì)胞在轉(zhuǎn)錄上非常相似，CD4和CD8的RNA丟失水平相當(dāng)高。這說(shuō)明了單獨(dú)從scRNA-seq數(shù)據(jù)中定義細(xì)微的免疫細(xì)胞差異的挑戰(zhàn)。
saveRDS(cbmc, file = "Analysis/cbmc.rds")
6.識(shí)別簇間差異表達(dá)的蛋白
對(duì)cluster進(jìn)行采樣，每個(gè)cluster最多300個(gè)細(xì)胞(使小cluster的熱圖更容易看到)。
cbmc.small <- subset(cbmc, downsample = 300)
找到所有cluster的蛋白質(zhì)標(biāo)記，并繪制一個(gè)熱圖。
adt.markers <- FindAllMarkers(cbmc.small, assay = "ADT", only.pos = TRUE)
DoHeatmap(cbmc.small, features = unique(adt.markers$gene), assay = "ADT", angle = 90) + NoLegend()
你可以看到我們的未知細(xì)胞同時(shí)表達(dá)骨髓和淋巴標(biāo)記(在RNA水平上也是如此)。它們可能是應(yīng)該被丟棄的細(xì)胞團(tuán)塊(多胞胎)。我們現(xiàn)在也要移除老鼠細(xì)胞。
cbmc <- subset(cbmc, idents = c("Multiplets", "Mouse"), invert = TRUE)
7.直接在蛋白質(zhì)水平上聚類(lèi)
還可以直接在CITE-seq數(shù)據(jù)上運(yùn)行降維和基于圖的聚類(lèi)。
#因?yàn)槲覀儗V泛地使用ADT數(shù)據(jù)，所以我們將把默認(rèn)assay改為“ADT”assay。這將導(dǎo)致所有函數(shù)默認(rèn)使用ADT數(shù)據(jù)，而不是要求我們每次都指定它
DefaultAssay(cbmc) <- "ADT"
cbmc <- RunPCA(cbmc, features = rownames(cbmc), reduction.name = "pca_adt", reduction.key = "pca_adt_",
verbose = FALSE)
DimPlot(cbmc, reduction = "pca_adt")
在我們?cè)贏DT級(jí)別上重新聚類(lèi)數(shù)據(jù)之前，我們將稍后隱藏RNA集群id。
cbmc[["rnaClusterID"]] <- Idents(cbmc)
現(xiàn)在，我們僅在ADT(蛋白質(zhì))水平上使用PCA重新計(jì)算tSNE。
cbmc <- RunTSNE(cbmc, dims = 1:10, reduction = "pca_adt", reduction.key = "adtTSNE_", reduction.name = "tsne_adt")
cbmc <- FindNeighbors(cbmc, features = rownames(cbmc), dims = NULL)
cbmc <- FindClusters(cbmc, resolution = 0.2, graph.name = "ADT_snn")
我們可以比較RNA和蛋白質(zhì)的聚類(lèi)，并使用它來(lái)注釋蛋白質(zhì)聚類(lèi)。
(當(dāng)然，我們也可以使用findmarker)
clustering.table <- table(Idents(cbmc), cbmc$rnaClusterID)
clustering.table
Memory CD4 T CD14+ Mono Naive CD4 T NK B CD8 T CD16+ Mono T/Mono doublets CD34+ Eryth Mk DC
0 1754 0 1217 29 0 27 0 5 2 4 24 1
1 0 2189 0 4 0 0 30 1 1 5 25 55
2 3 0 2 890 3 1 0 0 1 3 7 2
3 0 4 0 2 319 0 2 0 2 2 3 0
4 24 0 18 4 1 243 0 0 0 1 2 0
5 1 27 4 157 2 2 10 56 0 9 16 6
6 4 5 0 1 0 0 0 1 113 81 16 5
7 4 59 4 0 0 0 9 117 0 0 2 0
8 0 9 0 2 0 0 179 0 0 0 1 0
9 0 0 1 0 0 0 0 0 0 0 0 1
10 1 0 2 0 25 0 0 2 0 0 0 0
pDCs
0 2
1 0
2 1
3 0
4 0
5 2
6 0
7 1
8 0
9 43
10 0
根據(jù)上面的結(jié)果對(duì)蛋白質(zhì)水平的細(xì)胞進(jìn)行細(xì)胞類(lèi)型鑒定
new.cluster.ids <- c("CD4 T", "CD14+ Mono", "NK", "B", "CD8 T", "NK", "CD34+", "T/Mono doublets",
"CD16+ Mono", "pDCs", "B")
names(new.cluster.ids) <- levels(cbmc)
cbmc <- RenameIdents(cbmc, new.cluster.ids)
tsne_rnaClusters <- DimPlot(cbmc, reduction = "tsne_adt", group.by = "rnaClusterID") + NoLegend()
tsne_rnaClusters <- tsne_rnaClusters + ggtitle("Clustering based on scRNA-seq") + theme(plot.title = element_text(hjust = 0.5))
tsne_rnaClusters <- LabelClusters(plot = tsne_rnaClusters, id = "rnaClusterID", size = 4)
tsne_adtClusters <- DimPlot(cbmc, reduction = "tsne_adt", pt.size = 0.5) + NoLegend()
tsne_adtClusters <- tsne_adtClusters + ggtitle("Clustering based on ADT signal") + theme(plot.title = element_text(hjust = 0.5))
tsne_adtClusters <- LabelClusters(plot = tsne_adtClusters, id = "ident", size = 4)
注意:為了進(jìn)行比較，RNA和蛋白的聚類(lèi)都是在使用ADT距離矩陣生成的tSNE上顯示的。
wrap_plots(list(tsne_rnaClusters, tsne_adtClusters), ncol = 2)
基于adt的聚類(lèi)產(chǎn)生了類(lèi)似的結(jié)果，但有一些區(qū)別：
- 基于CD4、CD8、CD14和CD45RA的強(qiáng)大ADT數(shù)據(jù)，CD4/CD8 T細(xì)胞群的聚類(lèi)得到了改善
- 然而，一些ADT數(shù)據(jù)中不包含有很好的區(qū)別蛋白標(biāo)記的簇(如Mk/Ery/DC)會(huì)失去分離
- 也可以在RNA水平上使用findmarker驗(yàn)證這一點(diǎn)
tcells <- subset(cbmc, idents = c("CD4 T", "CD8 T"))
FeatureScatter(tcells, feature1 = "CD4", feature2 = "CD8")
RidgePlot(cbmc, features = c("adt_CD11c", "adt_CD8", "adt_CD16", "adt_CD4", "adt_CD19", "adt_CD14"),
ncol = 2)
保存結(jié)果
saveRDS(cbmc, file = "Analysis/cbmc_multimodal.rds")
