本文首發(fā)于我的個人博客兑凿， http://xuzhougeng.top/

往期回顧:

• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第一章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第二章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第三章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第四章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第五章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第六章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第七章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第八章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第九章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第十章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第十一章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第十二章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第十三章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第十四章)
• 使用ArchR分析單細(xì)胞ATAC-seq數(shù)據(jù)(第十五章)

第16章 使用ArchR進行軌跡分析

ArchR在ArchRProject中低維子空間對細(xì)胞進行排序來創(chuàng)建細(xì)胞軌跡坑夯，實現(xiàn)偽時間中細(xì)胞排序。我們之前已經(jīng)在二維UMAP子空間中執(zhí)行了這種排序馁害，但是ArchR改進了這種方法类垦，使其能夠在N維子空間（即LSI）內(nèi)進行對齊狈邑。

• 首先，ArchR需要用戶提供一個軌跡骨架描述細(xì)胞分組/聚類的大致分化方向蚤认。例如米苹，用戶提供的細(xì)胞編號為HSC, GMP, Monocyte，這表示細(xì)胞以干細(xì)胞作為起點砰琢，然后是祖細(xì)胞驱入，最后是分化細(xì)胞。這種分化方向信息來自于已有的生物學(xué)相關(guān)的細(xì)胞軌跡發(fā)育的背景知識氯析。
• 其次，ArchR計算每個聚類在N維空間的平均坐標(biāo)莺褒，對于聚類里的每個細(xì)胞掩缓，計算其相對于平均指標(biāo)的歐幾里得距離，最后只保留前5%的細(xì)胞遵岩。
• 然后你辣，ArchR根據(jù)軌跡計算cluster i里的每個細(xì)胞和cluster i+1平均坐標(biāo)的距離，根據(jù)距離計算擬時間向量尘执。之后依次遍歷所有的cluster執(zhí)行運算舍哄。根據(jù)每個細(xì)胞相對于細(xì)胞分組/聚類平均坐標(biāo)的歐幾里得距離，ArchR能夠為每個細(xì)胞確定N維坐標(biāo)和它的擬時間值從而成為軌跡的一部分誊锭。
• 然后表悬，ArchR使用smooth.spline函數(shù)根據(jù)擬時間值講連續(xù)的軌跡擬合到N維坐標(biāo)中。
• 接著丧靡，ArchR根據(jù)歐氏距離最近細(xì)胞沿著流形依次將所有細(xì)胞對齊到軌跡上蟆沫。
• 最后籽暇，ArchR將對其后的值縮放到100，并將其保存在ArchRProject中用于下游分析饭庞。

ArchR可以根據(jù)Arrow文件里的特征創(chuàng)建擬時間趨勢矩陣戒悠。例如，ArchR可以分析隨著擬時間發(fā)生變化的TF deviation, 基因得分舟山，整合基因表達(dá)量绸狐，從而鑒定隨著細(xì)胞軌跡動態(tài)變化的調(diào)節(jié)因子或調(diào)控元素。

• 首先累盗，ARchR將細(xì)胞沿著細(xì)胞軌跡進行分組寒矿，組數(shù)由用戶定義（默認(rèn)是1/100）
• 然后，ArchR使用用戶定義平滑滑窗(默認(rèn)是9/100)對矩陣的特征進行平滑處理(調(diào)用data.table::frollmean函數(shù))
• 最后幅骄，ArchR會返回一個SummarizedExperiment對象劫窒，它是一個平滑后的擬時間 X 特征矩陣，用于下游分析拆座。

ArchR還可以根據(jù)以下信息對任意兩個平滑后擬時間 X 特征矩陣進行關(guān)聯(lián)分析主巍，例如匹配的命名(如chromVAR TF deviations的正向調(diào)控因子和基因得分/整合譜)，之前章節(jié)提到的利用低重疊細(xì)胞聚集的基因組位置重疊方法(例如 peak-to-gene 關(guān)聯(lián))挪凑。最終孕索，ArchR能為細(xì)胞軌跡的整合分析提供幫助，解釋多組學(xué)之間動態(tài)調(diào)控關(guān)系躏碳。

16.1 髓系軌跡-單核細(xì)胞分化

這一節(jié)搞旭，我們會創(chuàng)建一個細(xì)胞估計用以模擬HSC發(fā)育成分化后單核細(xì)胞的過程。 在開始之前菇绵，我們先來檢查保存在cellColData中的聚類和之前定義的細(xì)胞類型肄渗，分別是"Clusters"和"Clusters2"。我們通過UMAP來進行展示咬最，從中找到我們感興趣的細(xì)胞類型翎嫡。

p1 <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "Clusters", embedding = "UMAP")
p2 <- plotEmbedding(ArchRProj = projHeme5, colorBy = "cellColData", name = "Clusters2", embedding = "UMAP")
ggAlignPlots(p1, p2, type = "h")

