系列回顾：
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基因
ArchR官网教程学习笔记8:定义与scRNA-seq一致的聚类
ArchR官网教程学习笔记9:ArchR的伪批量重复
ArchR官网教程学习笔记10:ArchR的call peak
在前面讨论基因评分的章节中，我们已经介绍了标记特征的识别。可以使用相同的函数(getMarkerFeatures())来标识ArchRProject中存储的任何矩阵中的标记特性。标记特征是特定细胞群所特有的特征。这些对于理解特定cluster或细胞类型非常有用。在本章中，我们将讨论如何使用这个功能来识别Marker峰。
（一）鉴定Marker峰
通常，我们感兴趣的是想知道哪些峰是单个cluster或一小群cluster所特有的。在ArchR中，我们可以使用addMarkerFeatures()函数和useMatrix = "PeakMatrix"以无监督的方式完成此任务。
首先，让我们回顾一下我们在projHeme5中研究的细胞类型和它们的相对比例。
#Our scRNA labels
> table(projHeme5$Clusters2)
B CD4.M CD4.N CLP Erythroid GMP Mono
432 639 1279 384 732 1201 2579
NK pDC PreB Progenitor
874 300 358 1473
现在，我们准备使用useMatrix = "PeakMatrix"调用addMarkerFeatures()函数来识别Marker峰。此外，我们通过设置bias参数来考虑TSS富集和每个细胞唯一片段的数量，让ArchR考虑细胞群之间数据质量的差异。
> markersPeaks <- getMarkerFeatures(
ArchRProj = projHeme5,
useMatrix = "PeakMatrix",
groupBy = "Clusters2",
bias = c("TSSEnrichment", "log10(nFrags)"),
testMethod = "wilcoxon"
)
这个对象会通过getMarkerFeatures()返回一个SummarizedExperiment，其中包含一些不同的assays:
> markersPeaks
class: SummarizedExperiment
dim: 140865 11
metadata(2): MatchInfo Params
assays(7): Log2FC Mean ... AUC MeanBGD
rownames(140865): 1 2 ... 140864 140865
rowData names(4): seqnames idx start end
colnames(11): B CD4.M ... PreB Progenitor
colData names(0):
我们可以使用getMarkers()函数检索这个SummarizedExperiment中我们感兴趣的特定片段。此函数的默认行为是返回一个DataFrame对象列表，每个对象是一个细胞群。
> markerList <- getMarkers(markersPeaks, cutOff = "FDR <= 0.01 & Log2FC >= 1")
> markerList
List of length 11
names(11): B CD4.M CD4.N CLP Erythroid GMP Mono NK pDC PreB Progenitor
如果你对某一个细胞群的Marker峰感兴趣，可以单独查看：
> markerList$Erythroid
DataFrame with 2656 rows and 7 columns
seqnames idx start end Log2FC FDR MeanDiff
<Rle> <array> <array> <array> <numeric> <numeric> <numeric>
87111 chr22 1219 30129813 30130313 4.28186227428957 1.83893337171705e-15 1.07953356543409
6670 chr1 6670 110407019 110407519 8.67005424991154 2.96546706376336e-13 0.820868679383739
58340 chr17 6183 73631581 73632081 8.61339322841304 3.53910479048555e-13 0.789176607419672
15006 chr10 1557 30440804 30441304 5.96401972746272 1.43825817219516e-12 0.758436365804101
2590 chr1 2590 27869127 27869627 6.24899843746968 3.09707434920801e-12 0.892089356271252
... ... ... ... ... ... ... ...
104472 chr5 770 32714018 32714518 1.82287547506236 0.0094547934235557 0.126001924074063
125339 chr7 6487 150724392 150724892 2.78589698640991 0.0094547934235557 0.111227644844517
51924 chr16 5054 89045463 89045963 1.75678742468666 0.00948266858646735 0.20433103914127
120570 chr7 1718 30185784 30186284 1.58859209325531 0.00955616601235496 0.326414734260411
83256 chr20 3287 50157615 50158115 1.35032538225079 0.00960141941952872 0.462686019918627
我们可以使用getMarkers()来代替DataFrame对象列表，通过设置returnGR = TRUE来返回GRangesList对象:
> markerList <- getMarkers(markersPeaks, cutOff = "FDR <= 0.01 & Log2FC >= 1", returnGR = TRUE)
> markerList
List of length 11
names(11): B CD4.M CD4.N CLP Erythroid GMP Mono NK pDC PreB Progenitor
这个GRangesList对象也可以查看某一个细胞群的marker峰：
> markerList$Erythroid
GRanges object with 2656 ranges and 3 metadata columns:
seqnames ranges strand | Log2FC FDR MeanDiff
<Rle> <IRanges> <Rle> | <numeric> <numeric> <numeric>
[1] chr22 30129813-30130313 * | 4.28186227428957 1.83893337171705e-15 1.07953356543409
[2] chr1 110407019-110407519 * | 8.67005424991154 2.96546706376336e-13 0.820868679383739
[3] chr17 73631581-73632081 * | 8.61339322841304 3.53910479048555e-13 0.789176607419672
[4] chr10 30440804-30441304 * | 5.96401972746272 1.43825817219516e-12 0.758436365804101
[5] chr1 27869127-27869627 * | 6.24899843746968 3.09707434920801e-12 0.892089356271252
... ... ... ... . ... ... ...
[2652] chr5 32714018-32714518 * | 1.82287547506236 0.0094547934235557 0.126001924074063
[2653] chr7 150724392-150724892 * | 2.78589698640991 0.0094547934235557 0.111227644844517
[2654] chr16 89045463-89045963 * | 1.75678742468666 0.00948266858646735 0.20433103914127
[2655] chr7 30185784-30186284 * | 1.58859209325531 0.00955616601235496 0.326414734260411
[2656] chr20 50157615-50158115 * | 1.35032538225079 0.00960141941952872 0.462686019918627
-------
seqinfo: 23 sequences from an unspecified genome; no seqlengths
（二）绘制Marker峰
（１）热图
我们可以使用markerHeatmap()函数将这些Marker峰(或getMarkerFeatures()输出的任何特征)可视化为一个热图：
> heatmapPeaks <- markerHeatmap(
seMarker = markersPeaks,
cutOff = "FDR <= 0.1 & Log2FC >= 0.5",
transpose = TRUE
)
> reorder_name = c("Mono","Erythroid","GMP","Progenitor","NK","CD4.M","CD4.N","B","PreB","CLP","pDC")
> draw(heatmapPeaks, row_order= reorder_name,heatmap_legend_side = "bot", annotation_legend_side = "bot")
（２）MA Plot和火山图
除了绘制热图，我们还可以为任何单个细胞群绘制一个MA plot或火山图。我们使用markerPlot()函数。对于MA图，我们指定plotAs =“MA”。这里，我们通过name参数指定“Erythroid”细胞群。
> pma <- markerPlot(seMarker = markersPeaks, name = "Erythroid", cutOff = "FDR <= 0.1 & Log2FC >= 1", plotAs = "MA")
> pma
> pv <- markerPlot(seMarker = markersPeaks, name = "Erythroid", cutOff = "FDR <= 0.1 & Log2FC >= 1", plotAs = "Volcano")
> pv
（３）在Browser Tracks展示Marker峰
此外，通过将相关的峰区域传递给plotBrowserTrack()函数中的特性参数，我们可以看到这些峰区域在browser tracks上的覆盖。这将增加一个额外的BED形式的marker峰在ArchR轨迹图的底部。这里我们通过geneSymbol参数绘制GATA1基因:
> p <- plotBrowserTrack(
ArchRProj = projHeme5,
groupBy = "Clusters2",
geneSymbol = c("GATA1"),
features = getMarkers(markersPeaks, cutOff = "FDR <= 0.1 & Log2FC >= 1", returnGR = TRUE)["Erythroid"],
upstream = 50000,
downstream = 50000
)
> grid::grid.newpage()
> grid::grid.draw(p$GATA1)
> plotPDF(p, name = "Plot-Tracks-With-Features", width = 5, height = 5, ArchRProj = projHeme5, addDOC = FALSE)
（三）比较两个细胞群
Marker特征识别是一种非常特殊的差异测试。不过，ArchR还使用相同的getMarkerFeatures()函数支持标准差异测试。这里的小窍门是将useGroups设置为两个细胞组中的一个，而将bgdGroups设置为另一个细胞组。这将在提供的两个组之间执行差异测试。在所有这些差异测试中，传递给useGroups的组中较高的峰将具有正的fold change值，而传递给bgdGroups的组中较高的峰值将具有负的fold change值。
在这里，我们在“Erythroid”细胞组和“Progenitor”细胞组之间进行了配对试验:
> markerTest <- getMarkerFeatures(
ArchRProj = projHeme5,
useMatrix = "PeakMatrix",
groupBy = "Clusters2",
testMethod = "wilcoxon",
bias = c("TSSEnrichment", "log10(nFrags)"),
useGroups = "Erythroid",
bgdGroups = "Progenitor"
)
画MA plot：
> pma <- markerPlot(seMarker = markerTest, name = "Erythroid", cutOff = "FDR <= 0.1 & abs(Log2FC) >= 1", plotAs = "MA")
> pma
画火山图：
> pv <- markerPlot(seMarker = markerTest, name = "Erythroid", cutOff = "FDR <= 0.1 & abs(Log2FC) >= 1", plotAs = "Volcano")
> pv
下一章我们将继续这个差异分析，来寻找差异可接近peaks的motif富集。
