前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >比较不同的对单细胞转录组数据寻找差异基因的方法

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

作者头像
生信技能树
发布2018-03-09 10:47:14
8.8K0
发布2018-03-09 10:47:14
举报
文章被收录于专栏:生信技能树

背景介绍

如果是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)

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

代码语言:javascript
复制
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")

加载代码如下:

代码语言:javascript
复制
library(ROCR)
library(edgeR)
library(DESeq2)

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

加载测试数据

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

代码语言:javascript
复制
## 读取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);
代码语言:javascript
复制
## [1] 16026   576
代码语言:javascript
复制
dim(norm);
代码语言:javascript
复制
## [1] 16026   576
代码语言:javascript
复制
dim(DE);
代码语言:javascript
复制
## [1] 1083    1
代码语言:javascript
复制
dim(notDE);
代码语言:javascript
复制
## [1] 10897     1
代码语言:javascript
复制
table(group)
代码语言:javascript
复制
## group
## NA19098 NA19101 NA19239 
##       0     288     288

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

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

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

差异分析方法大全

Kolmogorov-Smirnov test

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

代码语言:javascript
复制
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) 
代码语言:javascript
复制
## [1] 5095
代码语言:javascript
复制
# Number of KS-DE genes
sum(GroundTruth$DE %in% sigDE) 
代码语言:javascript
复制
## [1] 792
代码语言:javascript
复制
# Number of KS-DE genes that are true DE genes
sum(GroundTruth$notDE %in% sigDE)
代码语言:javascript
复制
## [1] 3190
代码语言:javascript
复制
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))
代码语言:javascript
复制
## 0.7346939 0.2944706
代码语言:javascript
复制
ks_pVals=pVals

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

再看看ROC和RUC

代码语言:javascript
复制
# 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)
代码语言:javascript
复制
aucObj <- ROCR::performance(pred, "auc")
aucObj@y.values[[1]] # AUC
代码语言:javascript
复制
## [1] 0.7954796

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

代码语言:javascript
复制
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的差异。

代码语言:javascript
复制
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)
代码语言:javascript
复制
## 0.8376623 0.3729346
代码语言:javascript
复制
DE_Quality_AUC(pVals) 
代码语言:javascript
复制
## [1] 0.8320326

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

edgeR

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

代码语言:javascript
复制
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)
代码语言:javascript
复制
## 0.8692022 0.3948121
代码语言:javascript
复制
DE_Quality_AUC(pVals) 
代码语言:javascript
复制
## [1] 0.8477189

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

Monocle

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

代码语言:javascript
复制
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 分布模型

代码语言:javascript
复制
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 分布模型

代码语言:javascript
复制
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 分布模型里面。

代码语言:javascript
复制
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)
代码语言:javascript
复制
save(SCDE_pVals,BPSC_pVals,MAST_pVals,monocle_pVals,edgeR_pVals,Wilcox_pVals,ks_pVals,file = 'DEG_results.Rdata')

统计学基础

负二项分布 negative binomial model

代码语言:javascript
复制
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

代码语言:javascript
复制
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

代码语言:javascript
复制
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

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2018-01-21,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 背景介绍
  • 加载测试数据
  • 差异分析方法大全
    • Kolmogorov-Smirnov test
      • Wilcox/Mann-Whitney-U Test
        • edgeR
          • Monocle
            • MAST
              • BPSC
                • SCDE
                • 统计学基础
                  • 负二项分布 negative binomial model
                    • zero-inflated negative binomial models
                      • Poisson-Beta distribution
                      领券
                      问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档