1.tradeSeq包
rm(list=ls())options(stringsAsFactors = F)library(Seurat)library(gplots)library(ggplot2)library(monocle)library(dplyr)load(file = 'output_of_monocle.Rdata')cds=my_cds_subsetcolnames(pData(cds))plot_cell_trajectory(cds)library(tradeSeq)library(RColorBrewer)library(SingleCellExperiment)library(monocle)#因为加载tradeSeq包后反而没有extract_monocle_info这个函数，因此需要我们source一下source('extract_monocle_info.R')info <- extract_monocle_info(cds)infosce <- fitGAM(counts = Biobase::exprs(cds),cellWeights = info$cellWeights,pseudotime = info$pseudotime)table(rowData(sce)$tradeSeq$converged)assoRes <- associationTest(sce)head(assoRes)# Discovering progenitor marker genesstartRes <- startVsEndTest(sce)oStart <- order(startRes$waldStat, decreasing = TRUE)sigGeneStart <- names(sce)[oStart[3]]counts= Biobase::exprs(cds)plotSmoothers(sce,counts,gene = sigGeneStart)# Comparing specific pseudotime values within a lineagecustomRes <- startVsEndTest(sce, pseudotimeValues = c(0.1, 0.8))endRes <- diffEndTest(sce)o <- order(endRes$waldStat, decreasing = TRUE)sigGene <- names(sce)[o[1]]plotSmoothers(sce, counts, sigGene)# Discovering genes with different expression patternspatternRes <- patternTest(sce)oPat <- order(patternRes$waldStat, decreasing = TRUE)head(rownames(patternRes)[oPat])plotSmoothers(sce,counts,gene = rownames(patternRes)[oPat][4])
2.弹弓函数slingshot
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(gplots)
library(ggplot2)
library(monocle)
library(dplyr)
load(file = 'output_of_monocle.Rdata')
cds=my_cds_subset
phe=pData(cds)
counts = Biobase::exprs(cds)
dim(counts) #[1] 1722 642
library(destiny)
destinyObj <- as.ExpressionSet(as.data.frame(t(counts)))
destinyObj$condition <- factor(phe$celltype)
sigma=50
# 类似于 PCA 分析
dm <- DiffusionMap(destinyObj,
sigma,
rotate = TRUE)
# 也可以看主成分
plot(eigenvalues(dm),
ylim = 0:1,
pch = 20,
xlab ='Diffusion component (DC)',
ylab ='Eigenvalue')
# 使用 Slingshot 来进行 pseudotime 分析，基于 diffusion map的主要成分
# 前面的图看的 DC3到DC4之间有一个拐点，所以我们考虑前面的4个diffusion map的主要成分
table(dm$condition)
library(slingshot)
crv <- slingshot(
dm@eigenvectors[,1:15],
dm$condition,
start.clus = 'FCGR3A+ Mono',
end.clus='CD14+ Mono',
maxit=100000,
# shrink.method="cosine" ,
shrink.method="tricube"
)
crv
sce=crv
sce
# 下面的代码要求 slingshot_2.0.0 以上版本
# https://bioconductor.org/packages/release/bioc/vignettes/slingshot/inst/doc/vignette.html
summary(sce$slingPseudotime_1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 8.631 21.121 21.414 34.363 43.185
library(grDevices)
colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plotcol <- colors[cut(sce$slingPseudotime_1, breaks=100)]
plot(reducedDims(sce)$PCA, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(sce), lwd=2, col='black')
3.tradeSeq和弹弓函数结合
#tradeSeq和弹弓函数结合
rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(gplots)
library(ggplot2)
library(monocle)
library(dplyr)
load(file = 'output_of_monocle.Rdata')
cds=my_cds_subset
phe=pData(cds)
counts = Biobase::exprs(cds)
sce <- SingleCellExperiment(assays = List(counts = counts))
sce
# 官网的示例数据，让你熟悉后面的要求
library(slingshot, quietly = FALSE)
data("slingshotExample")
rd <- slingshotExample$rd #rd就是降维后的两个主成分，可以是PCA,UMAP 等等，也可以是DiffusionMap class {destiny}
cl <- slingshotExample$cl
dim(rd) # data representing cells in a reduced dimensional space
head(rd)
table(cl)
length(cl) # vector of cluster labels
## ----genefilt-----------------------------------------------------------------
# filter genes down to potential cell-type markers
# at least M (15) reads in at least N (15) cells
geneFilter <- apply(assays(sce)$counts,1,function(x){
sum(x >= 3) >= 10
})
sce <- sce[geneFilter, ]
sce
## ----norm---------------------------------------------------------------------
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(sce)$norm <- FQnorm(assays(sce)$counts)
pca <- prcomp(t(log1p(assays(sce)$norm)), scale. = FALSE)
rd1 <- pca$x[,1:2]
plot(rd1,
col = rgb(0,0,0,.5),
pch=16,
asp = 1)
## ----umap, cache=TRUE---------------------------------------------------------
library(uwot)
rd2 <- uwot::umap(t(log1p(assays(sce)$norm)))
colnames(rd2) <- c('UMAP1', 'UMAP2')
plot(rd2,
col = rgb(0,0,0,.5),
pch=16,
asp = 1)
#可以看到pca和umap其实都是可以自己制作，并不一定要是seurat这样的一站式R包
## ----add_RDs, cache=TRUE------------------------------------------------------
reducedDims(sce) <- SimpleList(PCA = rd1, UMAP = rd2)
# 聚类也可以自己选择R 包
## ----clustering_mclust--------------------------------------------------------
library(mclust, quietly = TRUE)
cl1 <- Mclust(rd1)$classification
table(cl1)
colData(sce)$GMM <- cl1
library(RColorBrewer)
plot(rd1,
col = brewer.pal(9,"Set1")[cl1],
pch=16,
asp = 1)
cl2 <- kmeans(rd1, centers = 4)$cluster
table(cl2)
colData(sce)$kmeans <- cl2
plot(rd1,
col = brewer.pal(9,"Set1")[cl2],
pch=16,
asp = 1)
table(cl1,cl2)
#降维聚类分群细节，但是不同流程理论上是有可比性
#真正的函数在这里slingshot
sce <- slingshot(sce,
clusterLabels = 'GMM',
reducedDim = 'PCA')
summary(sce$slingPseudotime_1)
library(grDevices)
colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plotcol <- colors[cut(sce$slingPseudotime_1, breaks=100)]
plot(reducedDims(sce)$PCA,
col = plotcol,
pch=16,
asp = 1)
lines(SlingshotDataSet(sce),
lwd=2,
col='black')
## ----plot_curve_2-------------------------------------------------------------
plot(reducedDims(sce)$PCA,
col = brewer.pal(9,'Set1')[sce$GMM],
pch=16,
asp = 1)
lines(SlingshotDataSet(sce),
lwd=2,
type = 'lineages',
col = 'black')
## ----tradeseq, eval=FALSE-----------------------------------------------------
library(tradeSeq)
sce <- fitGAM(sce)
# test for dynamic expression
ATres <- associationTest(sce)
topgenes <- rownames(ATres[order(ATres$pvalue), ])[1:250];topgenes
pst.ord <- order(sce$slingPseudotime_1, na.last = NA)
heatdata <- assays(sce)$counts[topgenes, pst.ord]
heatclus <- sce$GMM[pst.ord]
heatmap(log1p(heatdata),
Colv = NA,
ColSideColors = brewer.pal(9,"Set1")[heatclus])
pst.ord <- order(sce$slingPseudotime_1, na.last = NA)
heatdata <- assays(sce)$counts[topgenes, pst.ord]
heatclus <- sce$GMM[pst.ord]
heatmap(log1p(heatdata),
Colv = NA,
ColSideColors = brewer.pal(9,"Set1")[heatclus])
rd=rd2
cl=cl2
lin1 <- getLineages(rd, cl, start.clus = '1')
lin1
plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16)
lines(SlingshotDataSet(lin1), lwd = 3, col = 'black')
## ----lines_sup_end------------------------------------------------------------
lin2 <- getLineages(rd, cl, start.clus= '1', end.clus = '3')
plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16)
lines(SlingshotDataSet(lin2), lwd = 3, col = 'black', show.constraints = TRUE)
## ----curves-------------------------------------------------------------------
crv1 <- getCurves(lin1)
crv1
plot(rd, col = brewer.pal(9,"Set1")[cl], asp = 1, pch = 16)
lines(SlingshotDataSet(crv1), lwd = 3, col = 'black')
## ----sling_approxpoints-------------------------------------------------------
sce5 <- slingshot(sce, clusterLabels = 'GMM', reducedDim = 'PCA',
approx_points = 5)
colors <- colorRampPalette(brewer.pal(11,'Spectral')[-6])(100)
plotcol <- colors[cut(sce5$slingPseudotime_1, breaks=100)]
plot(reducedDims(sce5)$PCA, col = plotcol, pch=16, asp = 1)
lines(SlingshotDataSet(sce5), lwd=2, col='black')
## ----sling_omega--------------------------------------------------------------
rd2 <- rbind(rd, cbind(rd[,2]-12, rd[,1]-6))
cl2 <- c(cl, cl + 10)
pto2 <- slingshot(rd2, cl2, omega = TRUE, start.clus = c(1,11))
plot(rd2, pch=16, asp = 1,
col = c(brewer.pal(9,"Set1"), brewer.pal(8,"Set2"))[cl2])
lines(SlingshotDataSet(pto2), type = 'l', lwd=2, col='black')
## ----sling_multtraj-----------------------------------------------------------
plot(rd2, pch=16, asp = 1,
col = c(brewer.pal(9,"Set1"), brewer.pal(8,"Set2"))[cl2])
lines(SlingshotDataSet(pto2), lwd=2, col='black')
## ----session------------------------------------------------------------------
sessionInfo()
若有收获，就点个赞吧
0 人点赞
