jamesbang

V1

2022/10/16阅读:20主题:雁栖湖

🤩 scRNA-seq | 吐血整理的单细胞入门教程(差异分析)(十三)

1写在前面

完成了聚类后,我们就要进行差异分析,寻找差异基因了。🥳
由于scRNAseq高维数据,而且并没有明确的,你可以选择之前介绍的SC3包等,先进行聚类,然后确定了后,进行比较,或者采用生物学分组进行比较。😘

本期我们介绍一下常用的一些差异分析方法,再比较各种方法的准确性。🤒

2用到的包

rm(list = ls())
library(scRNA.seq.funcs)
library(edgeR)
#library(monocle)
library(MAST)
library(ROCR)

3示例数据

1️⃣ 由于数据量比较大,这里我们只比较一下NA19101NA19239的差异基因。🤨

molecules <- read.table("./molecules.txt", sep = "\t")
anno <- read.table("./annotation.txt", sep = "\t", header = TRUE)
keep <- anno[,1] == "NA19101" | anno[,1] == "NA19239"
data <- molecules[,keep]
group <- anno[keep,1]
batch <- anno[keep,4]
# 过滤一下
gkeep <- rowSums(data > 0) > 5;
counts <- data[gkeep,]
# Library size normalization
lib_size = colSums(counts)
norm <- t(t(counts)/lib_size * median(lib_size))
# Variant of CPM for datasets with library sizes of fewer than 1 mil molecules
norm
norm

2️⃣ 这里我们再提供一个事先计算好的差异基因数据,方便后面对不同方法的评估,TPs为真实的差异基因TNs为真实的无差异基因

DE <- read.table("./TPs.txt")
notDE <- read.table("./TNs.txt")

GroundTruth <- list(
DE = as.character(unlist(DE)),
notDE = as.character(unlist(notDE))
)

4Kolmogorov-Smirnov Test

4.1 差异分析

我们先介绍一种非参数检验的方法,即Kolmogorov-Smirnov test (KS-test)。😉

pVals <- apply(
norm, 1, function(x) {
ks.test(
x[group == "NA19101"],
x[group == "NA19239"]
)$p.value
}
)
# 校正
pVals <- p.adjust(pVals, method = "fdr")

4.2 准确性评估

这里我们使用KS-test的方法得到了5095个差异基因。🤯

sigDE <- names(pVals)[pVals < 0.05]
length(sigDE)
5095
5095

接着我们可以看下有多少个真的差异基因在这个KS-test里,也就是真阳性792个。🫠

sum(GroundTruth$DE %in% sigDE) 
792
792

我们再看一下无差异基因,也就是假阳性3190个。😤

sum(GroundTruth$notDE %in% sigDE)
3190
3190

4.3 准确性可视化

我们用利用ROC曲线进行一下可视化。

我们先自定义一个函数,称为DE_Quality_AUC,后面就不用总是计算了。😉

DE_Quality_AUC <- function(pVals) {
pVals <- pVals[names(pVals) %in% GroundTruth$DE |
names(pVals) %in% GroundTruth$notDE]
truth <- rep(1, times = length(pVals));
truth[names(pVals) %in% GroundTruth$DE] = 0;
pred <- ROCR::prediction(pVals, truth)
perf <- ROCR::performance(pred, "tpr", "fpr")
ROCR::plot(perf)
aucObj <- ROCR::performance(pred, "auc")
return(aucObj@y.values[[1]])
}

DE_Quality_AUC(pVals)
0.7954796
0.7954796

AUC = 0.7954796

5Wilcox/Mann-Whitney-U Test

Wilcox-rank-sum test也是一种非参数检验,看一下它的表现吧。

pVals <- apply(
norm, 1, function(x) {
wilcox.test(
x[group == "NA19101"],
x[group == "NA19239"]
)$p.value
}
)
# multiple testing correction
pVals <- p.adjust(pVals, method = "fdr")
DE_Quality_AUC(pVals)
0.8320326
0.8320326

