使用ArchR分析单细胞ATAC-seq数据(第八章)

第8章: 使用scRNA-seq定义cluster类型

除了使用基因得分定义细胞类群以外,ArchR还能整合scRNA-seq数据。通过将scATAC-seq数据里的基因得分矩阵和scRNA-seq数据的基因表达量矩阵进行对比,ArchR就能将scATAC-seq的细胞比对到scRNA-seq的细胞,实现两种数据的整合。之后,我们借助scRNA-seq数据已经定义的细胞类群,或者整合后的scRNA-seq的基因表达量来注释细胞类群。在代码内部,我们调用了Seurat::FindTransferAnchors(),从而实现两种数据集之间的比较。当然,ArchR不是简单地对函数进行封装,而是在此基础上通过对数据的拆分,实现并行计算,从而能将该流程应用到更大规模的细胞中。

img

该整合过程实际上会寻找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下载数据

if(!file.exists("scRNA-Hematopoiesis-Granja-2019.rds")){
    download.file(
        url = "https://jeffgranja.s3.amazonaws.com/ArchR/TestData/scRNA-Hematopoiesis-Granja-2019.rds",
        destfile = "scRNA-Hematopoiesis-Granja-2019.rds"
    )
}

下载之后,我们使用readRDS进行读取,并查看该对象

seRNA <- readRDS("scRNA-Hematopoiesis-Granja-2019.rds")
seRNA
## class: RangedSummarizedExperiment
## dim: 20287 35582
## metadata(0):
## assays(1): counts
## rownames(20287): FAM138A OR4F5 … S100B PRMT2
## rowData names(3): gene_name gene_id exonLength
## colnames(35582): CD34_32_R5:AAACCTGAGTATCGAA-1
## CD34_32_R5:AAACCTGAGTCGTTTG-1 …
## BMMC_10x_GREENLEAF_REP2:TTTGTTGCATGTGTCA-1
## BMMC_10x_GREENLEAF_REP2:TTTGTTGCATTGAAAG-1
## colData names(10): Group nUMI_pre … BioClassification Barcode

从输出信息中,我们可以发现它里面有基因表达量的count矩阵和对应的元信息。

元信息列里的BioClassification记录着scRNA-seq数据中每个细胞对应的细胞类型分类

colnames(colData(seRNA))
# [1] "Group"             "nUMI_pre"          "nUMI"             
# [4] "nGene"             "initialClusters"   "UMAP1"            
# [7] "UMAP2"             "Clusters"          "BioClassification"
# [10] "Barcode"

使用table(),我们可以看到scRNA-seq细胞类型每一群的细胞数

table(colData(seRNA)$BioClassification)
#        01_HSC 02_Early.Eryth  03_Late.Eryth  04_Early.Baso    05_CMP.LMPP 
#          1425           1653            446            111           2260 
#      06_CLP.1         07_GMP    08_GMP.Neut         09_pDC         10_cDC 
#           903           2097           1050            544            325 
#11_CD14.Mono.1 12_CD14.Mono.2   13_CD16.Mono         14_Unk       15_CLP.2 
#          1800           4222            292            520            377 
#      16_Pre.B           17_B      18_Plasma       19_CD8.N      20_CD4.N1 
#           710           1711             62           1521           2470 
#     21_CD4.N2       22_CD4.M      23_CD8.EM      24_CD8.CM          25_NK 
#          2364           3539            796           2080           2143 
#        26_Unk 
#           161 

我们后续会用到两种整合方法。第一种是无约束整合,我们对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参数接受SeuratRangedSummarizedExperiment对象作为输入。因为第一轮是探索性质的无约束整合,因此,我们不会将结果保存在Arrow文件中(addToArrow=FALSE)。整合后的矩阵将会根据matrixName进行命名,存放在ArchRProject中。该函数的其他参数对应cellColData中列名用于存放额外的信息,nameCell存放scRNA-seq中匹配的细胞ID,nameGroup存放scRNA-seq细胞中的分组ID,nameScore存放跨平台整合得分。

projHeme2 <- addGeneIntegrationMatrix(
    ArchRProj = projHeme2, 
    useMatrix = "GeneScoreMatrix",
    matrixName = "GeneIntegrationMatrix",
    reducedDims = "IterativeLSI",
    seRNA = seRNA,
    addToArrow = FALSE,
    groupRNA = "BioClassification",
    nameCell = "predictedCell_Un",
    nameGroup = "predictedGroup_Un",
    nameScore = "predictedScore_Un"
)

