单细胞分析绕不过去的坎，三种方法帮你解决！总有一款适合你！

单细胞数据的差异分析

大家好，我是风风。欢迎来到风风的从零开始单细胞系列。我们接着我们未完成的系列内容，算一算还蛮多内容可以扩展。今天我们来看看单细胞数据表达差异的几种常见方法。

咦？单细胞也有表达差异分析吗？当然有！表达差异并不是bulk seq数据的专属，而且用的方法也大同小异。2019年BMC Bioinformatics曾发过一篇比较单细胞RNA测序数据差异基因表达分析工具的文章，题目为“Comparative analysis of differential gene expression analysis tools for single-cell RNA sequencing data”，文章比较了11种不同的分析方法的优劣，并将不同方法进行了总结归纳，不仅提供了每种方法的应用链接，连相应文献也在这篇文章的引文中找到，真乃一文走完所有的路，让别人无路可走啊！我摘录了几张比较重要的结果如下：

我把文章放在了资料下载区，感兴趣的话可以后台回复关键词后获取。

今天我们主要介绍三种常见单细胞数据差异分析方法：Wilcox Test、edgeRh和Monocle，并仿照这篇文章，使用ROC曲线对同一份数据用这三种方法的差异分析结果进行比较。

数据导入与处理

首先我们导入相应的数据进行分析：

rm(list = ls())options(stringsAsFactors = FALSE)set.seed(123456) # 设置随机种子，方便重复结果# 安装R包# install.packages("devtools")# BiocManager::install("scRNA.seq.funcs")# BiocManager::install("edgeR")# BiocManager::install("MAST")# BiocManager::install("monocle")# install.packages("ROCR")library(scRNA.seq.funcs)library(edgeR)## Loading required package: limmalibrary(monocle)library(ROCR)# NA19101和NA19239的单细胞数据进行分析DE <- read.table("ENSG2.txt")notDE <- read.table("ENSG1.txt")GroundTruth <- list( DE = as.character(unlist(DE)), notDE = as.character(unlist(notDE)))molecules <- read.table("matrix.txt", sep = "\t") # 加载矩阵molecules[1:6,1:6]## NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03 NA19098.r1.A04## ENSG00000237683 0 0 0 1## ENSG00000187634 0 0 0 0## ENSG00000188976 3 6 1 3## ENSG00000187961 0 0 0 0## ENSG00000187583 0 0 0 0## ENSG00000187642 0 0 0 0## NA19098.r1.A05 NA19098.r1.A06## ENSG00000237683 0 0## ENSG00000187634 0 0## ENSG00000188976 4 2## ENSG00000187961 0 0## ENSG00000187583 0 0## ENSG00000187642 0 0anno <- read.table("annotation.txt", sep = "\t", header = TRUE) # 加载注释信息head(anno)## individual replicate well batch sample_id## 1 NA19098 r1 A01 NA19098.r1 NA19098.r1.A01## 2 NA19098 r1 A02 NA19098.r1 NA19098.r1.A02## 3 NA19098 r1 A03 NA19098.r1 NA19098.r1.A03## 4 NA19098 r1 A04 NA19098.r1 NA19098.r1.A04## 5 NA19098 r1 A05 NA19098.r1 NA19098.r1.A05## 6 NA19098 r1 A06 NA19098.r1 NA19098.r1.A06keep <- anno[,1] == "NA19101" | anno[,1] == "NA19239" data <- molecules[,keep]group <- anno[keep,1]batch <- anno[keep,4]# 去除至少在6个细胞中没有表达的基因gkeep <- rowSums(data > 0) > 5;counts <- data[gkeep,]# 标准化lib_size = colSums(counts)norm <- t(t(counts)/lib_size * median(lib_size))

接下来我们来比较各种单细胞差异分析的方法。

Wilcox Test

Wilcox检验，又称为Wilcox-rank-sum检验，这是一种大家非常熟悉的非参数检验方法，一般用于检验一个组的数值是否大于或者小于另一个组的数值。因此，这种方法通常被认为是检验两组之间表达量中位数的差异：

# 定义函数进行ROC分析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([email protected][[1]])}

