在评估不同软件性能的时候，我们会需要模拟一些数据。由于模拟数据当中的情况是已知的，例如差异表达基因的数目。因此，通过比较不同软件在模拟数据上的效果，我们可以获得软件的量化性能指标，例如灵敏度、特异性和准确度等。

本文根据 DESeq2 文章中的方法记录如何进行简单的基于负二项分布（Negative Binomial distribution）模拟RNA-seq基因表达数据。

为了模拟基因表达数据，我们需要：

1. 从已知的测序数据中获取基因表达的均值与离散程度作为负二项分布的参数
2. 构建样本信息，包括基因数目、样本数目、差异表达程度和样本分组等
3. 基于负二项分布模拟基因表达数据

Mean-Dispersion Estimation

首先，我们从2010 Nature发布的RNA-seq数据中，获取真实基因表达的均值(mean)与离散程度(dispersion)，以此作为我们模拟表达的基因的均值与离散程度。

DESeq2paper 提供了该数据的 RData (https://www.huber.embl.de/DESeq2paper/data/pickrell_sumexp.RData)

但是，在实际操作过程中，发现从 DESeq2paper 主页上下载的数据在读入后会出错

load("data/pickrell_sumexp.RData")
pickrell

# class: SummarizedExperiment 
# Error in h(simpleError(msg, call)) : 
#   error in evaluating the argument 'x' in selecting a method for function 'nrow': no slot of name "elementMetadata" for this object of class "SummarizedExperiment"
# Error during wrapup: no slot of name "NAMES" for this object of class "SummarizedExperiment"
# Error: no more error handlers available (recursive errors?); invoking 'abort' restart

搜索发现 DESeq2 作者 Michael Love 本人也遇到这个问题

https://support.bioconductor.org/p/131099/

但他也没有修复这个官网上的数据，而是把修复的数据放到：https://github.com/mikelove/apeglmPaper/blob/master/sim/rdata/pickrell_sumexp.RData

下载该数据，放到项目目录的 data/ 文件夹下。

接下来，使用 DESeq2 中的函数从数据中获取基因表达的均值与离散程度。

pickrell 数据中包含69个样本共56299个基因的表达数据

# mean-dispersion estimation from pickrell
library(DESeq2)
# https://github.com/mikelove/apeglmPaper/blob/master/sim/rdata/pickrell_sumexp.RData
load("data/pickrell_sumexp.RData")
pickrell

## class: RangedSummarizedExperiment 
## dim: 56299 69 
## metadata(0):
## assays(1): counts
## rownames(56299): ENSG00000000003 ENSG00000000005 ... ENSG00000261841
##   ENSG00000261842
## rowData names(0):
## colnames: NULL
## colData names(1): rownames

数据的均值和离散程度分别存储在 mcols(ddspickrell) 的 baseMean 和 dispGeneEst 列。

ddspickrell <- DESeqDataSet(pickrell, ~1)
ddspickrell <- estimateSizeFactors(ddspickrell)
# estimate gene-wise dispersion
ddspickrell <- estimateDispersionsGeneEst(ddspickrell)
# estimate fitted dispersion
ddspickrell <- estimateDispersionsFit(ddspickrell)
head(mcols(ddspickrell))

## DataFrame with 6 rows and 6 columns
##                   baseMean     baseVar   allZero dispGeneEst dispGeneIter
##                  <numeric>   <numeric> <logical>   <numeric>    <numeric>
## ENSG00000000003   0.992996 4.91552e+00     FALSE  1.92047047            8
## ENSG00000000005   0.018963 1.22328e-02     FALSE  0.00000001           28
## ENSG00000000419 128.895411 1.11440e+03     FALSE  0.05723383           10
## ENSG00000000457  79.020253 5.77887e+02     FALSE  0.04260266           11
## ENSG00000000460  57.365964 6.06033e+02     FALSE  0.11135132           14
## ENSG00000000938 713.350827 1.31691e+05     FALSE  0.25701274            8
##                   dispFit
##                 <numeric>
## ENSG00000000003  2.029100
## ENSG00000000005 98.513543
## ENSG00000000419  0.165160
## ENSG00000000457  0.174293
## ENSG00000000460  0.183204
## ENSG00000000938  0.153303

提取离散程度大于 1e-06 的基因的均值与离散程度，并保存供后续分析所用

# keep gene with dispersion greater than 1e-06
# extract the gene-wise mean and dispersion
meanDispPairs <- mcols(ddspickrell)[which(mcols(ddspickrell)$dispGeneEst > 1e-06), c("baseMean", "dispGeneEst")]
names(meanDispPairs) <- c("mean", "disp")
# save the mean-dispersion pairs for further usage
save(meanDispPairs, file="data/meanDispPairs.RData")

Simple Simulation

假设我们需要模拟一组数据，包含两组样本，每组4个重复，20000个基因，其中有20%的基因发生两倍差异表达。

# gene expression data simulation following negative binomial ditribution
load("data/meanDispPairs.RData")
# total number of genes
n <- 20000
# total sample size
m <- 8
# log fold changes used
es <- log2(2)
# formula
condition <- factor(rep(c("A","B"), each = m/2))
x <- model.matrix(~ condition)
# beta indicates the fold-change of genes
# 80% non-DE genes and 20% DE genes with up- or down-regulated
beta <- c(rep(0, n * 8/10), sample(c(-es,es), n * 2/10, TRUE))  

set.seed(123)
sf=rep(1,m)
# randomly sample mean-dispersion from real gene expression dataset
idx <- sample(nrow(meanDispPairs), n, replace=TRUE)
mu0 <- meanDispPairs[idx,1]
disp <- meanDispPairs[idx,2]
betafull <- cbind(log2(mu0), beta)
# betafull: [n,2]; t(x): [2,m]
# mu: [n,m]; mu matrix contains the mean of each genes for each sample
mu <- 2^(betafull %*% t(x))
head(mu)

##           1      2      3      4      5      6      7      8
## [1,]   2.55   2.55   2.55   2.55   2.55   2.55   2.55   2.55
## [2,]   0.90   0.90   0.90   0.90   0.90   0.90   0.90   0.90
## [3,]  79.34  79.34  79.34  79.34  79.34  79.34  79.34  79.34
## [4,]   4.81   4.81   4.81   4.81   4.81   4.81   4.81   4.81
## [5,] 368.22 368.22 368.22 368.22 368.22 368.22 368.22 368.22
## [6,] 719.99 719.99 719.99 719.99 719.99 719.99 719.99 719.99

# replicate mu matrix m times
muMat <- matrix(rep(mu, times=m) * rep(sf, each=n), ncol=m)
# n*m values to be simulated
mat <- matrix(rnbinom(n*m, mu=muMat, size=1/disp), ncol=m)
head(mat)

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,]    3    2    2    4    1    1    3    3
## [2,]    0    1    0    0    0    0    1    2
## [3,]   89   82  103  137   50   82   62   79
## [4,]    6    2    3    4    3    5    3    1
## [5,]  519  490  213  321  306  256  460  375
## [6,]  580 1459  829  855  508  340  647  410

