ArchR官网教程学习笔记7:ArchR的基因评分和Marker

系列回顾：
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能够稳定的call clusters，但它不可能事先知道每个cluster表示哪种细胞类型。这个任务通常留给手动注释，因为每个数据集都不同。

为了进行细胞类型注释，我们使用细胞类型特异性标记基因的先验知识，并通过基因评分从我们的染色质可接近数据中评估这些基因的基因表达。基因评分本质上是根据基因附近调控元件的可接近性对高表达基因的预测。为了创建这些基因得分，ArchR允许使用复杂的用户提供的自定义距离加权接近模型（distance-weighted accessibility models）。

（一）计算基因评分

在我们的发表的文章中，我们测试了50多个不同的基因评分模型，并确定了一类在各种测试条件下始终优于其他模型的一种模型。这类模型是ArchR中默认实现的，它有三个主要组成:
1.在整个基因body中的可接近性决定了基因得分。
2.一种指数加权函数，以距离依赖的方式解释假定的远端调节元件的活性。
3.基因边界，最大限度地减少不相关的调控元件对基因得分的影响。

那么ArchR是如何计算基因分数的呢?对于每一个染色体，ArchR创建一个用户定义的分片矩阵大小(默认是500bp)，这些分片与用户定义的基因窗口进行重叠(默认是基因两侧100bp)。然后计算每个分片到(开始或结束)基因body(上游或下游)或基因开始的距离。我们发现最佳的基因表达预测是基因区域的局部可接近性，包括启动子和gene body。如上所述，为了正确地考虑给定基因的远端可接近性，ArchR识别出位于基因窗口内且不跨越另一个基因区域的tiles子集。这种过滤允许包含远端调控元件，可以提高预测基因表达值的准确性，但排除了更可能与其他基因相关的调控元件(例如，附近基因的启动子)。然后使用用户定义的可接近性模型(默认为e(-abs(distance)/5000) + e-1)将每个分片到基因的距离转换为距离权重。当基因body被包含在基因区域(基于距离的权重是可能的最大权重)内时，我们发现大的基因会使整体基因评分产生偏差。在这些情况下，由于内含子和外显子的Tn5插入，总基因评分会有很大的差异。为了调整这些基因大小之间存在的差异，ArchR应用一个单独的weight：基因大小的倒数(1 /基因大小)，并按比例缩放这个权重，从1到一个用户定义的最大值(默认5)。更小的基因从而获得更大的相对权重，一定程度上减少了长度的影响。然后将对应的距离和基因大小权重乘以每个分片的Tn5插入数，并对基因窗口内的所有分片进行求和，同时仍然考虑上面所述的邻近基因区域。这个累加的可接近性是一个“基因得分”，并且对所有基因的深度标准化为用户定义的常量(默认值为10,000)。计算出的基因评分存储在相应的Arrow file中，以便进行下游分析。

为了说明ArchR基因评分模型是什么样的，我们提供了这个示例，展示了在整个基因区域应用的权重:

如果参数addGeneScoreMat设置为TRUE(默认)，那么在创建每个Arrow file是就会计算基因评分。或者，可以使用addGeneScoreMatrix()函数随时将基因评分添加到Arrow files中。一旦被计算，嵌入的单个细胞可以根据它们的基因评分着色，以帮助识别不同的细胞类型。我们将在本章剩下的部分阐明基因评分的应用。

值得注意的是，并不是所有的基因都在基因评分里表现得很好（behave well）。特别是，位于基因密集区域的基因可能会有问题。因此，最好检查所有的基因评分分析，通过查看测序tracks，这将在后面的章节中描述。

（二）鉴定Marker特征