Plot-UMAP-Clusters12-Combined

16.1.1 擬時間UMAP和特征單獨作圖

我們接下來會用到"Clusters2"里面定義的細(xì)胞類型。和之前提到的那樣永乌，我們會以干細(xì)胞("Progenitor")-骨髓祖細(xì)胞("GMP")-單核細(xì)胞("Mono")這個分化過程創(chuàng)建軌跡惑申。

第一步，我們需要以用一個向量記錄軌跡骨干翅雏，里面按順序記錄發(fā)育各過程中細(xì)胞分組的標(biāo)簽圈驼。

trajectory <- c("Progenitor", "GMP", "Mono")

第二步，我們使用addTrajectory()創(chuàng)建軌跡望几，并加入到ArchRProject中绩脆。我們將其命名為"MyeloidU"。該函數(shù)的作用就是在cellColData中新建一列，命名為"MyeloidU"衙伶，然后記錄每個細(xì)胞在軌跡中擬時間值祈坠。不在軌跡中的細(xì)胞記做NA

projHeme5 <- addTrajectory(
    ArchRProj = projHeme5, 
    name = "MyeloidU", 
    groupBy = "Clusters2",
    trajectory = trajectory, 
    embedding = "UMAP", 
    force = TRUE
)

我們可以檢查下該信息，你會發(fā)現(xiàn)每個細(xì)胞都有一個唯一值矢劲，值的取值范圍是0-100赦拘。我們將NA的細(xì)胞過濾掉，因為它們不屬于軌跡芬沉。

head(projHeme5$MyeloidU[!is.na(projHeme5$MyeloidU)])
# [1] 46.07479 53.75027 44.82834 43.18828 47.49617 43.21015

我們使用plotTrajectory()函數(shù)繪制軌跡躺同，它會用UMAP進行可視化，用擬時間值進行上色丸逸，箭頭標(biāo)識軌跡的發(fā)育方向蹋艺。非軌跡的細(xì)胞會以灰色標(biāo)識。在這個例子中黄刚，我們使用colorBy = "cellColData"和name="MyeloidU"讓ArchR使用cellColData里"MyeloidU"作為擬時間軌跡捎谨。參數(shù)中trajectory和name都是"MyeloidU"，可能不太容易理解憔维，ArchR根據(jù)trajectory來提取細(xì)胞涛救，根據(jù)name來對細(xì)胞進行著色。

p <- plotTrajectory(projHeme5, trajectory = "MyeloidU", colorBy = "cellColData", name = "MyeloidU")
p[[1]]

Plot-MyeloidU-Traj-UMAP

使用plotPDF()保存為可編輯的矢量圖

plotPDF(p, name = "Plot-MyeloidU-Traj-UMAP.pdf", ArchRProj = projHeme5, addDOC = FALSE, width = 5, height = 5)

我們還可以在UMAP嵌入圖中展示軌跡和其他特征业扒。這可以用來展示和我們軌跡有關(guān)的特征检吆。

如果你還沒有在projHeme5項目中加入推測權(quán)重，請先運行下面的代碼

projHeme5 <- addImputeWeights(projHeme5)

然后我們繪制"MyeloidU"軌跡程储，只不過顏色來自于"CEBPB"的基因得分蹭沛，這是一個已知的調(diào)控單核細(xì)胞功能的基因，會隨著發(fā)育而活躍章鲤。我們通過colorBy參數(shù)來指定表達(dá)量矩陣摊灭，使用name來指定特征

p1 <- plotTrajectory(projHeme5, trajectory = "MyeloidU", colorBy = "GeneScoreMatrix", name = "CEBPB", continuousSet = "horizonExtra")

重復(fù)上面的代碼，只不過這次使用的是GeneIntegrationMatrix作為基因表達(dá)量矩陣

p2 <- plotTrajectory(projHeme5, trajectory = "MyeloidU", colorBy = "GeneIntegrationMatrix", name = "CEBPB", continuousSet = "blueYellow")

plotTrajectory函數(shù)實際上返回了一系列圖败徊。列表中第一個圖是UMAP嵌入圖帚呼，根據(jù)函數(shù)調(diào)用進行上色。

通過比較基因得分和基因表達(dá)量的UMAP圖集嵌，不難發(fā)現(xiàn)CEBPB基因在單核細(xì)胞中特異性表達(dá)，出現(xiàn)在你時序軌跡的后半部分御毅。

ggAlignPlots(p1[[1]], p2[[1]], type = "h")

Plot-UMAP-CEBPB-Combined

plotTrajectory()函數(shù)返回的是一個點圖根欧，x軸是擬時間，y軸是對應(yīng)特征的值端蛆，此處是CEBPB的基因得分或基因表達(dá)量凤粗。下圖中的細(xì)胞以擬時間進行著色

