6.4 真实数据集中的DE

> library(scRNA.seq.funcs)
> library(edgeR)
> library(MAST)
> library(ROCR)
> set.seed(1)

6.4.1 简介
为了测试不同的单细胞差异表达方法，我们将使用Blischak数据集。在该实验中，除了单细胞数据之外，还生成了每个细胞系的Bulk RNA-seq数据。我们将使用标准方法在bulk数据上识别出的差异表达基因作为评估单细胞方法准确性的标准。为了节省时间，我们已预先计算了这些数据。可以运行以下命令来加载这些数据。

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

以上是对NA19101和NA19239个体的比较得出的标准。现在加载相应的单细胞数据：

> molecules <- read.table("data/tung/molecules.txt", sep = "\t")
> anno <- read.table("data/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

现在我们将比较各种单细胞DE方法。这里我们将仅运行以R包形式提供且运行速度相对较快的方法。
6.4.2 Kolmogorov-Smirnov检验
非参数检验最容易实现。最常用的非参数检验是Kolmogorov-Smirnov检验（KS检验），我们可以用它来比较两个个体中每个基因的分布。KS检验量化了两个群体中每个基因表达的经验累积分布之间的距离。它对均值表达和变异性的变化很敏感。但它假设数据是连续的，当数据包含大量相同值（例如零）时，其性能可能很差。KS检验的另一个问题是，它对大样本量非常敏感，因此即使差异幅度非常小，也可能变得显著。

双样本Kolmogorov–Smirnov统计量的图示。红线和蓝线分别对应一个经验分布函数，黑色箭头为双样本KS统计量。

运行检验：

> 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")

此代码将函数“应用”到表达矩阵、数据的每一行（由1指定）。在函数中我们只返回ks.test输出的p.value。现在我们可以考虑KS检验检测到了多少个真实阳性和阴性DE基因。
6.4.2.1 准确性评估

> 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
# Number of KS-DE genes that are truly not-DE

KS检验将更多的真实阴性基因（假阳性）鉴定为DE，而不是真实阳性基因（真阳性）。这可能是由于notDE基因的数量较多，因此我们通常将这些计数标准化为真实阳性率（TPR）、TP/(TP+FN)和假阳性率（FPR）、FP/(FP+TP)。

> 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

可以看到TPR远高于FPR，表明KS检验正在识别DE基因。
目前为止，我们只评估了单一显著性阈值下的性能。改变阈值并在一系列值上进行评估通常很有用。将其绘制为受试者工作特征曲线（ROC），并计算该曲线下的面积（AUC）作为准确度统计值。使用ROCR包来简化此绘制。

# 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)

KS检验的ROC曲线。

> aucObj <- ROCR::performance(pred, "auc")
> aucObj@y.values[[1]]  #AUC
[1] 0.7954796

最后，为了方便与其他DE方法进行比较，我们将此代码放入函数中，这样我们就不需要重复它了：

> 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]])
  }

6.4.3 Wilcox/Mann-Whitney-U检验
Wilcox秩和检验是另一种非参数检验，专门检验一组值是否大于/小于另一组值。因此，它通常被认为是两组间中位表达差异的检验；而KS检验对表达分布的任何变化都很敏感。

> 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)
[1] 0.8320326

Wilcox检验的ROC曲线。

6.4.4 edgeR
我们已经在6.6章使用edgeR进行差异表达分析。edgeR基于基因表达的负二项分布模型并使用广义线性模型（GLM）框架，使我们能够将批次等其他因素纳入模型中。

> 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)
[1] 0.8466764

edgeR的ROC曲线。

6.4.5 MAST
MAST基于零膨胀负二项分布模型。它使用结合了基因表达的离散和连续值的障碍度模型来检验差异表达。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(~ 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)
[1] 0.8284046

MAST的ROC曲线。

6.4.6 运行速度慢的方法（运行时间超过1小时）
以下方法目前运行速度太慢，鼓励您自行尝试。
6.4.7 BPSC
BPSC使用我们在上一章讨论过的单细胞基因表达的Poisson-Beta模型，并将其与我们在使用edgeR时遇到的广义线性模型相结合。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)

6.4.8 SCDE
SCDE是第一个单细胞DE方法，它使用贝叶斯统计将零膨胀负二项分布模型拟合到表达数据中。下面的用法检验不同群体中单个基因的平均表达的差异，但最近的版本包括检验一组基因——通常代表一条通路的平均表达或离散的差异的方法。

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)

往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（1）
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（2）
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（3）