第七章 差异表达

7.1 动机

有许多不同的方法可以探测差异表达基因。本章将探索这些不同的方法，而不是描述一个标准的方法。目的是给你一个这些方法现存的思想，这样你在自己的分析中可以自己选择合适的方法。

7.1.1 逐个基因分析的方法

截至目前的实践都是逐个基因的进行差异分析，这忽视了基因之间的依赖关系，比如它们之间在调节模型上的相互作用。很明显这并不满足我们的需要，而且随着学习的更多，这一点也会改变。本章中我们将会覆盖一个基因一个基因分析的方法。

7.1.2 非特异性过滤

大多数的微阵列包含的探针对应着比差异表达更多的基因。实际上，均一化的一个基本假设是大部分基因都不是差异表达基因。为了减少单个基因假设检验的强大的多样性带来的损失，我们建议进行一些形式的非特性过滤。

非特异性意思是说这一步过滤没有使用检测的RNA样本的参数或者控制条件。它的目的是移除在任何对比下都没有出现差异表达的探针集。我们发现这在变异基础上选择基因时非常有用。只有在样本间表现出任何差异的基因才有可能在感兴趣的组间的差异表达。

7.1.3 倍数变化和t检验

选择基因最简单的方法是倍数变化准则。这可能在只能获得少数几个复制本时可行。一个只是基于倍数变化的分析，在存在生物或者实验差异的情况下，排除了存在观察值变化的可能，而这种可能型在不同的基因间是不一样的。这是用统计检验来进行差异表达分析的主要原因。

总的来说，一个人会看到一个基因在不同条件下表达水平分布的的各种差异。通常会考虑位置参数（比如：平均值或者中位数）。这使得我们使用t检验和它的各种变异。通常我们也有很好的理由去去考虑关于分布的其它性质，比如受试者工作特征曲线
下的区域面积作为分类阈值。

或许你可以区分参数检验和非参数检验，比如t检验和Mann-Whitney检验或者排列检验。如果近似满足基本的模型假设，比如t检验的正态分布，参数检验将会产生更有效的效果。非参数检验不需要满足数据分布的严格假设。在很多微阵列研究中，一个小的样本量常常导致非参数检验的效度不足。对这种情况的几个程序性的解决方法是进行参数检验，但是谨慎的使用p值结果来对基因差异表达的证据进行排序。

7.2 非特异性过滤

我们现在加载第一章中用过的数据集，你可以发现此处我们用到的对急性淋巴白血病数据的一个综合描述。

library("ALL")
data("ALL")

首先我们创建一个B系肿瘤的样本列表。

bcell = grep("^B", as.character(ALL$BT))

BCR/ABL易位t(9;22)(q34;q11)，常被叫做费城染色体，产生一个由BCR和ABL1组成的融合基因，在急淋白血病中相对比较常见且和治疗相关。这里我们关注携带此易位突变的样本且和非发现常见细胞遗传突变的样本（NEG）进行对比。

moltyp = which(as.character(ALL$mol.biol)%in% c("NEG", "BCR/ABL"))

现在我们创建一个新的对象ALL_bcrneg，它只包含满足这两个条件的样本。

ALL_bcrneg = ALL[, intersect(bcell, moltyp)]
ALL_bcrneg$mol.biol = factor(ALL_bcrneg$mol.biol)

第二行语句会清除ALL_bcrneg$mol.biol这一因子变量中的空水平变量。

现在如果我们基于差异进行过滤，我们首先需要确保这些变异不是由它们对基因的平均表达水平所主导的。如果是，那么基于变异进行选择会与基于表达绝对i水平进行选择相混杂。因为探针序列特异性的背景和增益因子影响，我们有充分的理由不去使用绝对水平进行基因筛选。为了检验相关性，我们对行均值和行标准差进行画图，并对他们的回归进行平滑估计。

library("vsn")
meanSdPlot(ALL_bcrneg)

练习 7.1

• a. 评价该散点图，你认为平均值和标准差之间是弱相关吗？
• b. 查看meanSdPlot的帮助文档，rank参数的作用是什么？
  我们假设关联不强，下一步是去除低差异的探针集。在下面的代码中，我们将会去除最低的80%的探针集。我们选择这么大一部分是想要减小接下来的计算长度。这一部分最好的选择取决于探针设计和生物样本，但在实际中这会更低。

sds = esApply(ALL, 1, sd)
sel = (sds > quantile(sds, 0.8))
ALLset1 = ALL_bcrneg[sel, ]

这种方法一个潜在的不利之处会发生在我们感兴趣实验因素下的一个分组的样本量较少时。在这种情况下，一个在该组和其它组间差异表达的基因也许没有较大的标准差。那你应该怎么处理这种情况？

在这里你也许想看一下数据热图看它们是否存在明显的模式。?heatmap可以查询该函数的主页。

差异基因##