ggAlignPlots(p1[[2]], p2[[2]], type = "h")

Plot-UMAP-CEBPB-Combined2

16.1.2 擬時間熱圖

我們可以用熱圖的形式對特征隨著擬時間變化進行可視化。首先我們用getTrajectory()函數(shù)從ArchRProject中提取感興趣的軌跡，它會返回一個SummarizedExperiment對象嫌拣。通過設(shè)置useMatrix參數(shù)柔袁，我們可以用motif，基因得分异逐，基因表達(dá)量捶索，peak開放性來創(chuàng)建擬時間熱圖

trajMM  <- getTrajectory(ArchRProj = projHeme5, name = "MyeloidU", useMatrix = "MotifMatrix", log2Norm = FALSE)

接著將SummarizedExperiment傳給plotTrajectoryHeatmap()函數(shù)

p1 <- plotTrajectoryHeatmap(trajMM, pal = paletteContinuous(set = "solarExtra"))

Plot-MyeloidU-Traj-Heatmaps

同樣的分析可以使用基因得分繪制擬時間熱圖，只需要設(shè)置useMatrix = "GeneScoreMatrix"

trajGSM <- getTrajectory(ArchRProj = projHeme5, name = "MyeloidU", useMatrix = "GeneScoreMatrix", log2Norm = TRUE)
p2 <- trajectoryHeatmap(trajGSM,  pal = paletteContinuous(set = "horizonExtra"))

Plot-MyeloidU-Traj-Heatmaps_2

同樣的灰瞻，還可以繪制基因表達(dá)量的擬時間熱圖腥例，設(shè)置useMatrix = "GeneIntegrationMatrix".

trajGIM <- getTrajectory(ArchRProj = projHeme5, name = "MyeloidU", useMatrix = "GeneIntegrationMatrix", log2Norm = FALSE)
p3 <- plotTrajectoryHeatmap(trajGIM,  pal = paletteContinuous(set = "blueYellow"))

Plot-MyeloidU-Traj-Heatmaps_3

最后，我們可以通過設(shè)置useMatrix = "PeakMatrix"繪制peak開放性的擬時間熱圖

trajPM  <- getTrajectory(ArchRProj = projHeme5, name = "MyeloidU", useMatrix = "PeakMatrix", log2Norm = TRUE)
p4 <- plotTrajectoryHeatmap(trajPM, pal = paletteContinuous(set = "solarExtra"))

Plot-MyeloidU-Traj-Heatmaps_4

使用plotPDF()保存可編輯的矢量圖

plotPDF(p1, p2, p3, p4, name = "Plot-MyeloidU-Traj-Heatmaps.pdf", ArchRProj = projHeme5, addDOC = FALSE, width = 6, height = 8)

16.1.3 整合擬時間分析

我們也可以進行整合分析酝润，例如整合基因得分/基因表達(dá)量和motif開發(fā)狀態(tài)燎竖，鑒定隨擬時間變化的正向TF調(diào)控因子。這是一個非常有效的手段要销，例如鑒定分化相關(guān)的驅(qū)動因子构回。我們可以用correlateTrajectories()函數(shù)進行該分析，它接受兩個SummarizedExperiment對象疏咐，該對象可以用getTrajectories()函數(shù)獲取纤掸。

首先，讓我們找到開放狀態(tài)隨擬時間和TF基因的基因得分相關(guān)的motif

corGSM_MM <- correlateTrajectories(trajGSM, trajMM)

correlateTrajectories()函數(shù)的主要輸出是一個列表凳鬓，第一個是DataFrame對象茁肠。DataFrame的列名分別是idx1, matchname1, name1和VarAssay1, 對應(yīng)索引，匹配名缩举，原名垦梆，傳遞給correlateTrajectories()函數(shù)的第一個軌跡的特征的方差分位數(shù)(variance quantile)。方差分位數(shù)是給定特征的標(biāo)準(zhǔn)化值仅孩，能用于比較不同來源assay里的相關(guān)值托猩。該DataFrame包括所有符合correlateTrajectories()函數(shù)里閾值的特征。

corGSM_MM[[1]]
## DataFrame with 36 rows and 12 columns
## idx1 idx2 matchname1 matchname2 name1 name2
##
## 1 82 1081 PRDM16 PRDM16 chr1:PRDM16 z:PRDM16_211
## … … … … … … …
## 36 22097 1565 RXRA RXRA chr9:RXRA z:RXRA_695
## Correlation VarAssay1 VarAssay2 TStat
##
## 1 0.859588836579798 0.999783802481948 0.860344827586207 16.6530781620713
## … … … … …
## 36 0.59509316489084 0.88640982401522 0.890229885057471 7.33039570155334
## Pval FDR
##
## 1 2.47466579579935e-30 4.04855324192774e-27
## … … …
## 36 6.60626612679289e-11 1.25672690505037e-09

