第八節(jié)：擬時(shí)序細(xì)胞軌跡推斷

在本節(jié)教程中腌乡，我們將學(xué)習(xí)如何通過(guò)擬時(shí)序分析推斷細(xì)胞分化軌跡炼团。slingshot包可以對(duì)單細(xì)胞RNA-seq數(shù)據(jù)進(jìn)行細(xì)胞分化譜系構(gòu)建和偽時(shí)間推斷，它利用細(xì)胞聚類(lèi)簇和空間降維信息，以無(wú)監(jiān)督或半監(jiān)督的方式學(xué)習(xí)細(xì)胞聚類(lèi)群之間的關(guān)系诞丽，揭示細(xì)胞聚類(lèi)簇之間的全局結(jié)構(gòu)眉菱，并將該結(jié)構(gòu)轉(zhuǎn)換為由一維變量表示的平滑譜系迹栓，稱(chēng)之為“偽時(shí)間”。

運(yùn)行slingshot的至少需要兩個(gè)輸入文件：即細(xì)胞在降維空間中的坐標(biāo)矩陣和細(xì)胞聚類(lèi)群的標(biāo)簽向量俭缓。通過(guò)這兩個(gè)輸入文件克伊，我們可以：

• 使用getLineages函數(shù)在細(xì)胞聚類(lèi)群上構(gòu)建最小生成樹(shù)（MST），確定細(xì)胞全局的譜系結(jié)構(gòu)华坦；
• 利用getCurves函數(shù)擬合主曲線來(lái)構(gòu)造平滑譜系愿吹，并推斷偽時(shí)間變量；
• 使用內(nèi)置的可視化工具評(píng)估不同步驟的分析結(jié)果惜姐。

安裝并加載所需的R包

BiocManager::install("slingshot")
library(slingshot)
library(SingleCellExperiment)

構(gòu)建測(cè)試數(shù)據(jù)集

這里犁跪，我們模擬了兩個(gè)不同形式的數(shù)據(jù)集進(jìn)行譜系推斷分析椿息。
接下來(lái)，我們構(gòu)建第一個(gè)測(cè)試數(shù)據(jù)集坷衍，用于表示單個(gè)譜系寝优，其中三分之一的基因與過(guò)渡相關(guān)。我們將該數(shù)據(jù)集包含在一個(gè)SingleCellExperiment對(duì)象中枫耳，并將其用于演示相應(yīng)的分析流程乏矾。

# generate synthetic count data representing a single lineage
# 構(gòu)建單細(xì)胞表達(dá)矩陣
means <- rbind(
    # non-DE genes
    matrix(rep(rep(c(0.1,0.5,1,2,3), each = 300),100),
        ncol = 300, byrow = TRUE),
    # early deactivation
    matrix(rep(exp(atan( ((300:1)-200)/50 )),50), ncol = 300, byrow = TRUE),
    # late deactivation
    matrix(rep(exp(atan( ((300:1)-100)/50 )),50), ncol = 300, byrow = TRUE),
    # early activation
    matrix(rep(exp(atan( ((1:300)-100)/50 )),50), ncol = 300, byrow = TRUE),
    # late activation
    matrix(rep(exp(atan( ((1:300)-200)/50 )),50), ncol = 300, byrow = TRUE),
    # transient
    matrix(rep(exp(atan( c((1:100)/33, rep(3,100), (100:1)/33) )),50), 
        ncol = 300, byrow = TRUE)
)

counts <- apply(means,2,function(cell_means){
    total <- rnbinom(1, mu = 7500, size = 4)
    rmultinom(1, total, cell_means)
})
rownames(counts) <- paste0('G',1:750)
colnames(counts) <- paste0('c',1:300)
counts[1:5,1:5]
#   c1 c2 c3 c4 c5
#G1  0  2  0  0  0
#G2  5  2  2  2  9
#G3  5  9  3  5  8
#G4  6  7  5 12 23
#G5 14 21 12 17 20
dim(counts)
#[1] 750 300

# 將表達(dá)矩陣轉(zhuǎn)換為SingleCellExperiment對(duì)象
sim <- SingleCellExperiment(assays = List(counts = counts))
sim
#class: SingleCellExperiment 
#dim: 750 300 
#metadata(0):
#assays(1): counts
#rownames(750): G1 G2 ... G749 G750
#rowData names(0):
#colnames(300): c1 c2 ... c299 c300
#colData names(0):
#reducedDimNames(0):
#spikeNames(0):

