比较不同的对单细胞转录组数据寻找差异基因的方法

背景介绍

如果是bulk RNA-seq,那么现在最流行的就是DESeq2 和 edgeR啦,而且有很多经过了RT-qPCR 验证过的真实测序数据可以来评价不同的差异基因算法的表现。

对单细胞测序数据来说,通常需要先聚类之后把细胞群体进行分组,然后来比较不同的组的差异表达情况。当然,也有不少单细胞测序实验设计本身就有时间点,不同个体来源,不同培养条件这样的分组!

同时还有不少方法是不需要预先分类的,因为分类本身就会引入偏差。

跟bulk RNA-seq不一样的地方是,scRNA-seq通常涉及到的样本数量更多。这时候可以使用非参检验算法,比如Kolmogorov-Smirnov test (KS-test) 等等。

下面用一个测试数据来评价一下不同的算法的表现。处理同样的表达矩阵得到差异结果跟已知的差异结果进行比较看看overlap怎么样。评价指标主要是:

  • 召回率,True positive rate (TPR), TP/(TP + FN)
  • 准确率,False positive rate (FPR), FP/(FP+TP)
  • receiver-operating-characteristic curve (ROC)
  • area under this curve (AUC)

所以需要安装并且加载一些包,安装代码如下;

install.packages('ROCR')
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R")
biocLite("MAST") 
biocLite("scde") 
install.packages("devtools")
library("devtools")
install_github("BPSC","nghiavtr")
library("BPSC")

加载代码如下:

library(ROCR)
library(edgeR)
library(DESeq2)

library(scde)
library(BPSC)
library(MAST)
library(monocle) 

加载测试数据

这里选取的是芝加哥大学Yoav Gilad lab实验的Tung et al 2017的单细胞测序文章的数据

## 读取tung文章的数据,生成测试数据,这个代码不需要运行。
if(F){

  DE <- read.table("tung/TPs.txt")
  notDE <- read.table("tung/TNs.txt")
  GroundTruth <- list(DE=as.character(unlist(DE)), notDE=as.character(unlist(notDE)))

  molecules <- read.table("tung/molecules.txt", sep = "\t")
  anno <- read.table("tung/annotation.txt", sep = "\t", header = TRUE)
  keep <- anno[,1] == "NA19101" | anno[,1] == "NA19239"
  data <- molecules[,keep]
  group <- anno[keep,1]
  batch <- anno[keep,4]
  # remove genes that aren't expressed in at least 6 cells
  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
  rm(molecules)
  rm(data)
  save.image(file = 'scRNAseq_DEG_input.Rdata')
}
load(file = 'scRNAseq_DEG_input.Rdata')
# 我已经把测试数据保存为rdata数据格式了,直接加载。
dim(counts);
## [1] 16026   576
dim(norm);
## [1] 16026   576
dim(DE);
## [1] 1083    1
dim(notDE);
## [1] 10897     1
table(group)
## group
## NA19098 NA19101 NA19239 
##       0     288     288

可以看到这里需要选择的测试数据来源于2个人,每个人都有288个细胞的表达数据。

就是要对它们进行差异比较,而已知的1083个基因是确定显著差异的,另外10897个基因是确定不显著的。(首先,我们要假定这个是金标准!!!)

但是总共却是16026个基因,所以有一些基因是不确定显著与否的。

差异分析方法大全

Kolmogorov-Smirnov test

KS检验有两个弊端,首先是它假设基因表达量是连续的,如果有很多细胞表达量一致,比如都是0,表现就很差。其次它对大样本量太敏感了,可能其实差异并不大,但是样本数量很多,也会被认为是显著差异。

pVals <- apply(norm, 1, function(x) {
        ks.test(x[group =="NA19101"], 
                x[group=="NA19239"])$p.value
         })
# multiple testing correction
pVals <- p.adjust(pVals, method = "fdr")
sigDE <- names(pVals)[pVals < 0.05]
length(sigDE) 
## [1] 5095
# Number of KS-DE genes
sum(GroundTruth$DE %in% sigDE) 
## [1] 792
# Number of KS-DE genes that are true DE genes
sum(GroundTruth$notDE %in% sigDE)
## [1] 3190
tp <- sum(GroundTruth$DE %in% sigDE)
fp <- sum(GroundTruth$notDE %in% sigDE)
tn <- sum(GroundTruth$notDE %in% names(pVals)[pVals >= 0.05])
fn <- sum(GroundTruth$DE %in% names(pVals)[pVals >= 0.05])
tpr <- tp/(tp + fn)
fpr <- fp/(fp + tn)
cat(c(tpr, fpr))
## 0.7346939 0.2944706
ks_pVals=pVals