我們?nèi)缓髲?code>SummarizedExperiment提取相應(yīng)的軌跡辽慕，

trajGSM2 <- trajGSM[corGSM_MM[[1]]$name1, ]
trajMM2 <- trajMM[corGSM_MM[[1]]$name2, ]

為了更好的對這些特征排序京腥，我們創(chuàng)建新的軌跡，它的值來自于原先的軌跡的相乘溅蛉。這可以讓我們的熱圖根據(jù)行進行排序

trajCombined <- trajGSM2
assay(trajCombined) <- t(apply(assay(trajGSM2), 1, scale)) + t(apply(assay(trajMM2), 1, scale))

我們可以從plotTrajectoryHeatmap()函數(shù)返回的結(jié)果提取最優(yōu)的行序公浪。

combinedMat <- plotTrajectoryHeatmap(trajCombined, returnMat = TRUE, varCutOff = 0)
rowOrder <- match(rownames(combinedMat), rownames(trajGSM2))

有了這個之后，我們就可以創(chuàng)建配對熱圖船侧。

首先欠气，我們根據(jù)基因得分軌跡創(chuàng)建熱圖。在rowOrder參數(shù)中設(shè)置給定的行序

ht1 <- plotTrajectoryHeatmap(trajGSM2,  pal = paletteContinuous(set = "horizonExtra"),  varCutOff = 0, rowOrder = rowOrder)

然后镜撩，我們根據(jù)motif軌跡創(chuàng)建熱圖预柒。在rowOrder參數(shù)中設(shè)置給定的行序

ht2 <- plotTrajectoryHeatmap(trajMM2, pal = paletteContinuous(set = "solarExtra"), varCutOff = 0, rowOrder = rowOrder)

通過將兩個熱圖并排排列，我們可以看到兩個熱圖的行是匹配的。你可能會發(fā)現(xiàn)該分析同時捕獲了GATA3和GATA3-AS1（一個GATA3的反義轉(zhuǎn)錄本）宜鸯。這是由于特征匹配所導(dǎo)致的憔古，我們需要在后續(xù)分析時進行手動過濾。

ht1 + ht2

Plot-MyeloidU-GS-TF

我們可以重復(fù)以上的分析淋袖，但這次使用GeneIntegrationMatrix而非基因得分鸿市。因為這是相同的分析流程，因此代碼就不再過多解釋

corGIM_MM <- correlateTrajectories(trajGIM, trajMM)
trajGIM2 <- trajGIM[corGIM_MM[[1]]$name1, ]
trajMM2 <- trajMM[corGIM_MM[[1]]$name2, ]

trajCombined <- trajGIM2
assay(trajCombined) <- t(apply(assay(trajGIM2), 1, scale)) + t(apply(assay(trajMM2), 1, scale))

combinedMat <- plotTrajectoryHeatmap(trajCombined, returnMat = TRUE, varCutOff = 0)

rowOrder <- match(rownames(combinedMat), rownames(trajGIM2))

ht1 <- plotTrajectoryHeatmap(trajGIM2,  pal = paletteContinuous(set = "blueYellow"),  varCutOff = 0, rowOrder = rowOrder)
ht2 <- plotTrajectoryHeatmap(trajMM2, pal = paletteContinuous(set = "solarExtra"), varCutOff = 0, rowOrder = rowOrder)

ht1 + ht2

Plot-MyeloidU-GEx-T

16.2 淋巴軌跡-B細(xì)胞分化

第二個例子會根據(jù)祖細(xì)胞-淋巴祖細(xì)胞-前B細(xì)胞-分化后B細(xì)胞來構(gòu)建B細(xì)胞的分化軌跡适贸。由于該分析從本質(zhì)上是重復(fù)之前單核細(xì)胞軌跡分析灸芳，因此這里只介紹關(guān)鍵性的代碼不同。

第一處是trajectory變量的賦值拜姿。

trajectory <- c("Progenitor", "CLP", "PreB", "B")

第二處是使用addTrajectory()創(chuàng)建軌跡烙样，參數(shù)name=LymphoidU。

projHeme5 <- addTrajectory(
    ArchRProj = projHeme5, 
    name = "LymphoidU", 
    groupBy = "Clusters2",
    trajectory = trajectory, 
    embedding = "UMAP", 
    force = TRUE
)

第三處, 我們需要注意plotTrajectory里的參數(shù)trajectory和name需要改成"LymphoidU"

p <- plotTrajectory(projHeme5, trajectory = "LymphoidU", colorBy = "cellColData", name = "LymphoidU")

后續(xù)的代碼基本上可以照搬上一節(jié)蕊肥。