接下來(lái)，我們構(gòu)建第二個(gè)示例數(shù)據(jù)集迁杨，該數(shù)據(jù)集由降維空間坐標(biāo)矩陣（如通過(guò)PCA, ICA和diffusion map等獲得）和細(xì)胞聚類(lèi)群標(biāo)簽（如K-means聚類(lèi)等生成）組成钻心。該數(shù)據(jù)集表示一個(gè)分叉的軌跡，可以用來(lái)演示slingshot包的其他功能铅协。

# 加載slingshot包內(nèi)置數(shù)據(jù)
data("slingshotExample")
rd <- slingshotExample$rd
cl <- slingshotExample$cl

dim(rd) # data representing cells in a reduced dimensional space
## [1] 140   2
length(cl) # vector of cluster labels
## [1] 140

head(rd)
#           dim-1      dim-2
#cell-1 -4.437117  0.3690586
#cell-2 -3.837149 -0.1009774
#cell-3 -2.730936 -0.9770827
#cell-4 -4.219531 -0.7100939
#cell-5 -3.717672 -0.6452345
#cell-6 -3.053360 -0.2838916
head(cl)
#cell-1 cell-2 cell-3 cell-4 cell-5 cell-6 
#    1      1      1      1      1      1 

上游分析

基因過(guò)濾

要開(kāi)始對(duì)單個(gè)譜系數(shù)據(jù)集進(jìn)行分析捷沸，我們需要降低數(shù)據(jù)的維數(shù)，并過(guò)濾掉無(wú)信息的基因狐史。這將大大提高下游分析的速度痒给，同時(shí)將信息損失降至最低。

對(duì)于基因過(guò)濾步驟预皇，我們保留了至少在10個(gè)細(xì)胞中表達(dá)侈玄，且基因表達(dá)count數(shù)至少大于3的基因。

# filter genes down to potential cell-type markers
# at least M (15) reads in at least N (15) cells
geneFilter <- apply(assays(sim)$counts,1,function(x){
    sum(x >= 3) >= 10
})
sim <- sim[geneFilter, ]

數(shù)據(jù)標(biāo)準(zhǔn)化

由于我們使用的是模擬數(shù)據(jù)吟温，因此無(wú)需擔(dān)心批處理效應(yīng)或其他潛在的混雜因素序仙。在這里，我們將使用分位數(shù)歸一化的方法對(duì)表達(dá)數(shù)據(jù)進(jìn)行歸一化處理鲁豪，它可以使每個(gè)細(xì)胞具有相同的表達(dá)值分布潘悼。

FQnorm <- function(counts){
    rk <- apply(counts,2,rank,ties.method='min')
    counts.sort <- apply(counts,2,sort)
    refdist <- apply(counts.sort,1,median)
    norm <- apply(rk,2,function(r){ refdist[r] })
    rownames(norm) <- rownames(counts)
    return(norm)
}

assays(sim)$norm <- FQnorm(assays(sim)$counts)

數(shù)據(jù)降維

slingshot基本假設(shè)是，轉(zhuǎn)錄相似的細(xì)胞在某些降維空間中會(huì)彼此靠近爬橡。由于我們?cè)跇?gòu)造譜系和測(cè)量偽時(shí)間時(shí)使用歐幾里得距離治唤，因此對(duì)數(shù)據(jù)進(jìn)行降維處理很重要。

在這里糙申，我們將演示兩種降維方法：主成分分析（PCA）和均勻流形逼近和投影（UMAP宾添，通過(guò)uwot包調(diào)用）。

# 使用prcomp函數(shù)進(jìn)行PCA降維
pca <- prcomp(t(log1p(assays(sim)$norm)), scale. = FALSE)
rd1 <- pca$x[,1:2]
head(rd1)
 #        PC1       PC2
#c1 -15.55731 -6.218455
#c2 -14.48979 -7.060021
#c3 -16.10326 -7.630399
#c4 -15.59147 -7.443861
#c5 -15.86607 -6.180724
#c6 -15.78163 -6.733366