可以看到KS检验判断的显著差异基因实在是太多了,高达5095个。所以它能找回来792个真正的差异基因。但是却找到了3190个假阳性。所以计算得到召回率73.46%,但是准确率只有29.44%,这个表现不佳。

再看看ROC和RUC

# Only consider genes for which we know the ground truth
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")
aucObj@y.values[[1]] # AUC
## [1] 0.7954796

把这两个评价分析包装成函数,后面可以直接使用!

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_rate <- function(sigDE) {
  (length(sigDE) )
  # Number of KS-DE genes
  ( sum(GroundTruth$DE %in% sigDE) )
  # Number of KS-DE genes that are true DE genes
  (sum(GroundTruth$notDE %in% sigDE))
  tp <- sum(GroundTruth$DE %in% sigDE)
  fp <- sum(GroundTruth$notDE %in% sigDE)
  tn <- sum(GroundTruth$notDE %in% names(pVals)[pVals >= 0.05])
  fn <- sum(GroundTruth$DE %in% names(pVals)[pVals >= 0.05])
  tpr <- tp/(tp + fn)
  fpr <- fp/(fp + tn)
  cat(c(tpr, fpr))
}

Wilcox/Mann-Whitney-U Test

也是一种非参检验,通常比较两个组数据的median的差异。

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")
sigDE <- names(pVals)[pVals < 0.05]
Wilcox_pVals=pVals
DE_Quality_rate(sigDE)
## 0.8376623 0.3729346
DE_Quality_AUC(pVals) 
## [1] 0.8320326

召回率是81.9%,准确率是31.9%,这个表现不佳。

edgeR

edgeR包在bulk RNA-seq测序领域应用很广泛,基于负二项分布模型,应用了 generalized linear model (GLM) 算法

library(edgeR)
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")
sigDE <- names(pVals)[pVals < 0.05]
edgeR_pVals=pVals
DE_Quality_rate(sigDE)
## 0.8692022 0.3948121
DE_Quality_AUC(pVals) 
## [1] 0.8477189

召回率是86.9%,准确率是39.4%,表现越来越好了。

Monocle

monocle不仅仅针对于基于read counts的表达矩阵,还可以是已经被各种normalization的表达矩阵,比如基于RPKM/TPM等等,它会把被normalization的表达矩阵用normal/gaussian model (gaussianff()) 算法处理一下。差异分析的时候同样也是基于负二项分布模型,应用了 generalized linear model (GLM) 算法

library(monocle)
pd <- data.frame(group=group, batch=batch)
rownames(pd) <- colnames(counts)
pd <- new("AnnotatedDataFrame", data = pd)
## 针对于基于read counts的表达矩阵
Obj <- newCellDataSet(as.matrix(counts), phenoData=pd, 
        expressionFamily=negbinomial.size()) 
Obj <- estimateSizeFactors(Obj)
Obj <- estimateDispersions(Obj)
res <- differentialGeneTest(Obj,fullModelFormulaStr="~group")
pVals <- res[,3]
names(pVals) <- rownames(res)
pVals <- p.adjust(pVals, method = "fdr")
sigDE <- names(pVals)[pVals < 0.05]
monocle_pVals=pVals
DE_Quality_rate(sigDE)
DE_Quality_AUC(pVals)

monocle做差异分析的耗时非常夸张,召回率是84.2%,准确率是38.1%

MAST

MAST基于 zero-inflated negative binomial 分布模型

library(MAST)
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.SingleCellAssay(~cond + cngeneson, obj) 
summaryCond <- summary(zlmCond, doLRT="condNA19101")
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")
sigDE <- names(pVals)[pVals < 0.05]
MAST_pVals=pVals
DE_Quality_rate(sigDE)
DE_Quality_AUC(pVals)

召回率是82.8%,准确率是34.9.%

BPSC

这个用的是 Poisson-Beta 分布模型

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")
sigDE <- names(pVals)[pVals < 0.05]
BPSC_pVals=pVals
DE_Quality_rate(sigDE)
DE_Quality_AUC(pVals)

召回率是64.8%,准确率是30.7.%

SCDE

