首先交代一下用来说明这几个统计量的例子。这里会使用基因表达作为一个例子。假设我们有两组细胞：对照组和处理组。我们正在研究基因 A 在处理的条件下是否受到表达或没有表达。每组我们有 12 个重复。我们通常做的是取每组 12 个重复的平均值，并进行 t 检验以比较差异是否显着（假设正态分布）。然后我们得到一个 p 值，比如 p = 0.035。如果它小于 0.05（所设置的阈值），我们得出结论，在处理后基因 A 的表达发生了显着变化。好，问题来了，p value为0.035告诉了我们怎样的信息？
首先我们需要，p值是从Null假设开始的：
H0 ：处理后基因 A 的基因表达没有差异。
和一个替代假设：
H1：处理后，基因 A 的表达发生变化。
每个 P 值的定义都是从假设Null假设为真开始的。p值为0.035，表示在Null情况下，我们看到处理后基因表达差异的概率为0.035，非常低。如果我们选择 0.05 的显着水平，则我们拒绝Null假设并接受备选假设。所以，如果你不能说明你Null假设是什么，你就无法理解 P 值。
对于典型的基因组研究，我们要比较数千个基因。我们如何报告包含差异表达基因的基因列表？我们可以对每个单个基因进行 a 检验，如果 p 值小于 0.05，我们将报告它。然而，它会给我们带来很多假阳性的结果，因为我们没有考虑多次测试。
让我们开始使用同时检测数千个基因的数据集。
###这里我用到哈佛大学统计的一个数据集
install_github('genomicsclass/GSE5859Subset')
##加载对应的包和数据
library(devtools)
library(qvalue)
library(GSE5859Subset)
data(GSE5859Subset)
### 查看数据
dim(geneExpression)
# [1] 8793 24
查看可用的数据和对象:
geneExpression[1:6, 1:6]
## GSM136508.CEL.gz GSM136530.CEL.gz GSM136517.CEL.gz
## 1007_s_at 6.543954 6.401470 6.298943
## 1053_at 7.546708 7.263547 7.201699
## 117_at 5.402622 5.050546 5.024917
## 121_at 7.892544 7.707754 7.461886
## 1255_g_at 3.242779 3.222804 3.185605
## 1294_at 7.531754 7.090270 7.466018
## GSM136576.CEL.gz GSM136566.CEL.gz GSM136574.CEL.gz
## 1007_s_at 6.837899 6.470689 6.450220
## 1053_at 7.052761 6.980207 7.096195
## 117_at 5.304313 5.214149 5.173731
## 121_at 7.558130 7.819013 7.641136
## 1255_g_at 3.195363 3.251915 3.324934
## 1294_at 7.122145 7.058973 6.992396
dim(sampleInfo)
## [1] 24 4
head(sampleInfo)
## ethnicity date filename group
## 107 ASN 2005-06-23 GSM136508.CEL.gz 1
## 122 ASN 2005-06-27 GSM136530.CEL.gz 1
## 113 ASN 2005-06-27 GSM136517.CEL.gz 1
## 163 ASN 2005-10-28 GSM136576.CEL.gz 1
## 153 ASN 2005-10-07 GSM136566.CEL.gz 1
## 161 ASN 2005-10-07 GSM136574.CEL.gz 1
sampleInfo$filename
## [1] 'GSM136508.CEL.gz' 'GSM136530.CEL.gz' 'GSM136517.CEL.gz'
## [4] 'GSM136576.CEL.gz' 'GSM136566.CEL.gz' 'GSM136574.CEL.gz'
## [7] 'GSM136575.CEL.gz' 'GSM136569.CEL.gz' 'GSM136568.CEL.gz'
## [10] 'GSM136559.CEL.gz' 'GSM136565.CEL.gz' 'GSM136573.CEL.gz'
## [13] 'GSM136523.CEL.gz' 'GSM136509.CEL.gz' 'GSM136727.CEL.gz'
## [16] 'GSM136510.CEL.gz' 'GSM136515.CEL.gz' 'GSM136522.CEL.gz'
## [19] 'GSM136507.CEL.gz' 'GSM136524.CEL.gz' 'GSM136514.CEL.gz'
## [22] 'GSM136563.CEL.gz' 'GSM136564.CEL.gz' 'GSM136572.CEL.gz'
head(geneAnnotation)
## PROBEID CHR CHRLOC SYMBOL
## 1 1007_s_at chr6 30852327 DDR1
## 30 1053_at chr7 -73645832 RFC2
## 31 117_at chr1 161494036 HSPA6
## 32 121_at chr2 -113973574 PAX8
## 33 1255_g_at chr6 42123144 GUCA1A
## 34 1294_at chr3 -49842638 UBA7
让我们看看一个单一的基因:
g<- sampleInfo$group
e<- geneExpression[25,]
# t-test, expression should be normal distribution
qqnorm(e[g==1])
qqline(e[g==1])
qqnorm(e[g==0])
qqline(e[g==1])
对其进行t-test：
t.test(e[g==1], e[g==0])
##
## Welch Two Sample t-test
##
## data: e[g == 1] and e[g == 0]
## t = 0.28382, df = 21.217, p-value = 0.7793
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.1431452 0.1884244
## sample estimates:
## mean of x mean of y
## 10.52505 10.50241
对所有基因进行 t 检验
mytest<- function(x) t.test(x[g==1], x[g==0], var.equal=T)$p.value
pvals<- apply(geneExpression, 1, mytest)
sum(pvals< 0.05) # 统计p值小于0.05的基因个数
## [1] 1383
看看 p 值分布
hist(pvals)
用随机数据模拟多重比较
m<- nrow(geneExpression)
n<- ncol(geneExpression)
# generate random numbers
randomData<- matrix(rnorm(n*m), m, n)
nullpvalues<- apply(randomData, 1, mytest)
hist(nullpvalues)
将此直方图与上面的直方图进行比较。看到了什么？即使我们随机生成数据，仍然会看到一些 p值 小于 0.05！！我们随机生成数据，应该没有表达的基因。然而，我们看到了几乎均匀分布的 p 值。
p 值是随机变量。在数学上，可以证明在原假设（并且满足一些假设，在这种情况下，检验统计量 T 遵循标准正态分布）下，p 值遵循均匀 (0,1) 分布，这意味着 P(p < p1) = p1。这意味着看到小于 p1 的 p 值的概率等于 p1。在 100 次 t 检验中，在 null（对照和处理之间没有差异）下，我们将看到 1 次检验的 p 值小于 0.01。我们将看到 2 个 p 值小于 0.02 的测试等等。这解释了为什么我们在随机生成的数字中看到一些 p 值小于 0.05。
一种方法是使用 Bonferroni 校正,来校正假阳性：仅当 P 值小于 alpha（通常为 0.05）除以比较次数（p < alpha/m）时，才将特定比较定义为具有统计显着性) 。例如我们计算了 100 次 t 检验，得到了 100 个 p 值，我们只认为 p 值小于 0.05/100 的基因是显着的。这种方法非常保守，一般用于全基因组关联研究 (GWAS)。
或者，我们可以使用错误发现率 (FDR) 来报告显着表达基因。在我们这个例子中， FDR= 0.05。这意味着在所有发现中，预计有 5% 是假阳性。Benjamini & Hochberg（BH 方法）在 1995 年提出了一种控制 FDR 的方法：让 k 成为最大的 i，使得p(i) <= (i/m) * alpha，（m 是比较次数）然后拒绝 H(i) for i =1, 2, …k
此过程控制FDR的alpha 级别。该方法为每次比较设置不同的阈值 p 值。假设我们计算了 100 个 t 检验，得到了 100 个 p 值，我们希望控制 FDR = 0.05。然后我们将 p 值从小到大排序。如果 p(1) <= 1/100 * 0.05，则我们拒绝null假设并接受替代的假设。如果 p(2) < = 2/100 * 0.05，则我们拒绝 null 并接受替代的假设。
## order the pvals computed above and plot it.
alpha = 0.05
m = length(pvals)
#m is the number of 8793 comparisons
plot(x=seq(1,100), y=pvals[order(pvals)][1:100])
abline(a=0, b=alpha/m)
title('slop is alpha/m')
# let's zoom in to look at the first 15 p values from small to big
plot(x=seq(1,100), y=pvals[order(pvals)][1:100], xlim=c(1,15))
abline(a=0, b=alpha/m)
title('slop is alpha/m')
# we can see that the 14th p value is bigger than its own threshold
# which is computed by (0.05/m) * 14 = 7.960878e-05
# we will use p.adjust function and the method 'fdr' or 'BH' to
# correct the p value, what the p.adjust function does to to
# recalculate the p-value. ?p.adjust to see more
# p(i)<= (i/m) * alpha
# p(i) * m/i <= alpha
# we can then only accept the returned if p.adjust(pvals) <= alpha
# number of p values smaller than their own thresholds after controlling FDR=0.05
我们可以看到第 14 个 p 值大于它自己的阈值，由 (0.05/m) * 14 = 7.960878e-05 计算我们将使用 p.adjust 函数和方法“fdr”或“BH”来纠正p 值，p.adjust 函数的作用是重新计算 p 值。p(i)<= (i/m) * alpha p(i) * m/i <= alpha 如果 p.adjust(pvals) <= alpha 我们只能接受返回的 p 值.
sum( p.adjust(pvals, method='fdr') < 0.05 )
## [1] 13
13，和我们从图中看到的一样。
Storey 在 2002 年提出的另一种方法是 FDR 的直接方法：让 K 是最大的 i，使得 pi_0 * p(i) < (i/m) * alpha 然后拒绝 H(i) for i =1,2,…k pi_0 是对基因列表中原假设为真的比例的估计，范围从 0 到 1。所以当 pi_0 为 1 时，我们有 Benjamini & Hochberg 校正。这种方法不如 BH 方法保守。使用 bioconductor 包“qvalue”中的 qvalue 函数
sum( qvalue(pvals)$qvalues < 0.05)
## [1] 22
它是 22，不如 BH 方法保守。请注意，FDR 是基因列表的一个属性。q 值是为特定基因定义的：
“但如果你确实想为每个基因分配一个数字，你可以做的一件简单的事情就是,可以一个基因一个基因，然后决定最小FDR4，这将包括这个基因在列表中。一旦你这样做了，那么你就定义了一个 q 值。这是基因列表中经常报告的内容”
“为了定义 q 值，我们对通过 p 值测试的特征进行排序，然后计算具有最重要、两个最重要、三个最重要等的列表的 FDR……列表的
FDR，例如， m 个最重要的测试被定义为第 m 个最重要特征的 q 值。换句话说，特征的 q 值是包含该基因的最大列表的 FDR”
解读到这里，希望这篇文章能帮助你更好地理解 p 值、FDR 和 q 值。
文章整理翻译于该博客：
https://divingintogeneticsandgenomics.rbind.io/post/understanding-p-value-multiple-comparisons-fdr-and-q-value/
联系客服