除了使用相关标记基因的先验知识对clusters进行注释外，ArchR还能对任何给定细胞群(例如cluster)的标记特征进行无偏倚的识别。这些特征可以是任何东西——峰、基因(基于基因评分)或转录因子motif(基于chromVAR差异)。ArchR使用getMarkerFeatures()函数来实现这一点，该函数可以通过useMatrix参数接受任何矩阵作为输入数据，并且它使用groupBy参数鉴定所指示的组的惟一特性。如果useMatrix参数设置为“GeneScoreMatrix”，那么该功能将识别在每个细胞类型中唯一活跃的基因。这提供了一种非偏差的方式来观察哪些基因在每个聚类中被预测为活跃的，并有助于聚类注释。

如上所述，可以对Arrow files中存储的任何矩阵使用getMarkerFeatures()函数来鉴定特定于特定细胞群的特性。这是通过useMatrix参数完成的。例如，useMatrix = "TileMatrix"将识别特定细胞群高度特异性的基因组区域，useMatrix = "PeakMatrix"将识别特定细胞群高度特异性的峰。后面的章节将提供如何在其他特征类型上使用getMarkerFeatures()函数的示例。

Marker特征鉴定是如何发生的？

Marker特征识别过程依赖于为每个细胞群选择一组bias-matched的背景细胞。在所有特征中，将每个细胞群与它自己的背景细胞群进行比较，以确定给定的细胞群是否具有明显更高的可接近性。

这些背景细胞群的选择对于此过程的成功至关重要，并在用户提供的多维空间中使用getMarkerFeatures()的bias参数执行。对于细胞群中的每个细胞，ArchR在提供的多维空间中查找不属于给定细胞群的最近临近细胞，并将其添加到细胞的背景群中。通过这种方式，ArchR创建了一组与给定细胞群尽可能相似的偏倚匹配细胞，从而使得即使细胞群很小，也能更可靠地确定其显著性。

ArchR的方法是利用bias参数提供的所有维度，并对其值进行分位数标准化，从而在相同的相对尺度上分布每个维度的variance 。以下面的示意图为例，如果将参数TSS和log10(Num Fragments)提供给bias，分位数归一化值之前可能是这样的:

这里，y轴上的相对方差与x轴上的相对方差相比很小。如果我们对这些轴进行标准化，使它们的值范围从0到1，我们就使相对方差更加相等。重要的是，我们还大大的改变了最近临近细胞：

ArchR标准化了所有的维度，并在这个标准化的多维空间中使用欧几里得距离来找到最临近的细胞。

（三）鉴定Marker基因

为了根据基因评分来识别标记基因，我们调用getMarkerFeatures()函数，并使用useMatrix = "GeneScoreMatrix"。我们通过groupBy = "Clusters"来指定特定于cluster的特征，它告诉ArchR使用cellColData中的"Clusters"列对细胞群进行分层。

> markersGS <- getMarkerFeatures(
  ArchRProj = projHeme2, 
  useMatrix = "GeneScoreMatrix", 
  groupBy = "Clusters",
  bias = c("TSSEnrichment", "log10(nFrags)"),
  testMethod = "wilcoxon"
)

这个函数返回一个 SummarizedExperiment对象，其中包含识别出的标记特征的相关信息。这种类型的返回值在ArchR中很常见，也是ArchR支持下游数据分析的关键方式之一。SummarizedExperiment对象类似于一个矩阵，其中行代表感兴趣的特征(即基因)，列代表样本。一个SummarizedExperiment对象包含一个或多个assay，每个assay由一个类似矩阵的数值数据表示，metadata应用于assay矩阵的行或列。如果你需要更多的信息，请查看bioconductor page。

我们可以使用getMarkers()函数得到一系列的DataFrame对象，每一个就是我们的clusters，包含相关的marker特征:

> markerList <- getMarkers(markersGS, cutOff = "FDR <= 0.01 & Log2FC >= 1.25")
> markerList$C6 #查看第6个cluster的marker特征
DataFrame with 174 rows and 9 columns
      seqnames     start       end  strand      name     idx    Log2FC         FDR  MeanDiff
         <Rle>   <array>   <array> <array>   <array> <array> <numeric>   <numeric> <numeric>