plot(rd1, col = rgb(0,0,0,.5), pch=16, asp = 1)

image.png

library(uwot)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following object is masked from 'package:S4Vectors':
## 
##     expand
# 進(jìn)行UMAP降維
rd2 <- umap(t(log1p(assays(sim)$norm)))
colnames(rd2) <- c('UMAP1', 'UMAP2')
head(rd2)
#        UMAP1     UMAP2
#[1,] 6.357223 -4.358757
#[2,] 6.442034 -4.216842
#[3,] 6.478745 -4.421883
#[4,] 6.712623 -4.487934
#[5,] 6.497382 -4.276127
#[6,] 6.302136 -4.531643

plot(rd2, col = rgb(0,0,0,.5), pch=16, asp = 1)

image.png

接下來(lái)柜裸，我們可以將降維的結(jié)果添加到SingleCellExperiment對(duì)象中缕陕。

reducedDims(sim) <- SimpleList(PCA = rd1, UMAP = rd2)
sim
#class: SingleCellExperiment 
#dim: 734 300
#metadata(0):
#assays(2): counts norm
#rownames(734): G2 G3 ... G749 G750
#rowData names(0):
#colnames(300): c1 c2 ... c299 c300
#colData names(0): 
#reducedDimNames(2): PCA UMAP
#spikeNames(0):

細(xì)胞聚類(lèi)

slingshot的最終輸入是細(xì)胞聚類(lèi)簇的標(biāo)簽向量。如果未提供此參數(shù)疙挺，slingshot則將數(shù)據(jù)視為單個(gè)細(xì)胞簇并擬合標(biāo)準(zhǔn)主曲線扛邑。但是，我們建議即使在僅期望單一譜系的數(shù)據(jù)集中也應(yīng)該對(duì)細(xì)胞進(jìn)行聚類(lèi)分群铐然，因?yàn)樗梢詽撛诘匕l(fā)現(xiàn)新的分支事件蔬崩。

在此步驟中得到的聚類(lèi)群集將用于確定潛在分化譜系的全局結(jié)構(gòu)恶座。這與聚類(lèi)單細(xì)胞數(shù)據(jù)的典型目標(biāo)不同，后者的主要目的是識(shí)別數(shù)據(jù)集中存在的所有生物學(xué)相關(guān)的細(xì)胞類(lèi)型沥阳。例如跨琳，在確定全局譜系結(jié)構(gòu)時(shí)，無(wú)需區(qū)分未成熟神經(jīng)元和成熟神經(jīng)元桐罕，因?yàn)閮煞N細(xì)胞類(lèi)型可能會(huì)沿著譜系的同一段落入湾宙。

在這里，我們采用了兩種聚類(lèi)方法：高斯混合建模和 k-means聚類(lèi)冈绊，它們類(lèi)似地假設(shè)低維空間中的歐幾里得距離反映了細(xì)胞之間的生物學(xué)差。前者可以使用mclust包進(jìn)行實(shí)現(xiàn)（Scrucca et al埠啃，2016年）死宣，它可以提供一種基于貝葉斯信息準(zhǔn)則（BIC）的方法自動(dòng)確定細(xì)胞聚類(lèi)的簇?cái)?shù)。

library(mclust, quietly = TRUE)
## Package 'mclust' version 5.4.6
## Type 'citation("mclust")' for citing this R package in publications.
## 
## Attaching package: 'mclust'
## The following object is masked from 'package:mgcv':
## 
##     mvn
# 使用Mclust函數(shù)進(jìn)行細(xì)胞聚類(lèi)
cl1 <- Mclust(rd1)$classification
head(cl1)
#c1 c2 c3 c4 c5 c6 
# 1  1  1  1  1  1 
colData(sim)$GMM <- cl1
sim
#class: SingleCellExperiment 
#dim: 734 300 
#metadata(0):
#assays(2): counts norm
#rownames(734): G2 G3 ... G749 G750
#rowData names(0):
#colnames(300): c1 c2 ... c299 c300
#colData names(1): GMM
#reducedDimNames(2): PCA UMAP
#spikeNames(0):

# 可視化聚類(lèi)分群的結(jié)果
library(RColorBrewer)
plot(rd1, col = brewer.pal(9,"Set1")[cl1], pch=16, asp = 1)

image.png

