哎呀，终究还是绕不过差异表达分析。
单细胞测序，测得还是转录组数据，和bluk data的差别就是单细胞测序存在许多零值。
什么是零值？就是为零的值！（哈哈~废话啦）
就是有些很多基因没有表达，原因有二：1）细胞本身不表达某些基因；2）测序的过程当中没测到。
为此，利用传统的转录组数据分析软件（DEseq\edgeR）就不行。
就要用些专门针对单细胞转录组的分析软件。
这里先介绍SCDE，下一篇再介绍DEsingle。
1、安装
这个包真的是很迷。
官网的protocol写的完全不明白什么意思。
因为用的人很少，所以还是自行查阅。
安装方法：
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("scde")
然后，可以借鉴官网上提供的方法，进行升级。
require(devtools)
devtools::install_version('flexmix', '2.3-13')
devtools::install_github('hms-dbmi/scde', build_vignettes = FALSE)
升级完，貌似也没啥变化。
2、使用
还是按照官网的数据，看一下安装的有没有问题。
# 载入示例数据
data(es.mef.small)
# 利用行名MEF和ESC提取细胞
sg <- factor(gsub("(MEF|ESC).*", "\\1", colnames(es.mef.small)), levels = c("ESC", "MEF"))
# 换一个名字，你可以View看一下数据到底是什么样子。
names(sg) <- colnames(es.mef.small)
table(sg)
#结果展示
## sg
## ESC MEF
## 20 20
# clean up the dataset：对筛选得到的数据进行筛选，筛选条件如下。可以dim看一下数据集的大小。
cd <- clean.counts(es.mef.small, min.lib.size=1000, min.reads = 1, min.detected = 1)
接着就是构建和计算模型：
# calculate models
o.ifm <- scde.error.models(counts = cd, groups = sg, n.cores = 1, threshold.segmentation = TRUE, save.crossfit.plots = FALSE, save.model.plots = FALSE, verbose = 1)
在运行过程中这里会报错。
报错信息如下：
error info: Error in FUN(X[[i]], ...) :
trying to get slot "logLik" from an object of a basic class ("function") with no slots
第一次看到这个错误的时候，就一脸懵逼。
怎么办？这个错误看不懂呀。
老路子，网上爬资料，找解决办法。
奈何，无解。
那就先看看这个模型到底有没有吧。
data(o.ifm)
head(o.ifm)
事实证明，勇敢的往下走，还是很靠谱的。
这个错误并不影响模型本身的计算和结果，那这个错误有可能是这个包本身的bug。
暂时先不管这个报错了，继续。
# filter out cells that don't show positive correlation with
# the expected expression magnitudes (very poor fits)
valid.cells <- o.ifm$corr.a > 0
table(valid.cells)
#输出结果
## valid.cells
## TRUE
## 40
o.ifm <- o.ifm[valid.cells, ]
# estimate gene expression prior
o.prior <- scde.expression.prior(models = o.ifm, counts = cd, length.out = 400, show.plot = FALSE)
# define two groups of cells
groups <- factor(gsub("(MEF|ESC).*", "\\1", rownames(o.ifm)), levels = c("ESC", "MEF"))
names(groups) <- row.names(o.ifm)
# run differential expression tests on all genes.
ediff <- scde.expression.difference(o.ifm, cd, o.prior, groups = groups, n.randomizations = 100, n.cores = 1, verbose = 1)
head(ediff[order(ediff$Z, decreasing = TRUE), ])
write.table(ediff[order(abs(ediff$Z), decreasing = TRUE), ], file = "results.txt", row.names = TRUE, col.names = TRUE, sep = "\t", quote = FALSE)
........
下面的具体操作，就不一一展示了，在官网上很详细了。
官网链接：
http://hms-dbmi.github.io/scde/diffexp.html
以上。
网友评论