364       chr1  25291612  25226002       2     RUNX3     364   1.36522 6.35404e-42   1.92850
12597     chr2  87035519  87011728       2      CD8A     538   2.43849 3.16226e-38   1.38007
6856     chr14  99737822  99635625       2    BCL11B     595   1.28029 3.24134e-38   1.27864
19040     chr6 112194655 111981535       2       FYN     881   1.31830 6.27952e-38   1.93826
14759    chr22  37545962  37521880       2     IL2RB     301   1.95915 2.01864e-37   1.20794
...        ...       ...       ...     ...       ...     ...       ...         ...       ...
11796    chr19  53324922  53300661       2     ZNF28    1344   1.26685  0.00717431 0.0562258
14598    chr22  24322019  24313554       2       DDT     140   2.64274  0.00875990 0.0659032
14172    chr21  14778781  14778705       2 MIR3156-3       7   3.73790  0.00912059 0.0248831
6820     chr14  94577079  94583033       1     IFI27     559   1.93959  0.00949383 0.0597439
11785    chr19  52839498  52870376       1    ZNF610    1333   1.28602  0.00968046 0.0994133

为了同时可视化所有标记特征，我们可以使用markerHeatmap()函数创建一个热图，也可以通过labelMarkers参数提供一些标记基因来在热图上进行标记：

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

> heatmapGS <- markerHeatmap(
  seMarker = markersGS, 
  cutOff = "FDR <= 0.01 & Log2FC >= 1.25", 
  labelMarkers = markerGenes,
  transpose = TRUE
)

为了使热图展示出来，我们要使用ComplexHeatmap::draw()函数，因为heatmapGS对象实际上是一系列的热图：

> ComplexHeatmap::draw(heatmapGS, heatmap_legend_side = "bot", annotation_legend_side = "bot")

这里需要注意的是，这一步出图我在windows里是报错的，是没有图出来的，我查了一下网上，看来仍然是系统的问题。后来我切换成Linux系统，再运行一遍就没有错了。
https://github.com/GreenleafLab/ArchR/issues/269

这张图看起来有一些杂乱，不像官网里(https://www.archrproject.com/bookdown/identifying-marker-genes.html)的图看起来有秩序，我们可以把行的顺序重新排列一下，让它看起来更好看一些

> reorder_row = c(3,4,5,11,9,12,1,2,10,7,6,8)
> reorder_row
 [1]  3  4  5 11  9 12  1  2 10  7  6  8
#这里我在原代码里加入了row_order=reorder_row这个参数，让行的顺序根据我们的需求来展示
> ComplexHeatmap::draw(heatmapGS,row_order=reorder_row,heatmap_legend_side = "bot", annotation_legend_side = "bot")

这张图就跟官网的图更像一些了～

保存图片：

> plotPDF(heatmapGS, name = "GeneScores-Marker-Heatmap", width = 8, height = 6, ArchRProj = projHeme2, addDOC = FALSE)

（四）可视化Marker基因

如前所述，我们可以在我们的UMAP嵌入z中覆盖每个细胞的基因评分。使用plotembeds()函数中的colorBy和name参数实现的：

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

> p <- plotEmbedding(
  ArchRProj = projHeme2, 
  colorBy = "GeneScoreMatrix", 
  name = markerGenes, 
  embedding = "UMAP",
  quantCut = c(0.01, 0.95),
  imputeWeights = NULL
)

你可以单独画出某一个子集：

> p$CD14

如果要画出所有基因，我们可以使用cowplot在每一个单独的图里排列不同的marker基因：

> p2 <- lapply(p, 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),p2))

保存图片：

> plotPDF(plotList = p, 
    name = "Plot-UMAP-Marker-Genes-WO-Imputation.pdf", 
    ArchRProj = projHeme2, 
    addDOC = FALSE, width = 5, height = 5)