在Bioconductor中，genefilter包可以通过不同的过滤方法简单选择基因。另外，在一些其它的检验和对比中，我们已经发展出了快速的版本。这些包括rowttests，它对基因表达矩阵的每一行进行t检验；rowFtests使用F-tests；rowQ对每一行计算分位数。
首先也许最简单的是使用t-test。

library("genefilter")
tt = rowttests(ALLset1, "mol.biol")

查询rowttests主页可以看到返回值tt的四个元素的意义。许多实践者已经了解到小的p值并不总是反应出有较大变化的基因。让我们看看所谓的独木舟点图。

plot(tt$dm, -log10(tt$p.value), pch=".", xlab =expression(mean~log[2]~fold~change), ylab = expression(-log[10](p)))

练习 7.2
用t检验结果查看有多少探针对应着差异表达基因。

7.2 多检验

多检验已经受到大量关注。我们在multtest包中提高了一个简短的函数功能介绍(Pollard et al., 2005)。

multtest包中使用的算法依赖于样本的随机排列。排列序数由参数B提供。接下来我们调用mt.maxT进行排列检验，用Welch方法。

library("multtest")
cl = as.numeric(ALLset1$mol.biol=="BCR/ABL")
resT = mt.maxT(exprs(ALLset1), classlabel=cl, B=1000)
ord = order(resT$index) ## the original gene order
rawp = resT$rawp[ord] ## permutation p-values

图示显示了向量rawp的未调整的序列p值的直方图。小p值的高占比提示大部分的基因的确是组间差异表达的。

hist(rawp, breaks = 50, co l= "#B2DF8A")

为了控制整体误差率（FWER）—显著基因集中的至少一个假阳性发生的概率—mt.maxT使用基于排序的maxT程序（Westfall and Young (1993)）。我们获得调整后p值小于0.05的34个基因。

sum(resT$adjp<0.05)
[1] 34

这个数字和直方图中的最左侧的高度相比显示我们丢失了大量的差异表达基因。FWER是一个严格的标准，在一些微阵列实验中，只有少量基因在此意义上是显著的，即使有很多是真正差异表达的。一个更敏感的标准有假性发现率提供的，就是所谓的显著基因中的假阳性的期望占比。我们可以使用 Benjamini和Hochberg
(1995) 提供的程序multtest来控制FDR：

res = mt.rawp2adjp(rawp, proc = "BH")
sum(res$adjp[,"BH"]<0.05)
[1] 209

7.5 温和的统计检验和limma包

t检验可以由limma包中的函数完成。首先，我们需要定义设计矩阵。一个可能是代表所有样本的一个基因的log2值的均值的截距项（下面矩阵的第一列，由ls组成），在第二列中编码组间差异。

library("limma")
design = cbind(mean = 1, diff = cl)

接下来，用线性模型用lmFit函数对每一个基因进行适配，可以用eBayes函数进行标准错误的经验贝尔斯递归（ Smyth (2004)）。这需要所有基因的信息来对达到每一个基因差异的更稳定的评估。

fit = lmFit(exprs(ALLset1), design)
fit = eBayes(fit)

我们可以用函数topTable列出10个最显著的差异基因。p值最小的三个探针集都映射到ABL1基因，这是t(9;22)(q34;q11)易位引起的融合基因的一部分，这些基因都过表达并在急性淋巴瘤中作为一个强的原癌基因。

library("hgu95av2")
ALLset1Syms = unlist(mget(featureNames(ALLset1),env = hgu95av2SYMBOL))
topTable(fit, coef = "diff", adjust.method = "fdr",sort.by = "p", genelist = ALLset1Syms)

当你比较这些结果的p值和那些参数检验得到的p值，你会看到他们几乎是一样的：

plot(-log10(tt$p.value), -log10(fit$p.value[, "diff"]), xlab = "-log10(p) from two-sample t-test", ylab = "-log10(p) from moderated t-test (limma)", pch=".")
abline(c(0, 1), col = "red")

结果见图，因为样本量比较大，经验贝尔斯检验在样本量较小时非常有用。让我们为我们数据的每一个组的三个阵列画出亚组：

subs = c(35, 65, 75, 1, 69, 71)
ALLset2 = ALL_bcrneg[, subs]
table(ALLset2$mol.biol)
BCR/ABL NEG
3 3

我们重复之前用过的检验步骤。

tt2 = rowttests(ALLset2, "mol.biol")
fit2 = eBayes(lmFit(exprs(ALLset2), design=design[subs, ]))

plot(-log10(tt2$p.value), -log10(fit2$p.value[, "diff"]), xlab = "-log10(p) from two-sample t-test", ylab = "-log10(p) from moderated t-test (limma)", pch=".")
abline(c(0, 1), col = "red")

图见如下：

让我们看看在常规t检验中小p值但在温和检验中较大的基因。