SCDE是第一个特意针对单细胞转录组测序数据的差异分析而设计的,用贝叶斯统计方法把表达矩阵拟合到 zero-inflated negative binomial 分布模型里面。

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
pVals <- p.adjust(pVals, method = "fdr")
sigDE <- names(pVals)[pVals < 0.05]
SCDE_pVals=pVals
DE_Quality_rate(sigDE)
DE_Quality_AUC(pVals)
save(SCDE_pVals,BPSC_pVals,MAST_pVals,monocle_pVals,edgeR_pVals,Wilcox_pVals,ks_pVals,file = 'DEG_results.Rdata')

统计学基础

负二项分布 negative binomial model

set.seed(1)
hist(rnbinom(1000, mu=10, size=100), col="grey50", xlab="Read Counts", main="Negative Binomial")

这个是被应用的最广泛的转录组表达数据分布模型。但是对单细胞转录组测序数据来说,因为有很高的dropout情况,导致模型失准,所以就提出来了zero-inflated negative binomial models

zero-inflated negative binomial models

d = 0.5;
counts <- rnbinom(1000, mu=10, size=100);
counts[runif(1000) < d] = 0;
hist(counts, col="grey50", xlab="Read Counts", main="Zero-inflated NB")

就是在原始的负二项分布数据里面随机挑选一些低表达量基因,给它们人为赋值为0表达量值。

Poisson-Beta distribution

a = 0.1
b = 0.1
g = 100
lambdas = rbeta(1000, a, b)
counts = sapply(g*lambdas, function(l) {rpois(1, lambda=l)})
hist(counts, col="grey50", xlab="Read Counts", main="Poisson-Beta")

文中提到的测试数据:

http://www.biotrainee.com/jmzeng/scRNA/DESeq_table.rds

http://www.biotrainee.com/jmzeng/scRNA/pollen.rds

http://www.biotrainee.com/jmzeng/scRNA/tung_umi.rds

http://www.biotrainee.com/jmzeng/scRNA/usoskin1.rds

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-01-21

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏ATYUN订阅号

【学术】强化学习系列(下):贝尔曼方程

在前一篇文章中,我们学习了马尔可夫决策和强化学习框架的一些主要组成部分。在本文中,我们将建立在这一理论上,学习价值函数和贝尔曼方程。 回报和返还(return...

4157
来自专栏梦里茶室

TensorFlow 深度学习笔记 Logistic Classification

Logistic Classification Github工程地址:https://github.com/ahangchen/GDLnotes 欢迎sta...

18610
来自专栏人工智能LeadAI

机器学习算法集锦

摘要: 机器学习 机器学习(Machine Learning, ML)是一门多领域交叉学科,涉及概率论、统计学、逼近论、凸分析、算法复杂度理论等多门学科。专门研...

3295
来自专栏机器学习算法原理与实践

用scikit-learn研究局部线性嵌入(LLE)

    在局部线性嵌入(LLE)原理总结中,我们对流形学习中的局部线性嵌入(LLE)算法做了原理总结。这里我们就对scikit-learn中流形学习的一些算法做...

832
来自专栏机器学习算法与Python学习

常用机器学习与数据挖掘相关术语(该充充电了...)

Sampling(采样): Simple Random Sampling(简单随机采样), OfflineSampling(离线等可能K...

2827
来自专栏大数据挖掘DT机器学习

数据处理的统计学习(scikit-learn教程)

Scikit-learn 是一个紧密结合Python科学计算库(Numpy、Scipy、matplotlib),集成经典机器学习算法的Python模块。 一、统...

4615
来自专栏社区的朋友们

Kaggle 实战:Ghouls, Goblins, and Ghosts

本例使用R语言中的决策树以及随机森林package对kaggle的一个分类问题解题的全部过程。本文需要读者对机器学习中的决策树、随机森林的原理有所了解,并且知道...

7860
来自专栏PPV课数据科学社区

数据挖掘常用模型构建示例(R语言版)

Linear Regression library(MASS) lm_fit = lm(medv~poly(rm,2)+crim,data = Boston) ...

2705
来自专栏ATYUN订阅号

正则化贪心森林(RGF)的入门简介,含案例研究

作为曾参与机器学习竞赛的数据科学家,我一直在寻找“非主流”的算法。使用这些算法可能不会成为竞争的赢家。但是他们有独特的预测方式。这些算法可以在集成模型中使用,以...

4186
来自专栏深度学习自然语言处理

word2vec理论与实践

导读 本文简单的介绍了Google 于 2013 年开源推出的一个用于获取 word vector 的工具包(word2vec),并且简单的介绍了其中的两个训练...

2696

扫码关注云+社区