6.4 真實數(shù)據(jù)集中的DE

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

6.4.1 簡介
為了測試不同的單細(xì)胞差異表達(dá)方法，我們將使用Blischak數(shù)據(jù)集。在該實驗中，除了單細(xì)胞數(shù)據(jù)之外，還生成了每個細(xì)胞系的Bulk RNA-seq數(shù)據(jù)。我們將使用標(biāo)準(zhǔn)方法在bulk數(shù)據(jù)上識別出的差異表達(dá)基因作為評估單細(xì)胞方法準(zhǔn)確性的標(biāo)準(zhǔn)。為了節(jié)省時間，我們已預(yù)先計算了這些數(shù)據(jù)?？梢赃\行以下命令來加載這些數(shù)據(jù)。

> 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個體的比較得出的標(biāo)準(zhǔn)?，F(xiàn)在加載相應(yīng)的單細(xì)胞數(shù)據(jù)：

> 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

現(xiàn)在我們將比較各種單細(xì)胞DE方法。這里我們將僅運行以R包形式提供且運行速度相對較快的方法。
6.4.2 Kolmogorov-Smirnov檢驗
非參數(shù)檢驗最容易實現(xiàn)。最常用的非參數(shù)檢驗是Kolmogorov-Smirnov檢驗（KS檢驗），我們可以用它來比較兩個個體中每個基因的分布。KS檢驗量化了兩個群體中每個基因表達(dá)的經(jīng)驗累積分布之間的距離。它對均值表達(dá)和變異性的變化很敏感。但它假設(shè)數(shù)據(jù)是連續(xù)的，當(dāng)數(shù)據(jù)包含大量相同值（例如零）時，其性能可能很差。KS檢驗的另一個問題是，它對大樣本量非常敏感，因此即使差異幅度非常小，也可能變得顯著。

雙樣本Kolmogorov–Smirnov統(tǒng)計量的圖示。紅線和藍(lán)線分別對應(yīng)一個經(jīng)驗分布函數(shù)，黑色箭頭為雙樣本KS統(tǒng)計量。

運行檢驗：

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

此代碼將函數(shù)“應(yīng)用”到表達(dá)矩陣、數(shù)據(jù)的每一行（由1指定）。在函數(shù)中我們只返回ks.test輸出的p.value。現(xiàn)在我們可以考慮KS檢驗檢測到了多少個真實陽性和陰性DE基因。
6.4.2.1 準(zhǔn)確性評估

> 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基因的數(shù)量較多，因此我們通常將這些計數(shù)標(biāo)準(zhǔn)化為真實陽性率（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遠(yuǎn)高于FPR，表明KS檢驗正在識別DE基因。
目前為止，我們只評估了單一顯著性閾值下的性能。改變閾值并在一系列值上進(jìn)行評估通常很有用。將其繪制為受試者工作特征曲線（ROC），并計算該曲線下的面積（AUC）作為準(zhǔn)確度統(tǒng)計值。使用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方法進(jìn)行比較，我們將此代碼放入函數(shù)中，這樣我們就不需要重復(fù)它了：

> 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秩和檢驗是另一種非參數(shù)檢驗，專門檢驗一組值是否大于/小于另一組值。因此，它通常被認(rèn)為是兩組間中位表達(dá)差異的檢驗；而KS檢驗對表達(dá)分布的任何變化都很敏感。

> 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
我們已經(jīng)在6.6章使用edgeR進(jìn)行差異表達(dá)分析。edgeR基于基因表達(dá)的負(fù)二項分布模型并使用廣義線性模型（GLM）框架，使我們能夠?qū)⑴蔚绕渌蛩丶{入模型中。

> 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基于零膨脹負(fù)二項分布模型。它使用結(jié)合了基因表達(dá)的離散和連續(xù)值的障礙度模型來檢驗差異表達(dá)。MAST也使用線性建?？蚣軄韺崿F(xiàn)復(fù)雜的模型。

> 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使用我們在上一章討論過的單細(xì)胞基因表達(dá)的Poisson-Beta模型，并將其與我們在使用edgeR時遇到的廣義線性模型相結(jié)合。BPSC將一個或多個組與對照組進(jìn)行比較，可包括模型中的批次等其他因素。

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是第一個單細(xì)胞DE方法，它使用貝葉斯統(tǒng)計將零膨脹負(fù)二項分布模型擬合到表達(dá)數(shù)據(jù)中。下面的用法檢驗不同群體中單個基因的平均表達(dá)的差異，但最近的版本包括檢驗一組基因——通常代表一條通路的平均表達(dá)或離散的差異的方法。

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)

往期內(nèi)容：
重生之我在劍橋大學(xué)學(xué)習(xí)單細(xì)胞RNA-seq分析——6. 生物學(xué)分析（1）
重生之我在劍橋大學(xué)學(xué)習(xí)單細(xì)胞RNA-seq分析——6. 生物學(xué)分析（2）
重生之我在劍橋大學(xué)學(xué)習(xí)單細(xì)胞RNA-seq分析——6. 生物學(xué)分析（3）