第8章: 使用scRNA-seq定义cluster类型
除了使用基因得分定义细胞类群以外,ArchR还能整合scRNA-seq数据。通过将scATAC-seq数据里的基因得分矩阵和scRNA-seq数据的基因表达量矩阵进行对比,ArchR就能将scATAC-seq的细胞比对到scRNA-seq的细胞,实现两种数据的整合。之后,我们借助scRNA-seq数据已经定义的细胞类群,或者整合后的scRNA-seq的基因表达量来注释细胞类群。在代码内部,我们调用了Seurat::FindTransferAnchors()
,从而实现两种数据集之间的比较。当然,ArchR不是简单地对函数进行封装,而是在此基础上通过对数据的拆分,实现并行计算,从而能将该流程应用到更大规模的细胞中。
该整合过程实际上会寻找scATAC-seq和scRNA-seq两者中最相似的细胞,然后将scRNA-seq中对应的细胞表达量赋值给scATAAC-seq细胞。最后,scATAC-seq中每个细胞都会有基因表达量特征,能被用于下游分析。这一章会阐述如何利用该信息定义细胞类型,之后会介绍如何使用连接的scRNA-seq数据做更加复杂的分析,例如识别预测的顺式调控元件。我们相信由于多组学单细胞谱的商业化,这类整合分析将会越来越多。同时,在ArchR中使用公共数据里匹配细胞类型的scRNA-seq数据或者自己使用目标样本得到的scRNA-seq数据也能加强scATAC-seq分析。
8.1 scATAC-seq细胞和scRNA-seq细胞跨平台连接
为了能将我们教程中的scATAC-seq数据与其匹配的scRNA-seq数据进行整合,我们将使用 Granja* et al (2019) 里的造血细胞scRNA-seq数据。
scRNA-seq数据以 RangedSummarizedExperiment
对象保存,大小为111MB。此外,ArchR还接受未经修改的Seurat
对象作为整合流程的输入。我们使用download.file
下载数据
1 | if(!file.exists("scRNA-Hematopoiesis-Granja-2019.rds")){ |
下载之后,我们使用readRDS
进行读取,并查看该对象
1 | seRNA <- readRDS("scRNA-Hematopoiesis-Granja-2019.rds") |
从输出信息中,我们可以发现它里面有基因表达量的count矩阵和对应的元信息。
元信息列里的BioClassification
记录着scRNA-seq数据中每个细胞对应的细胞类型分类
1 | colnames(colData(seRNA)) |
使用table()
,我们可以看到scRNA-seq细胞类型每一群的细胞数
1 | table(colData(seRNA)$BioClassification) |
我们后续会用到两种整合方法。第一种是无约束整合,我们对scATAC-seq实验里的细胞不作任何假设,直接尝试将这些细胞和scRNA-seq实验里的任意细胞进行配对。这是一种初步可行方案,后续会根据这一步得到的结果,对整合步骤进行约束来提升跨平台配对的质量。第二种方法是约束整合,即利用对细胞类型的先验知识限制搜索范围。举个例子,如果我们知道scATAC-seq中的Cluster A,B,C对应着三种不同的T细胞,scRNA-seq中的Cluster X,Y对应着两种不同的T细胞,我们告诉ArchR只需要尝试将scATAC-seq中的Cluster A,B,C跟scRNA-seq中的Cluster X,Y进行配对。下面,我们将先以无约束整合初步地鉴定每一种聚类的类型,然后根据分析结果做更加细致的约束整合。
8.1.1 无约束整合
我们使用addGeneIntegrationMatrix()
对scATAC-seq和scRNA-seq数据进行整合。正如之前所提到的,该函数的seRNA
参数接受Seurat
或RangedSummarizedExperiment
对象作为输入。因为第一轮是探索性质的无约束整合,因此,我们不会将结果保存在Arrow文件中(addToArrow=FALSE
)。整合后的矩阵将会根据matrixName
进行命名,存放在ArchRProject
中。该函数的其他参数对应cellColData
中列名用于存放额外的信息,nameCell
存放scRNA-seq中匹配的细胞ID,nameGroup
存放scRNA-seq细胞中的分组ID,nameScore
存放跨平台整合得分。
1 | projHeme2 <- addGeneIntegrationMatrix( |
无约束整合的结果可能并不准确,但是为后续更加精细的约束分析奠定了基础。
8.1.2 约束整合
现在我们已经有了初步的无约束整合结果,我们就有了大致的细胞类型分布情况,接着就是优化整合结果。
因为我们教程里的数据来自于造血细胞,我们将会在理想状态下将类似的细胞整合在一起。首先,我们先确认scRNA-seq里细胞类型在我们的scATAC-seq聚类中分布情况。这一步的目标是使用无约束整合的方法找到scATAC-seq和scRNA-seq中和T细胞和NK细胞对应的聚类,后续会用到该信息进行约束整合。具体操作为,我们创建一个confusionMatrix
,并关注Cluster
和predictedGroup_Un
的交叉部分中scRNA-seq的细胞类型。
1 | cM <- as.matrix(confusionMatrix(projHeme2$Clusters, projHeme2$predictedGroup_Un)) |
输出信息如下,展示了12个scATAC-seq聚类中对应最优可能的scRNA-seq细胞类型。
1 | # preClust |
首先,我们检查在无约束整合中用到的scRNA-seq数据里细胞类型标签。
1 | unique(unique(projHeme2$predictedGroup_Un)) |
从上面的输出中,我们发现scRNA-seq数据中与NK细胞和T细胞对应的聚类是Cluster 19 - 25。
接着我们创建一个字符串模式用来表示这些聚类,后续的约束整合会用到,其中|
在正则表达式中表示”或”,我们之后使用grep
根据这些字符串模式从scATAC-seq提取和scRNA-seq对应的聚类。。
1 | #From scRNA |
其余的聚类就称之为”Non-T cell, Non-NK cell”(例如Cluster 1 - 18),也创建了对应的字符串模式
1 | cNonTNK <- paste0(c(paste0("0", 1:9), 10:13, 15:18), collapse="|") |
接着再用字符串模式在preClust
找到对应的scATAC-seq列名,然后使用列名从混合矩阵提取对应的列。
对于T细胞和NK细胞,scATAC-seq聚类ID就是C7, C8, C9
1 | #Assign scATAC to these categories |
对于” Non-T cells and Non-NK cells”, ID就是scATAC-seq聚类余下的部分
1 | clustNonTNK <- rownames(cM)[grep(cNonTNK, preClust)] |
接着在scRNA-seq中做相同的操作,筛选出相同的细胞类型。首先,我们鉴定scRNA-seq数据中T细胞和NK细胞
1 | #RNA get cells in these categories |
然后,鉴定scRNA-seq数据中”Non-T cell Non-NK cell cells”
1 | rnaNonTNK <- colnames(seRNA)[grep(cNonTNK, colData(seRNA)$BioClassification)] |
约束整合需要我们提供一个嵌套list。这是一个由多个SimpleList
对象组成的SimpleList
, 每一组对应一个约束情况。在案例中,我们有两个组,一个组称之为TNK
,包括两个平台的T和NK细胞,另一个组为NonTNK
,包括两个平台的”Non-T cell Non-NK cell”细胞。每个SimpleList
对象都包含两个细胞ID的向量,一个是ATAC,一个是RNA.
1 | groupList <- SimpleList( |
我们将该列表传递给addGeneIntegrationMatrix()
函数的groupList
参数。注意,我们依旧没有将结果添加到Arrow文件中 (addToArrow = FALSE
)。我们强烈建议,在保存到Arrow文件前先彻底的检查结果,看结果是否符合预期。在教程的下一节会介绍该操作。
1 | projHeme2 <- addGeneIntegrationMatrix( |
8.1.3 约束整合和无约束整合对比
正如之前所提到的,我们的scATAC-seq数据可以根据整合的scRNA-seq数据进行定义,并且有约束和无约束这两种方式。为了对两者进行对比,我们分别根据这两种整合结果对scATAC-seq的数据进行上色。
首先,使用ArchR内置的paletteDiscrete()
函数生成调色板
1 | pal <- paletteDiscrete(values = colData(seRNA)$BioClassification) |
在ArchR中,调色板本质上一个命名向量,每个十六进制编码的颜色对应着一个名字。
1 | pal |
我们在scATAC-seq数据根据无约束整合得到的scRNA-seq细胞类型进行可视化
1 | p1 <- plotEmbedding( |
1 | do.call(cowplot::plot_grid, c(list(ncol = 3), p2c)) |
同样,我们也可以根据约束整合得到scATAC-seq对应的scRNA-seq的细胞类型进行可视化
1 | p2 <- plotEmbedding( |
这两者的结果差异其实非常细微,主要是我们感兴趣的细胞类型原本就存在明显的差异。当然,仔细观察还能发现其中的不同之处,尤其是T细胞(Clusters 17-22)
我们用plotPDF()
函数保存该图矢量版本。
1 | plotPDF(p1,p2, name = "Plot-UMAP-RNA-Integration.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5) |
我们现在可以用saveArchRProject()
函数保存我们的projHeme2
对象。
1 | saveArchRProject(ArchRProj = projHeme2, outputDirectory = "Save-ProjHeme2", load = FALSE) |
8.2 为每个scATAC-seq细胞增加拟scRNA-seq谱
既然我们对scATAC-seq和scRNA-seq整合的结果感到满意,我们就能用addToArrow=TRUE
重新运行,在Arrow文件中添加相关联的表达量矩阵数据。根据之前所提到的,我们传入groupList
约束整合,在nameCell
,nameGroup
和nameScore
中加入列名。这些元信息列都会被添加到cellColData
中。
这里,我们新建了一个projHeme3
,用于后续教程。
1 | #~5 minutes |
现在,当我们使用getAvailableMatrices()
检查哪些矩阵可用时,我们会发现GeneIntegrationMatrix
已经被添加到Arrow文件中
1 | getAvailableMatrices(projHeme3) |
在新的GeneIntegrationMatrix
中,我们可以比较连接的基因表达量和根据基因得分推断的基因表达量
我们需要先确保在项目中加入了填充权重值(impute weights)
1 | projHeme3 <- addImputeWeights(projHeme3) |
现在,我们来生成一些UMAP图,里面的基因表达量值来自于GeneIntegrationMatrix
1 | markerGenes <- c( |
以及一些相同UMAP图,但是使用GeneScoreMatrix
里的基因得分值
1 | p2 <- plotEmbedding( |
最后用cowplot
将这些标记基因绘制在一起
1 | p1c <- lapply(p1, function(x){ |
1 | do.call(cowplot::plot_grid, c(list(ncol = 3), p2c)) |
和预期的一样,两个方法推测的基因表达量存在相似性,但并不相同。
使用plotPDF()
函数保存可编辑的矢量版。
1 | plotPDF(plotList = p1, |
8.3 使用scRNA-seq信息标记scATAC-seq聚类
现在,我们确定了scATAC-seq和scRNA-seq数据间的对应关系,我们就可以使用scRNA-seq数据中细胞类型对我们的scATAC-seq聚类进行定义。
首先,我们会在scATAC-seq和整合分析得到predictedGroup
之间构建一个混合矩阵
1 | cM <- confusionMatrix(projHeme3$Clusters, projHeme3$predictedGroup) |
接着,对于每一个scATAC-seq聚类,我们根据predictedGroup
确定最能定义聚类的细胞类型。
1 | labelNew <- colnames(cM)[apply(cM, 1, which.max)] |
接着,我们需要对新的聚类标签进行重命名,简化分类系统。对于每一个scRNA-seq的聚类,我们会重新定义标签,以便更好解释。
1 | remapClust <- c( |
接着,使用mapLables()
函数进行标签转换,将旧的标签映射到新的标签上。
1 | labelNew2 <- mapLabels(labelNew, oldLabels = names(remapClust), newLabels = remapClust) |
合并labelsOld
和labelsNew2
,我们现在可以用mapLables()
函数在cellColData
里新建聚类标签。
1 | projHeme3$Clusters2 <- mapLabels(projHeme3$Clusters, newLabels = labelNew2, oldLabels = labelOld) |
得到新的标签后,我们绘制最新的UMAP展示聚类结果。
1 | p1 <- plotEmbedding(projHeme3, colorBy = "cellColData", name = "Clusters2") |
如果被分析scATAC-seq数据对应的细胞系统也有scRNA-seq数据存在,那么这种分析范式将会特别有用。正如之前所说,这种scRNA-seq和scATAC-seq整合分析也为后续更加复杂的基因调控分析提供了出色的框架。
用plotPDF()
函数保存该图矢量版本。
1 | plotPDF(p1, name = "Plot-UMAP-Remap-Clusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5) |
使用saveArchRProject
保存我们最初的projHeme3.
1 | saveArchRProject(ArchRProj = projHeme3, outputDirectory = "Save-ProjHeme3", load = FALSE) |