无约束整合的结果可能并不准确,但是为后续更加精细的约束分析奠定了基础。

8.1.2 约束整合

现在我们已经有了初步的无约束整合结果,我们就有了大致的细胞类型分布情况,接着就是优化整合结果。

因为我们教程里的数据来自于造血细胞,我们将会在理想状态下将类似的细胞整合在一起。首先,我们先确认scRNA-seq里细胞类型在我们的scATAC-seq聚类中分布情况。这一步的目标是使用无约束整合的方法找到scATAC-seq和scRNA-seq中和T细胞和NK细胞对应的聚类,后续会用到该信息进行约束整合。具体操作为,我们创建一个confusionMatrix,并关注ClusterpredictedGroup_Un的交叉部分中scRNA-seq的细胞类型。

cM <- as.matrix(confusionMatrix(projHeme2$Clusters, projHeme2$predictedGroup_Un))
preClust <- colnames(cM)[apply(cM, 1 , which.max)]
cbind(preClust, rownames(cM)) #Assignments

输出信息如下,展示了12个scATAC-seq聚类中对应最优可能的scRNA-seq细胞类型。

#      preClust              
# [1,] "03_Late.Eryth"  "C10"
# [2,] "20_CD4.N1"      "C8" 
# [3,] "16_Pre.B"       "C3" 
# [4,] "08_GMP.Neut"    "C11"
# [5,] "17_B"           "C4" 
# [6,] "11_CD14.Mono.1" "C1" 
# [7,] "01_HSC"         "C12"
# [8,] "22_CD4.M"       "C9" 
# [9,] "09_pDC"         "C5" 
# [10,] "25_NK"          "C7" 
# [11,] "12_CD14.Mono.2" "C2" 
# [12,] "06_CLP.1"       "C6"

首先,我们检查在无约束整合中用到的scRNA-seq数据里细胞类型标签。

unique(unique(projHeme2$predictedGroup_Un))
# [1] "08_GMP.Neut"    "25_NK"          "16_Pre.B"       "06_CLP.1"      
# [5] "07_GMP"         "11_CD14.Mono.1" "04_Early.Baso"  "22_CD4.M"      
# [9] "03_Late.Eryth"  "05_CMP.LMPP"    "17_B"           "19_CD8.N"      
#[13] "09_pDC"         "13_CD16.Mono"   "23_CD8.EM"      "12_CD14.Mono.2"
#[17] "20_CD4.N1"      "02_Early.Eryth" "21_CD4.N2"      "24_CD8.CM"     
#[21] "01_HSC" 

从上面的输出中,我们发现scRNA-seq数据中与NK细胞和T细胞对应的聚类是Cluster 19 - 25。

接着我们创建一个字符串模式用来表示这些聚类,后续的约束整合会用到,其中|在正则表达式中表示"或",我们之后使用grep根据这些字符串模式从scATAC-seq提取和scRNA-seq对应的聚类。。

#From scRNA
cTNK <- paste0(paste0(19:25), collapse="|")
cTNK
# [1] "19|20|21|22|23|24|25"

其余的聚类就称之为"Non-T cell, Non-NK cell"(例如Cluster 1 - 18),也创建了对应的字符串模式

cNonTNK <- paste0(c(paste0("0", 1:9), 10:13, 15:18), collapse="|")
cNonTNK
#[1] "01|02|03|04|05|06|07|08|09|10|11|12|13|15|16|17|18"

接着再用字符串模式在preClust找到对应的scATAC-seq列名,然后使用列名从混合矩阵提取对应的列。

对于T细胞和NK细胞,scATAC-seq聚类ID就是C7, C8, C9

#Assign scATAC to these categories
clustTNK <- rownames(cM)[grep(cTNK, preClust)]
clustTNK
#[1] "C8" "C9" "C7"

对于" Non-T cells and Non-NK cells", ID就是scATAC-seq聚类余下的部分

clustNonTNK <- rownames(cM)[grep(cNonTNK, preClust)]
clustNonTNK
# [1] "C10" "C3"  "C11" "C4"  "C1"  "C12" "C5"  "C2"  "C6" 

接着在scRNA-seq中做相同的操作,筛选出相同的细胞类型。首先,我们鉴定scRNA-seq数据中T细胞和NK细胞

