系列回顾：
ArchR官网教程学习笔记1：Getting Started with ArchR
ArchR官网教程学习笔记2:基于ArchR推测Doublet
ArchR官网教程学习笔记3:创建ArchRProject
ArchR官网教程学习笔记4:ArchR的降维
ArchR官网教程学习笔记5:ArchR的聚类
ArchR官网教程学习笔记6:单细胞嵌入（Single-cell Embeddings）
ArchR官网教程学习笔记7:ArchR的基因评分和Marker基因
除了使用基因评分进行clusters鉴定外，ArchR还支持与scRNA-seq的整合。这有助于clusters鉴定分配，因为你可以直接使用在scRNA-seq空间中的clusters，或者在整合后使用基因表达测定。这种整合的工作方式是通过比对scATAC-seq基因评分矩阵和scRNA-seq基因表达矩阵，直接将scRNA-seq的细胞与scATAC-seq的细胞对应上。这种对齐是使用来自Seurat包的FindTransferAnchors()函数执行的，该函数允许你跨两个数据集比对数据。但是，为了适当地将此过程扩展到数十万个细胞，ArchR通过将总细胞划分为更小的细胞组并平行化分别执行比对。
实际上，对于scATAC-seq数据中的每个细胞，这个整合过程是在scRNA-seq数据中查找看起来最相似的细胞，并将该scRNA-seq细胞的基因表达数据分配给scATAC-seq细胞。最后，scATAC-seq空间中的每个细胞都被分配了一个可用于下游分析的基因表达特征。本章介绍了如何使用这些信息来分配clusters，而后面的章节则展示了如何将链接的scRNA-seq数据应用于更复杂的分析，比如识别预测的顺式调控元件。我们相信，随着多组学单细胞的商业化，这些综合分析将变得越来越重要。同时，使用匹配细胞类型的公开可用的scRNA-seq数据，或在你感兴趣的样本上生成的scRNA-seq数据，可以增强在ArchR中执行的scATAC-seq分析。
（一）跨平台将scATAC-seq细胞与scRNA-seq细胞联系
为了将我们的教程scATAC-seq数据与匹配的scRNA-seq数据进行整合，我们将使用Granja*等人(2019)从相同的造血细胞类型中得到的scRNA-seq数据。
我们已经把这个scRNA-seq数据存储为一个大小为111 MB的RangedSummarizedExperiment对象。但是，ArchR也接受未修改的Seurat对象作为整合工作流的输入。下载并检查这个对象，我们看到它有一个基因表达计数矩阵和相关的metadata:
> 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"
)
}
> 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
这个metadata包含了一列，叫BioClassification，包含了scRNA-seq数据集里每个细胞的细胞类型鉴定。
> colnames(colData(seRNA))
[1] "Group" "nUMI_pre" "nUMI"
[4] "nGene" "initialClusters" "UMAP1"
[7] "UMAP2" "Clusters" "BioClassification"
[10] "Barcode"
使用table()可以看一下每一个细胞类型里有多少个细胞：
> table(colData(seRNA)$BioClassification)
01_HSC 02_Early.Eryth 03_Late.Eryth 04_Early.Baso
1425 1653 446 111
05_CMP.LMPP 06_CLP.1 07_GMP 08_GMP.Neut
2260 903 2097 1050
09_pDC 10_cDC 11_CD14.Mono.1 12_CD14.Mono.2
544 325 1800 4222
13_CD16.Mono 14_Unk 15_CLP.2 16_Pre.B
292 520 377 710
17_B 18_Plasma 19_CD8.N 20_CD4.N1
1711 62 1521 2470
21_CD4.N2 22_CD4.M 23_CD8.EM 24_CD8.CM
2364 3539 796 2080
25_NK 26_Unk
2143 161
我们可以执行两种类型的整合方法。无限制整合是一种完全不可知的（agnostic）方法，它采用scATAC-seq实验中的所有细胞，并尝试把它们与scRNA-seq实验中的任何细胞进行比对。虽然这是一个可行的初步解决方案，但我们可以通过限制整合过程来提高跨平台比对的质量。为了执行限制整合，我们使用细胞类型的先验知识来限制比对的搜索空间。举个例子：如果我们知道cluster A, B和C 在scATAC-seq数据对应于三种不同T细胞群，并且我们知细胞群X和Y在 scRNA-seq数据对应于两个不同T细胞群，我们可以告诉ArchR特别尝试使细胞从scATAC-seq里的A,B和C细胞比对scRNA-seq里的X和Y 。下面我们举例说明这两种方法，首先执行无限制整合来获得初步的cluster鉴定，然后使用这些先验知识来执行更精细的限制整合。
（一）无限制整合
为了整合scATAC-seq和scRNA-seq，我们使用了addGeneIntegrationMatrix()函数。如上所述，该函数通过seRNA参数接受Seurat对象或RangedSummarizedExperiment对象。我们执行的第一轮合成是初步的无限制合成，我们不会将其存储在Arrow files中(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"
)
（二）限制整合
现在我们已经完成了初步的无限制整合，我们将识别通用的细胞类型来提供一个框架，从而进一步优化整合结果。
鉴于本教程的数据来自于造血细胞，我们限制整合会很理想，使相似的细胞类型相关联。首先，我们将鉴定从scRNA-seq数据中哪些细胞类型在scATAC-seq cluster中最富集。这样做的目的是使用无限制整合识别scATAC-seq和scRNA-seq数据中与T细胞和NK细胞对应的细胞，以便我们可以使用这些信息执行有限制整合。为此，我们将创建一个confusionMatrix，它查看Clusters和predictedGroup_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
preClust
[1,] "08_GMP.Neut" "C4"
[2,] "21_CD4.N2" "C7"
[3,] "16_Pre.B" "C9"
[4,] "17_B" "C10"
[5,] "11_CD14.Mono.1" "C1"
[6,] "01_HSC" "C5"
[7,] "02_Early.Eryth" "C3"
[8,] "22_CD4.M" "C8"
[9,] "09_pDC" "C11"
[10,] "25_NK" "C6"
[11,] "15_CLP.2" "C12"
[12,] "11_CD14.Mono.1" "C2"
上面的list显示了哪些scRNA-seq细胞类型在12个scATAC-seq 聚类里最富集。
首先，来看一下上面在无限制整合的scRNA-seq数据中的细胞类型标签：
> unique(unique(projHeme2$predictedGroup_Un))
[1] "08_GMP.Neut" "23_CD8.EM" "16_Pre.B" "25_NK"
[5] "17_B" "12_CD14.Mono.2" "11_CD14.Mono.1" "02_Early.Eryth"
[9] "03_Late.Eryth" "07_GMP" "22_CD4.M" "21_CD4.N2"
[13] "05_CMP.LMPP" "09_pDC" "20_CD4.N1" "15_CLP.2"
[17] "06_CLP.1" "01_HSC" "04_Early.Baso" "19_CD8.N"
在上面的列表中，我们可以看到scRNA-seq数据中的对应的T细胞和NK细胞是clusters19-25。我们将创建一个基于字符串表示法用于下游限制整合。
> cTNK <- paste0(paste0(19:25), collapse="|")
> cTNK
[1] "19|20|21|22|23|24|25"
然后我们可以用所有其他clusters，创建一个基于字符串的表示法表示所有“非T细胞，非NK细胞”的clusters。（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"
这些基于字符串的表示法是模式匹配字符串，我们将使用grep提取与这些scRNA-seq细胞类型对应的scATAC-seq clusters。字符串中的|充当or语句，因此我们在confusion matrix的preClust列中搜索与模式匹配字符串中提供的scRNA-seq cluster数字匹配的任何行。
对于T细胞和NK细胞，识别scATAC-seq cluster C6、C7和C8:
> clustTNK <- rownames(cM)[grep(cTNK, preClust)]
> clustTNK
[1] "C7" "C8" "C6"
对于非T细胞和非NK细胞：
> clustNonTNK <- rownames(cM)[grep(cNonTNK, preClust)]
> clustNonTNK
[1] "C4" "C9" "C10" "C1" "C5" "C3" "C11" "C12" "C2"
然后我们执行类似的操作来鉴定对应于这些相同细胞类型的scRNA-seq细胞。首先，我们在scRNA-seq数据中鉴定T细胞和NK细胞：
> 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数据中鉴定非T细胞和非NK细胞：
> 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"
为了准备限制整合，我们创建一个嵌套列表。这是一个包含多个SimpleList对象的SimpleList，每个对象对应我们想要限制的组。在这个例子中，我们有两组：一组称为TNK，在两个平台上识别T细胞和NK细胞，另一组称为NonTNK，在两个平台上识别非T细胞和NK细胞。每个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 files中(' addToArrow = FALSE ')。
#我们建议在将结果保存到Arrow files之前，根据你的期望检查整合的结果。我们将在下一部分说明这个过程。
> 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"
)
（三）比较无限制整合和限制整合
为了比较无限制和有限制的整合，我们将根据通过整合识别的scRNA-seq细胞类型为scATAC-seq数据中的细胞着色。为此，我们将使用内置的ArchR函数paletteDiscrete()创建一个调色板(palette)。
> pal <- paletteDiscrete(values = colData(seRNA)$BioClassification)
Length of unique values greater than palette, interpolating..
在ArchR中，调色板本质上是一个命名的向量，其中的值是十六进制代码，对应于与名称相关联的颜色。
> pal
01_HSC 02_Early.Eryth 03_Late.Eryth 04_Early.Baso
"#D51F26" "#502A59" "#235D55" "#3D6E57"
05_CMP.LMPP 06_CLP.1 07_GMP 08_GMP.Neut
"#8D2B8B" "#DE6C3E" "#F9B712" "#D8CE42"
09_pDC 10_cDC 11_CD14.Mono.1 12_CD14.Mono.2
"#8E9ACD" "#B774B1" "#D69FC8" "#C7C8DE"
13_CD16.Mono 14_Unk 15_CLP.2 16_Pre.B
"#8FD3D4" "#89C86E" "#CC9672" "#CF7E96"
17_B 18_Plasma 19_CD8.N 20_CD4.N1
"#A27AA4" "#CD4F32" "#6B977E" "#518AA3"
21_CD4.N2 22_CD4.M 23_CD8.EM 24_CD8.CM
"#5A5297" "#0F707D" "#5E2E32" "#A95A3C"
25_NK 26_Unk
"#B28D5C" "#3D3D3D"
现在，我们可以通过在无限制整合的基础上覆盖scATAC-seq数据上的scRNA-seq细胞类型来可视化整合：
> p1 <- plotEmbedding(
projHeme2,
colorBy = "cellColData",
name = "predictedGroup_Un",
pal = pal
)
> p1
同样的，我们也可以可视化通过限制整合的基础上覆盖scATAC-seq数据上的scRNA-seq细胞类型：
> p2 <- plotEmbedding(
projHeme2,
colorBy = "cellColData",
name = "predictedGroup_Co",
pal = pal
)
> p2
在本例中，无限制整合和限制整合之间的区别非常微妙，这主要是因为感兴趣的细胞类型已经非常不同了。然而，你应该注意到差异，特别是在T细胞簇(簇17-22)。
保存图片：
> plotPDF(p1,p2, name = "Plot-UMAP-RNA-Integration.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)
保存project:
> saveArchRProject(ArchRProj = projHeme2, outputDirectory = "D:/yanfang/ArchR/Save-ProjHeme2", load = FALSE)
（二）给scATAC-seq数据添加拟scRNA-seq表达谱
现在我们对scATAC-seq和scRNA-seq整合的结果感觉还可以，我们可以使用addToArrow = TRUE重新运行整合的过程，这次我们将链接的基因表达数据添加到Arrow files中。如前所述，我们通过groupList限制整合，并给cellColData的每个元数据nameCell、nameGroup和nameScore添加列名。在这里，我们创建了projHeme3，将在本教程中继续使用。
> 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 files里了：
> getAvailableMatrices(projHeme3)
[1] "GeneIntegrationMatrix" "GeneScoreMatrix" "TileMatrix"
有了这个新的GeneIntegrationMatrix，我们可以比较链接的基因表达和从基因评分中推测的基因表达水平。
首先，我们先要确保在project里添加了填充权重：
> 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把所有的marker基因画出来：
> 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))
> do.call(cowplot::plot_grid, c(list(ncol = 3), p2c))
和我们预期的一样，用两种方法推测的基因表达很相似，但不是完全一样。
保存：
> plotPDF(plotList = p1,
name = "Plot-UMAP-Marker-Genes-RNA-W-Imputation.pdf",
ArchRProj = projHeme3,
addDOC = FALSE, width = 5, height = 5)
（三）使用scRNA-seq信息给scATAC-seq clusters做标记
一旦我们确定了 scATAC-seq 和scRNA-seq的比对是正确的，我们可以使用scRNA-seq数据里的细胞类型给scATAC-seq clusters标记。
首先我们要先创建一个confusion matrix，包括scATAC-seq聚类，以及从整合分析中获得的predictedGroup对象：
> cM <- confusionMatrix(projHeme3$Clusters, projHeme3$predictedGroup)
> labelOld <- rownames(cM)
> labelOld
[1] "C4" "C7" "C9" "C10" "C1" "C5" "C3" "C8" "C11" "C6" "C12" "C2"
对于每一个scATAC-seq cluster，我们从predictedGroup里鉴定细胞类型：
> labelNew <- colnames(cM)[apply(cM, 1, which.max)]
> labelNew
[1] "08_GMP.Neut" "21_CD4.N2" "16_Pre.B" "17_B"
[5] "12_CD14.Mono.2" "01_HSC" "03_Late.Eryth" "22_CD4.M"
[9] "09_pDC" "22_CD4.M" "15_CLP.2" "12_CD14.Mono.2"
接下来，我们需要重新分类这些新的labels，弄一个更简单的分类。对于每个scRNA-seq cluster，我们将重新定义它们的标签，使得更加容易的去解读。
> 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"
)
#将新定义的标签与上面从`predictedGroup`里鉴定细胞类型做个匹配
#把匹配上的存到一个新的变量里
> remapClust <- remapClust[names(remapClust) %in% labelNew]
使用mapLabels()函数，转换我们的标签：
> labelNew2 <- mapLabels(labelNew, oldLabels = names(remapClust), newLabels = remapClust)
> labelNew2
[1] "GMP" "CD4.N" "PreB" "B" "Mono"
[6] "Progenitor" "Erythroid" "CD4.M" "pDC" "CD4.M"
[11] "CLP" "Mono"
结合labelsOld和labelsNew2，我们现在可以使用mapLabels()功能再创建新的cluster标签在 cellColData里：
> projHeme3$Clusters2 <- mapLabels(projHeme3$Clusters, newLabels = labelNew2, oldLabels = labelOld)
用新鉴定的cluster画UMAP：
> p1 <- plotEmbedding(projHeme3, colorBy = "cellColData", name = "Clusters2")
> p1
保存：
> plotPDF(p1, name = "Plot-UMAP-Remap-Clusters.pdf", ArchRProj = projHeme2, addDOC = FALSE, width = 5, height = 5)
保存project3:
> saveArchRProject(ArchRProj = projHeme3, outputDirectory = "Save-ProjHeme3", load = FALSE)