k-means聚類(lèi)沒(méi)有類(lèi)似的功能碴开，因此我們需要人為的指定聚類(lèi)的細(xì)胞簇?cái)?shù)毅该。如果細(xì)胞簇?cái)?shù)選擇的太低，我們可能會(huì)錯(cuò)過(guò)一個(gè)真正的分支事件潦牛；如果它太高或有很多小簇眶掌，我們可能會(huì)開(kāi)始看到虛假的分支事件。這里巴碗，我們選擇聚類(lèi)簇?cái)?shù)為4進(jìn)行k-means聚類(lèi)分析朴爬。

# 使用kmeans函數(shù)進(jìn)行細(xì)胞聚類(lèi)
cl2 <- kmeans(rd1, centers = 4)$cluster
head(cl2)
#c1 c2 c3 c4 c5 c6 
# 2  2  2  2  2  2 
colData(sim)$kmeans <- cl2

# 可視化聚類(lèi)分群的結(jié)果
plot(rd1, col = brewer.pal(9,"Set1")[cl2], pch=16, asp = 1)

image.png

使用slingshot構(gòu)建分化譜系并進(jìn)行擬時(shí)推斷

現(xiàn)在，我們已經(jīng)得到了slingshot分析所要用的兩個(gè)文件：細(xì)胞降維的空間坐標(biāo)矩陣和聚類(lèi)分群的標(biāo)簽向量橡淆。接下來(lái)召噩，我們就可以運(yùn)行slinghsot進(jìn)行偽時(shí)間推斷分析。這是一個(gè)兩步的過(guò)程逸爵，包括用基于細(xì)胞群集的最小生成樹(shù)（MST）識(shí)別全局譜系結(jié)構(gòu)具滴，以及擬合主要曲線來(lái)描述每個(gè)譜系。

這兩個(gè)步驟可以使用getLineages和getCurves函數(shù)分開(kāi)運(yùn)行师倔，也可以使用slingshot封裝函數(shù)一起運(yùn)行（推薦）构韵。這里，我們將使用封裝函數(shù)對(duì)單軌跡數(shù)據(jù)集進(jìn)行譜系分析趋艘，稍后將在分叉數(shù)據(jù)集上演示分步函數(shù)的用法疲恢。

為了使用PCA降維產(chǎn)生的空間坐標(biāo)信息和高斯混合建模識(shí)別的細(xì)胞聚類(lèi)簇標(biāo)簽來(lái)運(yùn)行slingshot，我們將執(zhí)行以下操作：

sim <- slingshot(sim, clusterLabels = 'GMM', reducedDim = 'PCA')
## Using full covariance matrix
sim
class: SingleCellExperiment 
dim: 734 300 
metadata(0):
assays(2): counts norm
rownames(734): G2 G3 ... G749 G750
rowData names(0):
colnames(300): c1 c2 ... c299 c300
colData names(3): GMM kmeans slingPseudotime_1
reducedDimNames(2): PCA UMAP
spikeNames(0):

如上所述致稀，如果未提供聚類(lèi)分群信息冈闭，slingshot則假定所有細(xì)胞都是同一個(gè)聚類(lèi)群，并且將構(gòu)建一條擬合曲線抖单。如果未提供降維空間坐標(biāo)信息萎攒，slingshot將使用由返回的列表的第一個(gè)元素reducedDims中的信息遇八。

輸出的結(jié)果是一個(gè)包含了slingshot結(jié)果的SingleCellExperiment對(duì)象。大多數(shù)的輸出以列表的形式添加到元數(shù)據(jù)中耍休，并可以通過(guò)metadata(sce)$slingshot來(lái)進(jìn)行訪問(wèn)刃永。此外，所有推斷的偽時(shí)間變量（每個(gè)譜系一個(gè)）都添加到colData中羊精。如果要提取slingshot單個(gè)對(duì)象中的所有結(jié)果斯够，我們可以使用SlingshotDataSet函數(shù)。這對(duì)于可視化很有用喧锦，因?yàn)镾lingshotDataSet函數(shù)中包含了幾種對(duì)象繪制方法 读规。下面，我們對(duì)單軌跡數(shù)據(jù)譜系推斷的結(jié)果進(jìn)行可視化處理燃少，并用偽時(shí)間對(duì)細(xì)胞進(jìn)行著色束亏。