#RNA get cells in these categories
rnaTNK <- colnames(seRNA)[grep(cTNK, colData(seRNA)$BioClassification)]
head(rnaTNK)
#[1] "PBMC_10x_GREENLEAF_REP1:AAACCCAGTCGTCATA-1"
#[2] "PBMC_10x_GREENLEAF_REP1:AAACCCATCCGATGTA-1"
#[3] "PBMC_10x_GREENLEAF_REP1:AAACCCATCTCAACGA-1"
#[4] "PBMC_10x_GREENLEAF_REP1:AAACCCATCTCTCGAC-1"
#[5] "PBMC_10x_GREENLEAF_REP1:AAACGAACAATCGTCA-1"
#[6] "PBMC_10x_GREENLEAF_REP1:AAACGAACACGATTCA-1"

然后,鉴定scRNA-seq数据中"Non-T cell Non-NK cell cells"

rnaNonTNK <- colnames(seRNA)[grep(cNonTNK, colData(seRNA)$BioClassification)]
head(rnaNonTNK)
#[1] "CD34_32_R5:AAACCTGAGTATCGAA-1" "CD34_32_R5:AAACCTGAGTCGTTTG-1"
#[3] "CD34_32_R5:AAACCTGGTTCCACAA-1" "CD34_32_R5:AAACGGGAGCTTCGCG-1"
#[5] "CD34_32_R5:AAACGGGAGGGAGTAA-1" "CD34_32_R5:AAACGGGAGTTACGGG-1"

约束整合需要我们提供一个嵌套list。这是一个由多个SimpleList对象组成的SimpleList, 每一组对应一个约束情况。在案例中,我们有两个组,一个组称之为TNK,包括两个平台的T和NK细胞,另一个组为NonTNK,包括两个平台的"Non-T cell Non-NK cell"细胞。每个SimpleList对象都包含两个细胞ID的向量,一个是ATAC,一个是RNA.

groupList <- SimpleList(
    TNK = SimpleList(
        ATAC = projHeme2$cellNames[projHeme2$Clusters %in% clustTNK],
        RNA = rnaTNK
    ),
    NonTNK = SimpleList(
        ATAC = projHeme2$cellNames[projHeme2$Clusters %in% clustNonTNK],
        RNA = rnaNonTNK
    )    
)

我们将该列表传递给addGeneIntegrationMatrix()函数的groupList参数。注意,我们依旧没有将结果添加到Arrow文件中 (addToArrow = FALSE)。我们强烈建议,在保存到Arrow文件前先彻底的检查结果,看结果是否符合预期。在教程的下一节会介绍该操作。

projHeme2 <- addGeneIntegrationMatrix(
    ArchRProj = projHeme2, 
    useMatrix = "GeneScoreMatrix",
    matrixName = "GeneIntegrationMatrix",
    reducedDims = "IterativeLSI",
    seRNA = seRNA,
    addToArrow = FALSE, 
    groupList = groupList,
    groupRNA = "BioClassification",
    nameCell = "predictedCell_Co",
    nameGroup = "predictedGroup_Co",
    nameScore = "predictedScore_Co"
)

8.1.3 约束整合和无约束整合对比

正如之前所提到的,我们的scATAC-seq数据可以根据整合的scRNA-seq数据进行定义,并且有约束和无约束这两种方式。为了对两者进行对比,我们分别根据这两种整合结果对scATAC-seq的数据进行上色。

首先,使用ArchR内置的paletteDiscrete()函数生成调色板

pal <- paletteDiscrete(values = colData(seRNA)$BioClassification)

在ArchR中,调色板本质上一个命名向量,每个十六进制编码的颜色对应着一个名字。

pal
#        01_HSC 02_Early.Eryth  03_Late.Eryth  04_Early.Baso    05_CMP.LMPP 
#     "#D51F26"      "#502A59"      "#235D55"      "#3D6E57"      "#8D2B8B" 
#      06_CLP.1         07_GMP    08_GMP.Neut         09_pDC         10_cDC 
#     "#DE6C3E"      "#F9B712"      "#D8CE42"      "#8E9ACD"      "#B774B1" 
#11_CD14.Mono.1 12_CD14.Mono.2   13_CD16.Mono         14_Unk       15_CLP.2 
#     "#D69FC8"      "#C7C8DE"      "#8FD3D4"      "#89C86E"      "#CC9672" 
#      16_Pre.B           17_B      18_Plasma       19_CD8.N      20_CD4.N1 
#     "#CF7E96"      "#A27AA4"      "#CD4F32"      "#6B977E"      "#518AA3" 
#     21_CD4.N2       22_CD4.M      23_CD8.EM      24_CD8.CM          25_NK 
#     "#5A5297"      "#0F707D"      "#5E2E32"      "#A95A3C"      "#B28D5C" 
#        26_Unk 
#     "#3D3D3D"