# ROC曲线可以用来评估阈值的显著性效果，帮助我们调整阈值，可以计算出ROC曲线下的面积也就是AUC值进行准确性统计

# 接下来进行Wilcox Test差异分析pVals <- apply( norm, 1, function(x){ wilcox.test( x[group == "NA19101"], x[group == "NA19239"] )$p.value })# 多重检验矫正获取FDR值pVals <- p.adjust(pVals, method = "fdr")# 绘制ROC曲线DE_Quality_AUC(pVals) # 0.8320326

## [1] 0.8320326

edgeR

edgeR是基于基因表达的负二项式模型，并使用广义线性模型（GLM）框架，这得其他因素（如批次）可以纳入模型中，edgeR也是常见分析二代测序数据差异表达的方法：

dge <- DGEList( counts = counts, norm.factors = rep(1, length(counts[1,])), group = group)group_edgeR <- factor(group) # 设置组别design <- model.matrix(~ group_edgeR)# edgeR差异分析dge <- estimateDisp(dge, design = design, trend.method = "none")fit <- glmFit(dge, design) # 广义线性模型res <- glmLRT(fit)pVals <- res$table[,4] # 获取p值names(pVals) <- rownames(res$table)pVals <- p.adjust(pVals, method = "fdr") # 获取FDR值DE_Quality_AUC(pVals) # 0.8466764

## [1] 0.8466764

Monocle

Monocle是一种灵活的差异分析方法，可以选择使用几种不同的模型。对于计数数据，Monocle推荐使用负二项式模型（negbinomial size）。对于归一化数据，它建议将其进行对数转换，然后使用正态分布法（gaussianff）。与edgeR类似，这个方法使用GLM框架，所以理论上可以考虑到批次效应，但是实际上对于含有批次效应的数据一般不使用这种方法：

pd <- data.frame(group = group, batch = batch)rownames(pd) <- colnames(counts)pd <- new("AnnotatedDataFrame", data = pd)Obj <- newCellDataSet( as.matrix(counts), phenoData = pd, expressionFamily = negbinomial.size())## Warning in newCellDataSet(as.matrix(counts), phenoData = pd, expressionFamily =## negbinomial.size()): Warning: featureData must contain a column verbatim named## 'gene_short_name' for certain functions## Warning in newCellDataSet(as.matrix(counts), phenoData = pd, expressionFamily =## negbinomial.size()): Warning: featureData must contain a column verbatim named## 'gene_short_name' for certain functions## Warning in newCellDataSet(as.matrix(counts), phenoData = pd, expressionFamily =## negbinomial.size()): Warning: featureData must contain a column verbatim named## 'gene_short_name' for certain functionsObj## CellDataSet (storageMode: environment)## assayData: 16026 features, 576 samples ## element names: exprs ## protocolData: none## phenoData## sampleNames: NA19101.r1.A01 NA19101.r1.A02 ... NA19239.r3.H12 (576## total)## varLabels: group batch Size_Factor## varMetadata: labelDescription## featureData: none## experimentData: use 'experimentData(object)'## Annotation:Obj <- estimateSizeFactors(Obj)Obj <- estimateDispersions(Obj)## Warning: `group_by_()` was deprecated in dplyr 0.7.0.## Please use `group_by()` instead.## See vignette('programming') for more help## Warning: `select_()` was deprecated in dplyr 0.7.0.## Please use `select()` instead.## Warning in log(ifelse(y == 0, 1, y/mu)): 产生了NaNs## Warning: step size truncated due to divergence## Removing 169 outliersres <- differentialGeneTest(Obj, fullModelFormulaStr = "~group")head(res)## status family pval qval## ENSG00000237683 OK negbinomial.size 0.470353620 0.59881531## ENSG00000187634 OK negbinomial.size 0.414632560 0.54771690## ENSG00000188976 OK negbinomial.size 0.003590441 0.01114692## ENSG00000187961 OK negbinomial.size 0.199242717 0.31316828## ENSG00000187583 OK negbinomial.size 0.517081439 0.64213468## ENSG00000187642 OK negbinomial.size 0.040460145 0.08679083pVals <- res[,3]names(pVals) <- rownames(res)pVals <- p.adjust(pVals, method = "fdr")DE_Quality_AUC(pVals) # 0.8252662