系列回顧：
ArchR官網(wǎng)教程學習筆記1：Getting Started with ArchR
ArchR官網(wǎng)教程學習筆記2:基于ArchR推測Doublet
ArchR官網(wǎng)教程學習筆記3:創(chuàng)建ArchRProject
ArchR官網(wǎng)教程學習筆記4:ArchR的降維
ArchR官網(wǎng)教程學習筆記5:ArchR的聚類
ArchR官網(wǎng)教程學習筆記6:單細胞嵌入（Single-cell Embeddings）
雖然ArchR能夠穩(wěn)定的call clusters，但它不可能事先知道每個cluster表示哪種細胞類型挠羔。這個任務通常留給手動注釋独旷，因為每個數(shù)據(jù)集都不同牍氛。
為了進行細胞類型注釋窃诉，我們使用細胞類型特異性標記基因的先驗知識漆枚，并通過基因評分從我們的染色質(zhì)可接近數(shù)據(jù)中評估這些基因的基因表達盔夜÷沉ぃ基因評分本質(zhì)上是根據(jù)基因附近調(diào)控元件的可接近性對高表達基因的預測。為了創(chuàng)建這些基因得分旱眯，ArchR允許使用復雜的用戶提供的自定義距離加權接近模型（distance-weighted accessibility models）晨川。
（一）計算基因評分
在我們的發(fā)表的文章中证九，我們測試了50多個不同的基因評分模型，并確定了一類在各種測試條件下始終優(yōu)于其他模型的一種模型共虑。這類模型是ArchR中默認實現(xiàn)的愧怜，它有三個主要組成:
1.在整個基因body中的可接近性決定了基因得分。
2.一種指數(shù)加權函數(shù)妈拌，以距離依賴的方式解釋假定的遠端調(diào)節(jié)元件的活性拥坛。
3.基因邊界，最大限度地減少不相關的調(diào)控元件對基因得分的影響尘分。
那么ArchR是如何計算基因分數(shù)的呢?對于每一個染色體猜惋，ArchR創(chuàng)建一個用戶定義的分片矩陣大小(默認是500bp)，這些分片與用戶定義的基因窗口進行重疊(默認是基因兩側(cè)100bp)培愁。然后計算每個分片到(開始或結束)基因body(上游或下游)或基因開始的距離著摔。我們發(fā)現(xiàn)最佳的基因表達預測是基因區(qū)域的局部可接近性，包括啟動子和gene body定续。如上所述谍咆，為了正確地考慮給定基因的遠端可接近性，ArchR識別出位于基因窗口內(nèi)且不跨越另一個基因區(qū)域的tiles子集私股。這種過濾允許包含遠端調(diào)控元件摹察，可以提高預測基因表達值的準確性，但排除了更可能與其他基因相關的調(diào)控元件(例如庇茫，附近基因的啟動子)。然后使用用戶定義的可接近性模型(默認為e(-abs(distance)/5000) + e-1)將每個分片到基因的距離轉(zhuǎn)換為距離權重螃成。當基因body被包含在基因區(qū)域(基于距離的權重是可能的最大權重)內(nèi)時旦签，我們發(fā)現(xiàn)大的基因會使整體基因評分產(chǎn)生偏差。在這些情況下寸宏，由于內(nèi)含子和外顯子的Tn5插入宁炫，總基因評分會有很大的差異。為了調(diào)整這些基因大小之間存在的差異氮凝，ArchR應用一個單獨的weight：基因大小的倒數(shù)(1 /基因大小)羔巢，并按比例縮放這個權重，從1到一個用戶定義的最大值(默認5)罩阵。更小的基因從而獲得更大的相對權重竿秆，一定程度上減少了長度的影響。然后將對應的距離和基因大小權重乘以每個分片的Tn5插入數(shù)稿壁，并對基因窗口內(nèi)的所有分片進行求和幽钢，同時仍然考慮上面所述的鄰近基因區(qū)域。這個累加的可接近性是一個“基因得分”傅是，并且對所有基因的深度標準化為用戶定義的常量(默認值為10,000)匪燕。計算出的基因評分存儲在相應的Arrow file中蕾羊，以便進行下游分析。
為了說明ArchR基因評分模型是什么樣的帽驯，我們提供了這個示例龟再，展示了在整個基因區(qū)域應用的權重:
如果參數(shù)addGeneScoreMat設置為TRUE(默認)，那么在創(chuàng)建每個Arrow file是就會計算基因評分尼变±眨或者，可以使用addGeneScoreMatrix()函數(shù)隨時將基因評分添加到Arrow files中享甸。一旦被計算截碴，嵌入的單個細胞可以根據(jù)它們的基因評分著色，以幫助識別不同的細胞類型蛉威。我們將在本章剩下的部分闡明基因評分的應用日丹。
值得注意的是，并不是所有的基因都在基因評分里表現(xiàn)得很好（behave well）蚯嫌。特別是哲虾，位于基因密集區(qū)域的基因可能會有問題。因此择示，最好檢查所有的基因評分分析束凑，通過查看測序tracks，這將在后面的章節(jié)中描述栅盲。
（二）鑒定Marker特征
除了使用相關標記基因的先驗知識對clusters進行注釋外汪诉，ArchR還能對任何給定細胞群(例如cluster)的標記特征進行無偏倚的識別。這些特征可以是任何東西——峰谈秫、基因(基于基因評分)或轉(zhuǎn)錄因子motif(基于chromVAR差異)扒寄。ArchR使用getMarkerFeatures()函數(shù)來實現(xiàn)這一點，該函數(shù)可以通過useMatrix參數(shù)接受任何矩陣作為輸入數(shù)據(jù)拟烫，并且它使用groupBy參數(shù)鑒定所指示的組的惟一特性该编。如果useMatrix參數(shù)設置為“GeneScoreMatrix”，那么該功能將識別在每個細胞類型中唯一活躍的基因硕淑。這提供了一種非偏差的方式來觀察哪些基因在每個聚類中被預測為活躍的课竣，并有助于聚類注釋。
如上所述置媳，可以對Arrow files中存儲的任何矩陣使用getMarkerFeatures()函數(shù)來鑒定特定于特定細胞群的特性于樟。這是通過useMatrix參數(shù)完成的。例如拇囊，useMatrix = "TileMatrix"將識別特定細胞群高度特異性的基因組區(qū)域隔披，useMatrix = "PeakMatrix"將識別特定細胞群高度特異性的峰。后面的章節(jié)將提供如何在其他特征類型上使用getMarkerFeatures()函數(shù)的示例寂拆。
Marker特征鑒定是如何發(fā)生的奢米？
Marker特征識別過程依賴于為每個細胞群選擇一組bias-matched的背景細胞抓韩。在所有特征中，將每個細胞群與它自己的背景細胞群進行比較鬓长，以確定給定的細胞群是否具有明顯更高的可接近性谒拴。
這些背景細胞群的選擇對于此過程的成功至關重要，并在用戶提供的多維空間中使用getMarkerFeatures()的bias參數(shù)執(zhí)行涉波。對于細胞群中的每個細胞英上，ArchR在提供的多維空間中查找不屬于給定細胞群的最近臨近細胞，并將其添加到細胞的背景群中啤覆。通過這種方式苍日，ArchR創(chuàng)建了一組與給定細胞群盡可能相似的偏倚匹配細胞，從而使得即使細胞群很小窗声，也能更可靠地確定其顯著性相恃。
ArchR的方法是利用bias參數(shù)提供的所有維度，并對其值進行分位數(shù)標準化笨觅，從而在相同的相對尺度上分布每個維度的variance 拦耐。以下面的示意圖為例，如果將參數(shù)TSS和log10(Num Fragments)提供給bias见剩，分位數(shù)歸一化值之前可能是這樣的:
這里杀糯，y軸上的相對方差與x軸上的相對方差相比很小。如果我們對這些軸進行標準化苍苞，使它們的值范圍從0到1固翰，我們就使相對方差更加相等。重要的是羹呵，我們還大大的改變了最近臨近細胞：
ArchR標準化了所有的維度骂际，并在這個標準化的多維空間中使用歐幾里得距離來找到最臨近的細胞。
（三）鑒定Marker基因
為了根據(jù)基因評分來識別標記基因担巩，我們調(diào)用getMarkerFeatures()函數(shù)方援，并使用useMatrix = "GeneScoreMatrix"没炒。我們通過groupBy = "Clusters"來指定特定于cluster的特征涛癌，它告訴ArchR使用cellColData中的"Clusters"列對細胞群進行分層。
> markersGS <- getMarkerFeatures(
ArchRProj = projHeme2,
useMatrix = "GeneScoreMatrix",
groupBy = "Clusters",
bias = c("TSSEnrichment", "log10(nFrags)"),
testMethod = "wilcoxon"
)
這個函數(shù)返回一個 SummarizedExperiment對象送火，其中包含識別出的標記特征的相關信息拳话。這種類型的返回值在ArchR中很常見，也是ArchR支持下游數(shù)據(jù)分析的關鍵方式之一种吸。SummarizedExperiment對象類似于一個矩陣弃衍，其中行代表感興趣的特征(即基因)，列代表樣本坚俗。一個SummarizedExperiment對象包含一個或多個assay镜盯，每個assay由一個類似矩陣的數(shù)值數(shù)據(jù)表示岸裙，metadata應用于assay矩陣的行或列。如果你需要更多的信息速缆，請查看bioconductor page降允。
我們可以使用getMarkers()函數(shù)得到一系列的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()函數(shù)創(chuàng)建一個熱圖，也可以通過labelMarkers參數(shù)提供一些標記基因來在熱圖上進行標記：
> 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()函數(shù)翅楼，因為heatmapGS對象實際上是一系列的熱圖：
> ComplexHeatmap::draw(heatmapGS, heatmap_legend_side = "bot", annotation_legend_side = "bot")
這里需要注意的是，這一步出圖我在windows里是報錯的真慢，是沒有圖出來的毅臊，我查了一下網(wǎng)上，看來仍然是系統(tǒng)的問題晤碘。后來我切換成Linux系統(tǒng)褂微，再運行一遍就沒有錯了。
https://github.com/GreenleafLab/ArchR/issues/269
> 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這個參數(shù)童社，讓行的順序根據(jù)我們的需求來展示
> 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()函數(shù)中的colorBy和name參數(shù)實現(xiàn)的：
> 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數(shù)據(jù)的稀疏性项栏。我們可以使用MAGIC來分配基因評分使得相鄰細胞的信號更加平滑。在我們的經(jīng)驗中蹬竖，這可以大大提升基因分數(shù)的可視化的解釋沼沈。首先，我們需要添加分配權重到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的基礎上瀏覽這些標記基因的局部染色質(zhì)可接近性旦装。為此页衙，我們使用plotBrowserTrack()函數(shù)，該函數(shù)將創(chuàng)建一個圖的列表，每個圖對應由markerGenes指定的每個基因店乐。此函數(shù)將在groupBy參數(shù)中為每個細胞群繪制單個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數(shù)據(jù)分析的一個挑戰(zhàn)是在組內(nèi)觀察染色質(zhì)可接近性的genome-track水平可視化。傳統(tǒng)上眨八，track可視化需要對scATAC-seq片段進行分組描滔，創(chuàng)建基因組覆蓋的bigwig，并將該track標準化以進行定量可視化踪古。通常含长，終端用戶使用基因組瀏覽器，例如 WashU Epigenome Browser伏穆、UCSC Genome Browser或 IGV browser來可視化這些測序tracks拘泞。這個過程涉及到使用多個軟件，對細胞組的任何更改或添加更多樣品都需要重新生成bigwig文件等枕扫，這可能會非常耗時陪腌。
基于這個原因，ArchR有一個基于Shiny的交互式基因組瀏覽器烟瞧，只需一行代碼ArchRBrowser(ArchRProj)就可以啟動它诗鸭。在Arrow files中實現(xiàn)的數(shù)據(jù)存儲策略允許這個交互式瀏覽器動態(tài)更改細胞組、分辨率和標準化参滴，從而支持實時track-level的可視化强岸。ArchR基因組瀏覽器還以PDF格式創(chuàng)建高質(zhì)量的矢量圖像供發(fā)表。另外砾赔，瀏覽器接受用戶提供的輸入文件蝌箍，比如GenomicRanges對象顯示特征，通過features參數(shù)或者基因組交互文件來定義協(xié)同可接近性暴心、峰to基因的鏈接妓盲，或者通過loop參數(shù)從染色質(zhì)構象數(shù)據(jù)中獲得的loops。對于loops专普，期望的格式是一個GRanges對象悯衬，它的起始位置表示一個loop錨點的中心位置，它的結束位置表示另一個loop錨點的中心位置檀夹。
要啟動本地交互式基因組瀏覽器筋粗，我們使用ArchRBrowser()函數(shù):
> ArchRBrowser(projHeme2)
這時會彈出一個新的窗口：
我們可以通過 “Gene Symbol”來選擇基因:
然后就生成了這個頁面：
你也可以選擇只展示某些clusters：
比如我把6,7,8取消掉，再點擊plot track刷新頁面：
最后保存你的圖：