简单可视化展示模拟数据。这里把不表达的基因去除后进行可视化，基因表达的直方图表明我们的数据整体还是符合负二项分布的

par(mfrow=c(2,4))
sapply(seq_len(ncol(mat)), function(i) hist(log(mat[rowSums(mat)>0,] + 1)[,i], main=paste('Col',i)))

par(mfrow=c(1,1))
boxplot(log(mat[rowSums(mat)>0,] + 1), main='Simple simulation')

与Pickrell数据进行比较。这里使用稍微严格一点的过滤条件

# pickrell dataset
counts.pick <- assay(pickrell)
counts.pick <- counts.pick[rowSums(counts.pick>0) > 10, ]
# demonstrate first 8 columns
par(mfrow=c(2,4))
sapply(1:8, function(i) hist(log(counts.pick[,i]+1), main=paste('Col',i)))

par(mfrow=c(1,1))
boxplot(log(counts.pick+1))

在Pickrell数据当中样本之间的整体表达情况会有一点波动，而我们模拟的数据没有。这可能是由于实验操作或者生物体之间的误差引起的，而模拟数据并不存在这种这种误差。我们也可以人为引入这种误差来更好地模拟真实数据。

Assessment

最后，我们使用 DESeq2 对模拟数据进行差异分析，测试模拟数据

# DE with DESeq2
e <- DESeqDataSetFromMatrix(mat,
                            colData = as.data.frame(condition),
                            design = ~condition)

## converting counts to integer mode

de <- DESeq(e)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

res <- results(de)
table(abs(res$log2FoldChange) >= 1)/nrow(res)

## 
## FALSE  TRUE 
## 0.710 0.232

table(abs(res$log2FoldChange) >= 1 & res$padj < 0.1)/nrow(res)

## 
## FALSE  TRUE 
## 0.753 0.048

summary(res)

## 
## out of 18828 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 663, 3.5%
## LFC < 0 (down)     : 640, 3.4%
## outliers [1]       : 62, 0.33%
## low counts [2]     : 6545, 35%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

使用模拟数据进行差异分析发现约23%的基因具有2倍以上的差异表达，但如果考虑 adjusted p-value < 0.1 的话，则只有5%左右的基因具有显著差异表达。同时，我们可以发现有37%的基因被认为是低表达（基因的平均counts < 5）的情况。

进一步，我们通过灵敏度（sensitivity）和特异性（specificity）评估 DESeq2 在模拟数据上的表现

Figure: Wiki

sensitivity, recall, hit rate, or true positive rate (TPR)

specificity, selectivity or true negative rate (TNR)

nonzero <- rowSums(counts(e)) > 0
# TPR = Sensitivity
sensidx <- abs(beta) > 0 & nonzero
sensitivity <- sum((res$padj<0.1)[sensidx], na.rm=TRUE)/nrow(res)
sensitivity

## [1] 0.0565

# FPR = 1 - Specificity
nonDEidx <- abs(beta) == 0 & nonzero
oneMinusSpecificity <- sum((res$padj<0.1)[nonDEidx], na.rm=TRUE)/nrow(res)
oneMinusSpecificity

## [1] 0.00865

以上只是单次模拟的结果，如果需要对软件的性能进行评估，我们应当进行多次模拟，然后再进行性能评估。

Ref:

DESeq2: https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8

DESeq2paper: https://www.huber.embl.de/DESeq2paper/

biostar handbook|How to simulate sequencing results: https://blog.actorsfit.com/a?ID=01000-a374125f-99d7-452a-99f9-61466d0a3ee6

完