3.12 细胞轨迹推断
刘小泽写于2020.7.16
最后更新于
刘小泽写于2020.7.16
最后更新于
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:4,1:4]
c1 c2 c3 c4
G1 0 1 2 0
G2 2 3 2 5
G3 3 8 5 4
G4 5 16 6 8
> dim(counts)
[1] 750 300
# 构建一个sce对象
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):
altExpNames(0):geneFilter <- apply(assays(sim)$counts,1,function(x){
sum(x >= 3) >= 10
})
sim <- sim[geneFilter, ]
# 过滤掉11个基因
> dim(sim)
[1] 739 300FQnorm <- 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)pca <- prcomp(t(log1p(assays(sim)$norm)), scale. = FALSE)
rd1 <- pca$x[,1:2]
plot(rd1, col = rgb(0,0,0,.5), pch=16, asp = 1)library(destiny, quietly = TRUE)
dm <- DiffusionMap(t(log1p(assays(sim)$norm)))
rd2 <- cbind(DC1 = dm$DC1, DC2 = dm$DC2)
plot(rd2, col = rgb(0,0,0,.5), pch=16, asp = 1)reducedDims(sim) <- SimpleList(PCA = rd1, DiffMap = rd2)library(mclust, quietly = TRUE)
#根据PCA结果
cl1 <- Mclust(rd1)$classification
colData(sim)$GMM <- cl1
library(RColorBrewer)
plot(rd1, col = brewer.pal(9,"Set1")[cl1], pch=16, asp = 1)cl2 <- kmeans(rd1, centers = 4)$cluster
colData(sim)$kmeans <- cl2
plot(rd1, col = brewer.pal(9,"Set1")[cl2], pch=16, asp = 1)sim <- slingshot(sim, clusterLabels = 'GMM', reducedDim = 'PCA')
summary(sim$slingPseudotime_1)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 8.633 21.118 21.415 34.367 43.186colors <- 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')library(TSCAN)
data(lpsdata)
procdata <- preprocess(lpsdata)lpsmclust <- exprmclust(procdata)
# 然后看下结果
plotmclust(lpsmclust)lpsorder <- TSCANorder(lpsmclust)diffval <- difftest(procdata,lpsorder)head(row.names(diffval)[diffval$qval < 0.05])# 以STAT2基因为例
STAT2expr <- log2(lpsdata["STAT2",]+1)
singlegeneplot(STAT2expr, TSCANorder(lpsmclust,flip=TRUE,orderonly=FALSE))ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)# 默认使用DDRTree的方法
cds <- reduceDimension(cds, max_components = 2,
method = 'DDRTree')cds <- orderCells(cds)plot_cell_trajectory(cds, color_by = "Biological_Condition")plot_genes_in_pseudotime(cds[cg,],
color_by = "Biological_Condition")