summary(sim$slingPseudotime_1)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   8.633  21.118  21.415  34.367  43.186
colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plotcol <- colors[cut(sim$slingPseudotime_1, breaks=100)]

plot(reducedDims(sim)$PCA, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(sim), lwd=2, col='black')

image.png

通過(guò)使用type參數(shù)，還可以看到基于細(xì)胞聚類(lèi)群的最小生成樹(shù)是如何推斷出譜系結(jié)構(gòu)的阵具。

plot(reducedDims(sim)$PCA, col = brewer.pal(9,'Set1')[sim$GMM], pch=16, asp = 1)
lines(SlingshotDataSet(sim), lwd=2, type = 'lineages', col = 'black')

image.png

下游分析

鑒定隨偽時(shí)間表達(dá)變化的基因

運(yùn)行完slingshot推斷出細(xì)胞分化譜系后碍遍，我們可以進(jìn)一步鑒定出隨偽時(shí)間或發(fā)育過(guò)程中表達(dá)變化的基因。我們將使用tradeSeq軟件包演示這種類(lèi)型的分析（Van den Berge等阳液，2020）怕敬。

對(duì)于每個(gè)基因，我們將使用負(fù)二項(xiàng)式噪聲分布來(lái)擬合通用加性模型（GAM）帘皿，以對(duì)基因表達(dá)與偽時(shí)間之間的（可能是非線性的）關(guān)系進(jìn)行建模东跪。然后，我們使用 associationTest函數(shù)來(lái)鑒定與偽時(shí)間存在顯著關(guān)聯(lián)的基因鹰溜。

BiocManager::install("tradeSeq")
library(tradeSeq)

# fit negative binomial GAM
sim <- fitGAM(sim)

# test for dynamic expression
ATres <- associationTest(sim)

然后越庇，我們可以根據(jù)p值挑選出變化最顯著的基因，并通過(guò)熱圖可視化其在整個(gè)偽時(shí)間或發(fā)育過(guò)程中的表達(dá)情況奉狈。這里卤唉，我們使用top250基因進(jìn)行展示。

topgenes <- rownames(ATres[order(ATres$pvalue), ])[1:250]
pst.ord <- order(sim$slingPseudotime_1, na.last = NA)
heatdata <- assays(sim)$counts[topgenes, pst.ord]
heatclus <- sim$GMM[pst.ord]

heatmap(log1p(heatdata), Colv = NA,
        ColSideColors = brewer.pal(9,"Set1")[heatclus])

image.png

分步運(yùn)行slingshot

接下來(lái)仁期，我們將重點(diǎn)介紹slingshot包的一些其他功能桑驱。我們使用該包內(nèi)置的slingshotExample數(shù)據(jù)集進(jìn)行說(shuō)明。此數(shù)據(jù)集中包含了細(xì)胞降維空間的低維表示跛蛋，以及k-means聚類(lèi)后的細(xì)胞簇標(biāo)簽熬的。

確定細(xì)胞分化譜系的全局結(jié)構(gòu)

首先，我們getLineages函數(shù)構(gòu)建細(xì)胞分化譜系的全局結(jié)構(gòu)赊级。該函數(shù)使用n × p的細(xì)胞降維空間矩陣和聚類(lèi)分群的標(biāo)簽向量n作為輸入押框，通過(guò)最小生成樹(shù)（MST）算法映射相鄰細(xì)胞群集之間的連接，并進(jìn)一步通過(guò)這些連接關(guān)系鑒定出細(xì)胞分化譜系的路徑理逊。該函數(shù)的輸出結(jié)果是一個(gè) SlingshotDataSet對(duì)象橡伞，其中包含了輸入文件盒揉，以及推斷出的MST（由鄰接矩陣表示）和分化譜系（細(xì)胞聚類(lèi)簇的有序向量）。

通過(guò)指定已知的起始簇和終點(diǎn)簇兑徘，我們還可以以完全無(wú)監(jiān)督的方式或以半監(jiān)督的方式執(zhí)行此分析刚盈。如果我們不指定起點(diǎn)，則slingshot根據(jù)簡(jiǎn)約性選擇一個(gè)起點(diǎn)挂脑，以使分割前譜系之間共享的細(xì)胞簇?cái)?shù)最大化藕漱。如果沒(méi)有分割或多個(gè)簇產(chǎn)生相同的簡(jiǎn)約分?jǐn)?shù)，則可以任意選擇起始簇崭闲。