g = which(tt2$p.value < 1e-4 & fit2$p.value[, "diff"] > 0.02)

我们用不同的符号和颜色表示样本类别的表达值。

sel = (ALLset2$mol.bio == "BCR/ABL")+1
col = c("black", "red")[sel]
pch = c(1,16)[sel]
plot(exprs(ALLset2)[g,], pch=pch, col=col, ylab="expression")

结果如图：

7.6 ROC进行基因筛选

在这一部分我们探查一个发现差异表达基因的方法，不同于之前学习过的假设检验，这里使用基于分类的方法。这个方法由Pepe提出（2003）。

其目的是发现可能作为潜在标记的基因，也就是单个的基因的表达水平能够在组间进行消除。我们可以把其中一个称作对照组另一个是疾病组。如果我们用x轴表示给定基因的表达水平，我们采用简单的分类法x >= θ的样本，预测在疾病组，而x<θ的样本在对照组。对于不同的θ的选择，分层法可由它的特异性测量，敏感度由疾病组层中的真实疾病样本的占比测量。敏感度点图p = 1 − 特异度被叫做ROC曲线，例子见图：

我们要去鉴定最有可能探测样本携带BCR/ABL易位样本的基因。这个可由ROC曲线下面积（AUC）来表达，或者更一般的，pAUC。pAUC标准，对于一个小的p值，如p = 0.2，比AUC相关性更好，因为对于实际的诊断性标记，在考虑敏感度之前，我们需要高的特异度，比如高于80%。

类似于rowttests，genefilter包包含一个进行行的pAUC统计的函数。这个函数用ExpressionSet对象、因子变量名字符串、数值标量p作为参数。另外，这个函数还有参数flip。如果是TRUE（默认值），那么每一个基因的分层规则x<θ和x>θ将会被检测，这两者的曲线下的区域面积也会返回。在我们的数据中，疾病层是随意的且不一定总是关联着高表达值。实际上我们希望发现每组中的高表达和低表达基因。你可以设定flip为FALSE，如果你希望用x>θ规则来区别疾病。对于更进一步的分层方法，你也许可以考虑ROC包提供的功能。

rocs = rowpAUCs(ALLset1, "mol.biol", p=0.2)

下一步我们选择pAUC最大值的探针集，然后画出相关的ROC曲线。

j = which.max(area(rocs))
plot(rocs[j], main = featureNames(ALLset1)[j])

结果如上图左部。
练习 7.3
画出探针1636_g_at的表达值，如上图右部。
练习 7.4
如果一个基因在组间不显示差异表达，那么它的ROC曲线应该是怎样的？“理想”的ROC曲线应该怎样？为什么？

7.7 当增加效益

在7.3和7.6部分，我们展示两个发现差异表达基因的两条不同标准：基因表达分布的区位变化，简单阈值分层的ROC曲线。我们现在考虑这些方法并着重样本大小的影响。

我们想看样本量怎么影响两种方法发现的差异表达基因的数量。在这里，我们在我们的数据中重复取样本子集，大小不等，对每一个子集进行差异表达计算，每一次计算差异表达基因的确定数目。
让我们先围绕rowttests建立一个封装函数。

nrsel.ttest = function(x, pthresh=0.05) {pval = rowttests(x, "mol.biol")$p.value
return(sum(pval < pthresh))
}

还有rowpAUCs。

nrsel.pAUC = function(x, pAUCthresh=2.5e-2) {pAUC = area(rowpAUCs(x,fac="mol.biol", p=0.1))
return(sum(pAUC > pAUCthresh))
}

注意阈值的选择总是随机的，而t检验的pthresh不能和rowpAUCs、pAUCthresh直接对比。

我们需要一个对数据集进行大小不定的重复取样并画出结果。代码是直接的，但是有些无聊，为了简化我们在包BiocCaseStudies中提供了这个函数，它是resample。你可以输入resample来产看它的代码。

因为接下来的计算需要一点时间，我们构建了一个只包含1000个探针的ALLset1的子集来进行证明。（如果你计算机性能强大，你可以在下面的代码中用ALLset1代替x来运行。）

library(BiocCaseStudies)
x = ALLset1[sample(nrow(ALLset1), 1000), ]
}

resample(x, "nrsel.ttest")
resample(x, "nrsel.pAUC")
}

一次用nrsel.ttest运行重复取样函数，一次用nrsel.PAUC运行。你也需要传递**xII。图示如下：

练习 7.5
两条标准哪一条看起来更合理？

进一步学习的主题：

• 在本章中你学习了探测差异表达基因的几种方法，其它选择感兴趣基因的方法有什么？
• 你在这里看到的不同方法得出的结果重叠结果如何？在那些例子中他们产生不同- 结果的原因是什么？
• 你怎么制定标准来决定选择哪个方法比其它方法更好？
• 不同的方法是怎么处理缺失值和离群值？