AUC = 0.8320326,表现要比KS-test好一些。

6edgeR

edgeR基于广义线性模型(generalized linear model, glm)进行差异基因分析,而且还可以在模型中纳入多种影响因素,比如batch等等。

dge <- DGEList(
counts = counts,
norm.factors = rep(1, length(counts[1,])),
group = group
)
group_edgeR <- factor(group)
design <- model.matrix(~ group_edgeR)
dge <- estimateDisp(dge, design = design, trend.method = "none")
fit <- glmFit(dge, design)
res <- glmLRT(fit)
pVals <- res$table[,4]
names(pVals) <- rownames(res$table)

pVals <- p.adjust(pVals, method = "fdr")
DE_Quality_AUC(pVals)
0.8466764
0.8466764

AUC = 0.8466764,表现更好一些了。🤩

7MAST

log_counts <- log(counts + 1) / log(2)
fData <- data.frame(names = rownames(log_counts))
rownames(fData) <- rownames(log_counts);
cData <- data.frame(cond = group)
rownames(cData) <- colnames(log_counts)

obj <- FromMatrix(as.matrix(log_counts), cData, fData)
colData(obj)$cngeneson <- scale(colSums(assay(obj) > 0))
cond <- factor(colData(obj)$cond)

# Model expression as function of condition & number of detected genes
zlmCond <- zlm(~ cond + cngeneson, obj)

summaryCond <- summary(zlmCond, doLRT = "condNA19239")
summaryDt <- summaryCond$datatable

summaryDt <- as.data.frame(summaryDt)
pVals <- unlist(summaryDt[summaryDt$component == "H",4]) # H = hurdle model
names(pVals) <- unlist(summaryDt[summaryDt$component == "H",1])
pVals <- p.adjust(pVals, method = "fdr")
DE_Quality_AUC(pVals)
0.8284046
0.8284046

AUC = 0.8284046

8其他方法

这里补充一下其他方法,用的比较少,主要因为太慢了,单次运行一般都在1h以上,这里只提供一下代码,有愿意试一下的小伙伴可以试一试。

8.1 BPSC

library(BPSC)
bpsc_data <- norm[,batch=="NA19101.r1" | batch=="NA19239.r1"]
bpsc_group = group[batch=="NA19101.r1" | batch=="NA19239.r1"]

control_cells <- which(bpsc_group == "NA19101")
design <- model.matrix(~bpsc_group)
coef=2 # group label
res=BPglm(data=bpsc_data, controlIds=control_cells, design=design, coef=coef,
estIntPar=FALSE, useParallel = FALSE)
pVals = res$PVAL
pVals <- p.adjust(pVals, method = "fdr")
DE_Quality_AUC(pVals)

8.2 SCDE

library(scde)
cnts <- apply(
counts,
2,
function(x) {
storage.mode(x) <- 'integer'
return(x)
}
)
names(group) <- 1:length(group)
colnames(cnts) <- 1:length(group)
o.ifm <- scde::scde.error.models(
counts = cnts,
groups = group,
n.cores = 1,
threshold.segmentation = TRUE,
save.crossfit.plots = FALSE,
save.model.plots = FALSE,
verbose = 0,
min.size.entries = 2
)
priors <- scde::scde.expression.prior(
models = o.ifm,
counts = cnts,
length.out = 400,
show.plot = FALSE
)
resSCDE <- scde::scde.expression.difference(
o.ifm,
cnts,
priors,
groups = group,
n.randomizations = 100,
n.cores = 1,
verbose = 0
)
# Convert Z-scores into 2-tailed p-values
pVals <- pnorm(abs(resSCDE$cZ), lower.tail = FALSE) * 2
DE_Quality_AUC(pVals)

蛋糕
最后祝大家早日不卷!~

需要示例数据的小伙伴,在公众号回复scRNAseq3获取吧!

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰

分类:

后端

标签:

后端

作者介绍

jamesbang
V1

wx🔍: Grassssss 卷起来了