我们在scATAC-seq数据根据无约束整合得到的scRNA-seq细胞类型进行可视化

p1 <- plotEmbedding(
    projHeme2, 
    colorBy = "cellColData", 
    name = "predictedGroup_Un", 
    pal = pal
)
p1
Plot-UMAP-RNA-Integration_1
do.call(cowplot::plot_grid, c(list(ncol = 3), p2c))

同样,我们也可以根据约束整合得到scATAC-seq对应的scRNA-seq的细胞类型进行可视化

p2 <- plotEmbedding(
    projHeme2, 
    colorBy = "cellColData", 
    name = "predictedGroup_Co", 
    pal = pal
)
p2
Plot-UMAP-RNA-Integration_2

这两者的结果差异其实非常细微,主要是我们感兴趣的细胞类型原本就存在明显的差异。当然,仔细观察还能发现其中的不同之处,尤其是T细胞(Clusters 17-22)

我们用plotPDF()函数保存该图矢量版本。

plotPDF(p1,p2, name = "Plot-UMAP-RNA-Integration.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)

我们现在可以用saveArchRProject()函数保存我们的projHeme2对象。

saveArchRProject(ArchRProj = projHeme2, outputDirectory = "Save-ProjHeme2", load = FALSE)

8.2 为每个scATAC-seq细胞增加拟scRNA-seq谱

既然我们对scATAC-seq和scRNA-seq整合的结果感到满意,我们就能用addToArrow=TRUE重新运行,在Arrow文件中添加相关联的表达量矩阵数据。根据之前所提到的,我们传入groupList约束整合,在nameCellnameGroupnameScore中加入列名。这些元信息列都会被添加到cellColData中。

这里,我们新建了一个projHeme3,用于后续教程。

#~5 minutes
projHeme3 <- addGeneIntegrationMatrix(
    ArchRProj = projHeme2, 
    useMatrix = "GeneScoreMatrix",
    matrixName = "GeneIntegrationMatrix",
    reducedDims = "IterativeLSI",
    seRNA = seRNA,
    addToArrow = TRUE,
    force= TRUE,
    groupList = groupList,
    groupRNA = "BioClassification",
    nameCell = "predictedCell",
    nameGroup = "predictedGroup",
    nameScore = "predictedScore"
)

现在,当我们使用getAvailableMatrices()检查哪些矩阵可用时,我们会发现GeneIntegrationMatrix已经被添加到Arrow文件中

getAvailableMatrices(projHeme3)
# [1] "GeneIntegrationMatrix" "GeneScoreMatrix" "TileMatrix"

在新的GeneIntegrationMatrix中,我们可以比较连接的基因表达量和根据基因得分推断的基因表达量

我们需要先确保在项目中加入了填充权重值(impute weights)

projHeme3 <- addImputeWeights(projHeme3)

现在,我们来生成一些UMAP图,里面的基因表达量值来自于GeneIntegrationMatrix

markerGenes  <- c(
    "CD34", #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", #B-Cell Trajectory
    "CD14", #Monocytes
    "CD3D", "CD8A", "TBX21", "IL7R" #TCells
  )

p1 <- plotEmbedding(
    ArchRProj = projHeme3, 
    colorBy = "GeneIntegrationMatrix", 
    name = markerGenes, 
    continuousSet = "horizonExtra",
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme3)
)

以及一些相同UMAP图,但是使用GeneScoreMatrix里的基因得分

p2 <- plotEmbedding(
    ArchRProj = projHeme3, 
    colorBy = "GeneScoreMatrix", 
    continuousSet = "horizonExtra",
    name = markerGenes, 
    embedding = "UMAP",
    imputeWeights = getImputeWeights(projHeme3)
)

最后用cowplot将这些标记基因绘制在一起

p1c <- lapply(p1, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})

p2c <- lapply(p2, function(x){
    x + guides(color = FALSE, fill = FALSE) + 
    theme_ArchR(baseSize = 6.5) +
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) +
    theme(
        axis.text.x=element_blank(), 
        axis.ticks.x=element_blank(), 
        axis.text.y=element_blank(), 
        axis.ticks.y=element_blank()
    )
})

