教程背景
这个例子演示了允许用户使用Seurat分析和探索多模式数据的新特性。虽然这只是一个初始版本，但我们很高兴在未来为多模态数据集发布重要的新功能。
在这里，我们分析了由CITE-seq产生的8,617个脐带血单核细胞(CBMCs)的数据集，我们同时测量了单细胞转录组和11个表面蛋白的表达，其水平通过DNA-barcode抗体来量化。首先，我们加载两个计数矩阵:一个用于RNA测量，另一个用于抗体衍生标签ADT。你可以在这里下载ADT文件和RNA文件
1.加载RNA UMI矩阵
注意，这个数据集还包含了大约5%的小鼠细胞，我们可以将其作为蛋白质测量的阴性对照。因此，基因表达矩阵在每个基因的开头都附加了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))
为了让生活更容易继续下去，我们将抛弃除前100个高度表达的老鼠基因外的所有基因，并从CITE-seq前缀中去除“HUMAN_”
# 这个函数有几个默认参数，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添加多模态数据时，可以使用重复的特性名称。每一组模态数据(例如。RNA、ADT等)储存在自己的Assay对象中。其中一个Assay对象被称为“default assay”，意思是它用于所有的分析和可视化。若要从非默认Assay中提取数据，可以指定一个与Assay链接的键，用于提取特征。要查看所有对象的所有键，请使用Key function。
最后，我们观察到CCR5、CCR7和CD10的低富集，因此将它们从矩阵中去除(可选)。(关键是这三个蛋白得低富集怎么看出来的？)
cbmc.adt <- cbmc.adt[setdiff(rownames(x = cbmc.adt), c("CCR5", "CCR7", "CD10")), ]
3.设置一个Seurat对象，并基于RNA表达对细胞进行聚类
下面的步骤表示基于scRNA-seq数据对pbmc进行快速聚类。有关单个步骤或更高级选项的详细信息，请参见这里的PBMC集群指导教程。
cbmc <- CreateSeuratObject(counts = cbmc.rna) # 创建seurat对象
cbmc <- NormalizeData(cbmc) # 标准化
cbmc <- FindVariableFeatures(cbmc) # 选择高变基因
cbmc <- ScaleData(cbmc) # 归一化
cbmc <- RunPCA(cbmc, verbose = FALSE) # 运行PCA降维
ElbowPlot(cbmc, ndims = 50) # 确定取多少个PCs
这里教程说选择13个PCs，但是看代码选的是25个，蜜汁迷惑。
cbmc <- FindNeighbors(cbmc, dims = 1:25)
cbmc <- FindClusters(cbmc, resolution = 0.8)
cbmc <- RunTSNE(cbmc, dims = 1:25, method = "FIt-SNE")
# 找到定义每个cluster的标记，并使用它们对cluster进行注释，我们使用max.cells.per加快进程
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
注意，为了简单起见，我们将两个CD14+单核细胞簇(HLA-DR基因表达不同)和NK细胞簇(细胞周期阶段不同)合并。
这里直接进行了细胞类型注释，注释得过程？
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.将蛋白质表达水平添加到Seurat对象
Seurat v3.0允许您将来自多个分析的信息存储在同一个对象中，只要数据是多模态的(在同一组细胞中收集)。您可以使用SetAssayData和GetAssayData访问器函数来添加和从其他分析中获取数据。
我们将定义一个ADT试验，并为它存储原始计数。
如果你对这些数据内部如何存储感兴趣，您可以查看Assay类，它在object . R中定义。注意，所有的单细胞表达数据，包括RNA数据，仍然存储在Assay对象中，也可以使用GetAssayData访问。
cbmc[["ADT"]] <- CreateAssayObject(counts = cbmc.adt)
现在我们可以重复我们通常在RNA上运行的预处理(标准化和缩放)步骤，但是需要修改‘assay’的参数。对于CITE-seq数据，我们不建议使用典型的对数标准化。相反，我们使用中心对数比(CLR)归一化，独立计算每个特性。与原始版本相比，这是一个略微改进的过程，我们将很快发布更高级的CITE-seq规范化版本。
# 这里的标准化方法是CLR
cbmc <- NormalizeData(cbmc, assay = "ADT", normalization.method = "CLR")
cbmc <- ScaleData(cbmc, assay = "ADT")
5.观察RNA簇上的蛋白水平
您可以使用任何ADT标记的名称。“adt_CD4”)，在FetchData、FeaturePlot、RidgePlot、FeatureScatter、DoHeatmap或任何其他可视化特性中。
Q1：这里得ADTfeatures为什么突然多了三个字母adt_？
Q2：这里得蛋白和RNA对应得表达水平，从下图可以看出来什么？
# 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)
在这个图中，蛋白质(ADT)水平在上面，RNA水平在下面。
山脊图绘制
RidgePlot(cbmc, features = c("adt_CD3", "adt_CD11c", "adt_CD8", "adt_CD16"), ncol = 2)
绘制ADT散点图(如FACS的双轴图)。注意，如果需要的话，你甚至可以使用HoverLocator和FeatureLocator来“gate”细胞。
FeatureScatter(cbmc, feature1 = "adt_CD19", feature2 = "adt_CD3")
查看蛋白质和RNA之间的关系
FeatureScatter(cbmc, feature1 = "adt_CD3", feature2 = "CD3E")
查看T细胞中CD4和CD8的水平
tcells <- subset(cbmc, idents = c("Naive CD4 T", "Memory CD4 T", "CD8 T"))
FeatureScatter(tcells, feature1 = "adt_CD4", feature2 = "adt_CD8")
让我们看看原始的(非标准化的)ADT的count值。你可以看到这些值相当高，特别是与RNA值相比。这是由于细胞中蛋白质拷贝数显著增加，这大大减少了ADT数据中的“drop-out”。
对于drop-out如何理解？
FeatureScatter(tcells, feature1 = "adt_CD4", feature2 = "adt_CD8", slot = "counts")
如果你仔细观察，你会发现我们的CD8 T细胞群富含CD8 T细胞，但仍然包含许多CD4+ CD8- T细胞。这是因为初始CD4和CD8 T细胞在转录上非常相似，CD4和CD8的RNA丢失水平相当高。这说明了单独从scRNA-seq数据中定义细微的免疫细胞差异的挑战。
saveRDS(cbmc, file = "Analysis/cbmc.rds")
6.识别簇间差异表达的蛋白
对cluster进行采样，每个cluster最多300个细胞(使小cluster的热图更容易看到)。
cbmc.small <- subset(cbmc, downsample = 300)
找到所有cluster的蛋白质标记，并绘制一个热图。
adt.markers <- FindAllMarkers(cbmc.small, assay = "ADT", only.pos = TRUE)
DoHeatmap(cbmc.small, features = unique(adt.markers$gene), assay = "ADT", angle = 90) + NoLegend()
你可以看到我们的未知细胞同时表达骨髓和淋巴标记(在RNA水平上也是如此)。它们可能是应该被丢弃的细胞团块(多胞胎)。我们现在也要移除老鼠细胞。
cbmc <- subset(cbmc, idents = c("Multiplets", "Mouse"), invert = TRUE)
7.直接在蛋白质水平上聚类
还可以直接在CITE-seq数据上运行降维和基于图的聚类。
#因为我们将广泛地使用ADT数据，所以我们将把默认assay改为“ADT”assay。这将导致所有函数默认使用ADT数据，而不是要求我们每次都指定它
DefaultAssay(cbmc) <- "ADT"
cbmc <- RunPCA(cbmc, features = rownames(cbmc), reduction.name = "pca_adt", reduction.key = "pca_adt_",
verbose = FALSE)
DimPlot(cbmc, reduction = "pca_adt")
在我们在ADT级别上重新聚类数据之前，我们将稍后隐藏RNA集群id。
cbmc[["rnaClusterID"]] <- Idents(cbmc)
现在，我们仅在ADT(蛋白质)水平上使用PCA重新计算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和蛋白质的聚类，并使用它来注释蛋白质聚类。
(当然，我们也可以使用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
根据上面的结果对蛋白质水平的细胞进行细胞类型鉴定
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)
注意:为了进行比较，RNA和蛋白的聚类都是在使用ADT距离矩阵生成的tSNE上显示的。
wrap_plots(list(tsne_rnaClusters, tsne_adtClusters), ncol = 2)
基于adt的聚类产生了类似的结果，但有一些区别：
- 基于CD4、CD8、CD14和CD45RA的强大ADT数据，CD4/CD8 T细胞群的聚类得到了改善
- 然而，一些ADT数据中不包含有很好的区别蛋白标记的簇(如Mk/Ery/DC)会失去分离
- 也可以在RNA水平上使用findmarker验证这一点
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)
保存结果
saveRDS(cbmc, file = "Analysis/cbmc_multimodal.rds")