（五）使用MAGIC分配marker基因

在前面的部分，你可能注意到了，有些基因评分图看起来变化很大。这是因为scATAC-seq数据的稀疏性。我们可以使用MAGIC来分配基因评分使得相邻细胞的信号更加平滑。在我们的经验中，这可以大大提升基因分数的可视化的解释。首先，我们需要添加分配权重到ArchRProject：

> projHeme2 <- addImputeWeights(projHeme2)

当画基因评分的时候，这些分配权重可以被传递给plotEmbedding()和UMAP嵌入重合。

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

> p <- plotEmbedding(
  ArchRProj = projHeme2, 
  colorBy = "GeneScoreMatrix", 
  name = markerGenes, 
  embedding = "UMAP",
  imputeWeights = getImputeWeights(projHeme2)
)
#查看某一个基因
> p$CD14

画出所有marker基因：

> p2 <- lapply(p, 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),p2))
> plotPDF(plotList = p, 
        name = "Plot-UMAP-Marker-Genes-W-Imputation.pdf", 
        ArchRProj = projHeme2, 
        addDOC = FALSE, width = 5, height = 5)

（六）使用 ArchRBrowser画tracks

除了用UMAP覆盖图来绘制每个细胞的基因评分，我们还可以使用genome browser tracks在每个cluster的基础上浏览这些标记基因的局部染色质可接近性。为此，我们使用plotBrowserTrack()函数，该函数将创建一个图的列表，每个图对应由markerGenes指定的每个基因。此函数将在groupBy参数中为每个细胞群绘制单个track。

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

> p <- plotBrowserTrack(
  ArchRProj = projHeme2, 
  groupBy = "Clusters", 
  geneSymbol = markerGenes, 
  upstream = 50000,
  downstream = 50000
)
#画某一个marker基因的tracks
> grid::grid.newpage()
> grid::grid.draw(p$CD14)

保存：

> plotPDF(plotList = p, 
        name = "Plot-Tracks-Marker-Genes.pdf", 
        ArchRProj = projHeme2, 
        addDOC = FALSE, width = 5, height = 5)

（七）启动ArchRBrowser

对scATAC-seq数据分析的一个挑战是在组内观察染色质可接近性的genome-track水平可视化。传统上，track可视化需要对scATAC-seq片段进行分组，创建基因组覆盖的bigwig，并将该track标准化以进行定量可视化。通常，终端用户使用基因组浏览器，例如 WashU Epigenome Browser、UCSC Genome Browser或 IGV browser来可视化这些测序tracks。这个过程涉及到使用多个软件，对细胞组的任何更改或添加更多样品都需要重新生成bigwig文件等，这可能会非常耗时。

基于这个原因，ArchR有一个基于Shiny的交互式基因组浏览器，只需一行代码ArchRBrowser(ArchRProj)就可以启动它。在Arrow files中实现的数据存储策略允许这个交互式浏览器动态更改细胞组、分辨率和标准化，从而支持实时track-level的可视化。ArchR基因组浏览器还以PDF格式创建高质量的矢量图像供发表。另外，浏览器接受用户提供的输入文件，比如GenomicRanges对象显示特征，通过features参数或者基因组交互文件来定义协同可接近性、峰to基因的链接，或者通过loop参数从染色质构象数据中获得的loops。对于loops，期望的格式是一个GRanges对象，它的起始位置表示一个loop锚点的中心位置，它的结束位置表示另一个loop锚点的中心位置。

要启动本地交互式基因组浏览器，我们使用ArchRBrowser()函数:

> ArchRBrowser(projHeme2)

这时会弹出一个新的窗口：

我们可以通过 “Gene Symbol”来选择基因:

然后就生成了这个页面：

你也可以选择只展示某些clusters：

比如我把6,7,8取消掉，再点击plot track刷新页面：

最后保存你的图：