do.call(cowplot::plot_grid, c(list(ncol = 3), p1c))
Plot-UMAP-Markers-RNA-W-Imputation_1
do.call(cowplot::plot_grid, c(list(ncol = 3), p2c))
Plot-UMAP-Markers-RNA-W-Imputation_2

和预期的一样,两个方法推测的基因表达量存在相似性,但并不相同。

使用plotPDF()函数保存可编辑的矢量版。

plotPDF(plotList = p1, 
    name = "Plot-UMAP-Marker-Genes-RNA-W-Imputation.pdf", 
    ArchRProj = projHeme3, 
    addDOC = FALSE, width = 5, height = 5)

8.3 使用scRNA-seq信息标记scATAC-seq聚类

现在,我们确定了scATAC-seq和scRNA-seq数据间的对应关系,我们就可以使用scRNA-seq数据中细胞类型对我们的scATAC-seq聚类进行定义。

首先,我们会在scATAC-seq和整合分析得到predictedGroup之间构建一个混合矩阵

cM <- confusionMatrix(projHeme3$Clusters, projHeme3$predictedGroup)
labelOld <- rownames(cM)
labelOld
#  [1] "Cluster11" "Cluster2"  "Cluster12" "Cluster1"  "Cluster8"  "Cluster4" 
#  [7] "Cluster9"  "Cluster5"  "Cluster7"  "Cluster14" "Cluster3"  "Cluster10"
# [13] "Cluster6"  "Cluster13"

接着,对于每一个scATAC-seq聚类,我们根据predictedGroup确定最能定义聚类的细胞类型。

labelNew <- colnames(cM)[apply(cM, 1, which.max)]
labelNew

接着,我们需要对新的聚类标签进行重命名,简化分类系统。对于每一个scRNA-seq的聚类,我们会重新定义标签,以便更好解释。

remapClust <- c(
    "01_HSC" = "Progenitor",
    "02_Early.Eryth" = "Erythroid",
    "03_Late.Eryth" = "Erythroid",
    "04_Early.Baso" = "Basophil",
    "05_CMP.LMPP" = "Progenitor",
    "06_CLP.1" = "CLP",
    "07_GMP" = "GMP",
    "08_GMP.Neut" = "GMP",
    "09_pDC" = "pDC",
    "10_cDC" = "cDC",
    "11_CD14.Mono.1" = "Mono",
    "12_CD14.Mono.2" = "Mono",
    "13_CD16.Mono" = "Mono",
    "15_CLP.2" = "CLP",
    "16_Pre.B" = "PreB",
    "17_B" = "B",
    "18_Plasma" = "Plasma",
    "19_CD8.N" = "CD8.N",
    "20_CD4.N1" = "CD4.N",
    "21_CD4.N2" = "CD4.N",
    "22_CD4.M" = "CD4.M",
    "23_CD8.EM" = "CD8.EM",
    "24_CD8.CM" = "CD8.CM",
    "25_NK" = "NK"
)
remapClust <- remapClust[names(remapClust) %in% labelNew]

接着,使用mapLables()函数进行标签转换,将旧的标签映射到新的标签上。

labelNew2 <- mapLabels(labelNew, oldLabels = names(remapClust), newLabels = remapClust)
labelNew2
#  [1] "GMP"        "B"          "PreB"       "CD4.N"      "Mono"      
#  [6] "Erythroid"  "Progenitor" "CD4.M"      "pDC"        "NK"        
# [11] "CLP"        "Mono"

合并labelsOldlabelsNew2,我们现在可以用mapLables()函数在cellColData里新建聚类标签。

projHeme3$Clusters2 <- mapLabels(projHeme3$Clusters, newLabels = labelNew2, oldLabels = labelOld)

得到新的标签后,我们绘制最新的UMAP展示聚类结果。

p1 <- plotEmbedding(projHeme3, colorBy = "cellColData", name = "Clusters2")
p1

如果被分析scATAC-seq数据对应的细胞系统也有scRNA-seq数据存在,那么这种分析范式将会特别有用。正如之前所说,这种scRNA-seq和scATAC-seq整合分析也为后续更加复杂的基因调控分析提供了出色的框架。

plotPDF()函数保存该图矢量版本。

plotPDF(p1, name = "Plot-UMAP-Remap-Clusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)

使用saveArchRProject保存我们最初的projHeme3.

saveArchRProject(ArchRProj = projHeme3, outputDirectory = "Save-ProjHeme3", load = FALSE)

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×