對(duì)于我們的模擬數(shù)據(jù)肋联，這里選擇以“Cluster1”作為起始簇。但是刁俭，我們通常建議根據(jù)先驗(yàn)知識(shí)（樣本采集的時(shí)間或已知的基因標(biāo)記）確定初始簇的類(lèi)別牺蹄。

lin1 <- getLineages(rd, cl, start.clus = '1')
## Using full covariance matrix

lin1
## class: SlingshotDataSet 
## 
##  Samples Dimensions
##      140          2
## 
## lineages: 2 
## Lineage1: 1  2  3  5  
## Lineage2: 1  2  3  4  
## 
## curves: 0
plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16)
lines(lin1, lwd = 3, col = 'black')

image.png

在此步驟中，我們還可以指定已知的終點(diǎn)簇薄翅。當(dāng)構(gòu)造MST時(shí)，被指定為終端狀態(tài)的細(xì)胞簇將被約束為僅具有一個(gè)連接（即氓奈，它們必須是葉節(jié)點(diǎn)）翘魄。該約束可能會(huì)影響樹(shù)的其他部分的繪制方式，如以下示例中所示舀奶，在該示例中暑竟，我們將Cluster3指定為終點(diǎn)簇。

lin2 <- getLineages(rd, cl, start.clus= '1', end.clus = '3')
## Using full covariance matrix

plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16)
lines(lin2, lwd = 3, col = 'black', show.constraints = TRUE)

image.png

這種類(lèi)型的監(jiān)督對(duì)于確保結(jié)果與先前的生物學(xué)知識(shí)一致可能是有用的育勺。具體地但荤，它可以防止將已知的終末細(xì)胞命運(yùn)歸為暫時(shí)狀態(tài)。

我們還可以傳遞一些其他參數(shù)涧至，以實(shí)現(xiàn)對(duì)getLineages函數(shù)更好的控制：

• dist.fun is a function for computing distances between clusters. The default is squared distance between cluster centers normalized by their joint covariance matrix.
• omega is a granularity parameter, allowing the user to set an upper limit on connection distances. This can be useful for identifying outlier clusters which are not a part of any lineage or for separating distinct trajectories. This is typically a numeric argument, but a value of TRUE will indicate that a heuristic of 3 * (median edge length of unsupervised MST) should be used.

構(gòu)建平滑曲線確定細(xì)胞分化順序

為了構(gòu)建平滑的譜系曲線腹躁，我們使用getCurves函數(shù)對(duì)getLineages的結(jié)果進(jìn)行處理。該函數(shù)遵循類(lèi)似于（Hastie and Stuetzle 1989）中提出的主要曲線的迭代過(guò)程南蓬。如果只有一個(gè)分化譜系纺非，則生成的曲線只是經(jīng)過(guò)數(shù)據(jù)中心的主曲線，并進(jìn)行了一次調(diào)整：初始曲線是通過(guò)聚類(lèi)中心之間的線性連接而不是數(shù)據(jù)的第一個(gè)主成分而構(gòu)建的赘方。這種調(diào)整增加了穩(wěn)定性烧颖，并通常加快了算法的收斂速度。

當(dāng)有兩個(gè)或多個(gè)分化譜系時(shí)窄陡，我們向該算法添加一個(gè)額外的步驟：對(duì)共享細(xì)胞附近的曲線求平均炕淮。兩種譜系應(yīng)該在尚未分化的細(xì)胞上都很好地吻合，因此在每次迭代中跳夭，我們均對(duì)這些細(xì)胞附近的曲線進(jìn)行平均涂圆。這增加了算法的穩(wěn)定性们镜，并產(chǎn)生了平滑的分支譜系。

crv1 <- getCurves(lin1)
crv1
## class: SlingshotDataSet 
## 
##  Samples Dimensions
##      140          2
## 
## lineages: 2 
## Lineage1: 1  2  3  5  
## Lineage2: 1  2  3  4  
## 
## curves: 2 
## Curve1: Length: 15.045   Samples: 100.6
## Curve2: Length: 15.126   Samples: 103.5

plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16)
lines(crv1, lwd = 3, col = 'black')

image.png

getCurves函數(shù)的結(jié)果也是一個(gè)slingshotDataSet對(duì)象乘综，其中包含了分化譜系的主曲線和有關(guān)其擬合方式的其他信息憎账。我們可以使用pseudotime函數(shù)提取不同分化譜系中細(xì)胞的偽時(shí)間值，其中的NA表示未分配給特定譜系的細(xì)胞的值卡辰。使用curveWeights函數(shù)可以提取相似的權(quán)重矩陣胞皱，將每個(gè)細(xì)胞分配給一個(gè)或多個(gè)譜系。

可以使用curves函數(shù)訪問(wèn)曲線對(duì)象九妈，該函數(shù)將返回principal_curve對(duì)象的列表反砌。這些對(duì)象中包含以下信息：

• s: the matrix of points that make up the curve. These correspond to the orthogonal projections of the data points.
• ord: indices which can be used to put the cells along a curve in order based their projections.
  lambda: arclengths along the curve from the beginning to each cell’s projection. A matrix of these values is returned by the slingPseudotime function.
• dist_ind: the squared distances between data points and their projections onto the curve.
• dist: the sum of the squared projection distances.
• w: the vector of weights along this lineage. Cells that were unambiguously assigned to this lineage will have a weight of 1, while cells assigned to other lineages will have a weight of 0. It is possible for cells to have weights of 1 (or very close to 1) for multiple lineages, if they are positioned before a branching event. A matrix of these values is returned by the slingCurveWeights function.

在大型數(shù)據(jù)集上運(yùn)行slingshot

對(duì)于大型數(shù)據(jù)集，我們建議在運(yùn)行slingshot函數(shù)時(shí)使用approx_points參數(shù)萌朱。通過(guò)設(shè)置approx_points宴树，用戶可以指定曲線的分辨率（唯一點(diǎn)的數(shù)量），這可以大大縮短計(jì)算時(shí)間晶疼，并且精度損失最小酒贬。我們建議該值設(shè)置為100-200。

sim5 <- slingshot(sim, clusterLabels = 'GMM', reducedDim = 'PCA',
                   approx_points = 5)
## Using full covariance matrix

colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plotcol <- colors[cut(sim5$slingPseudotime_1, breaks=100)]

plot(reducedDims(sim5)$PCA, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(sim5), lwd=2, col='black')

image.png

sessionInfo()
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] mclust_5.4.6                uwot_0.1.8                 
##  [3] Matrix_1.2-18               SingleCellExperiment_1.12.0
##  [5] SummarizedExperiment_1.20.0 Biobase_2.50.0             
##  [7] GenomicRanges_1.42.0        GenomeInfoDb_1.26.0        
##  [9] IRanges_2.24.0              S4Vectors_0.28.0           
## [11] BiocGenerics_0.36.0         MatrixGenerics_1.2.0       
## [13] matrixStats_0.57.0          mgcv_1.8-33                
## [15] nlme_3.1-150                RColorBrewer_1.1-2         
## [17] slingshot_1.8.0             princurve_2.1.5            
## [19] BiocStyle_2.18.0           
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.5             RSpectra_0.16-0        compiler_4.0.3        
##  [4] BiocManager_1.30.10    XVector_0.30.0         bitops_1.0-6          
##  [7] tools_4.0.3            zlibbioc_1.36.0        digest_0.6.27         
## [10] evaluate_0.14          lattice_0.20-41        pkgconfig_2.0.3       
## [13] rlang_0.4.8            igraph_1.2.6           DelayedArray_0.16.0   
## [16] magick_2.5.0           yaml_2.2.1             xfun_0.18             
## [19] GenomeInfoDbData_1.2.4 stringr_1.4.0          knitr_1.30            
## [22] grid_4.0.3             rmarkdown_2.5          bookdown_0.21         
## [25] magrittr_1.5           codetools_0.2-16       splines_4.0.3         
## [28] htmltools_0.5.0        ape_5.4-1              stringi_1.5.3         
## [31] RCurl_1.98-1.2         FNN_1.1.3

參考來(lái)源：https://bioconductor.org/packages/release/bioc/vignettes/slingshot/inst/doc/vignette.html