额其实挺难的。

什么是QQ-plot

qqplot是通过比较两个概率分布的分位数来评价这俩分布的相似性的

然后统计学上常举的例子是标准正态分布
就假说有个标准正态分布
有个标准化了的待测的分布

分位数就是落在把概率分布分成几个等概率的部分的位置的数，比如二分位数就是中位数，第二四分位数也是中位数。
QQ-plot的x轴是这么来的：你从标准正态分布中拿出N个分位数，这N个分位数对应的值可以查表算到。
QQ-plot的y轴是这么来的：然后抱着这个N分位数算出来的标准值，去待测分布上查横轴。

N越大，那么结果越精确，假如是同分布，那么Y=X；
在标准正态分布的QQ-plot情况下：
假如是左偏（左偏是峰靠右，偏度是偏离对称性程度的意思）分布，凹曲线；
假如是右偏分布，凸曲线；

GWAS的QQ-plot

GWAS里面的QQ-plot的横坐标是log(p-value)，p-value这东西在零假设下是个均匀分布

p-value 在零假设下是(0, 1)均匀分布

p-value是谁的p-value？

p-value肯定是构造统计量的p-value，而且是用β构建的统计量，然后这个统计量还服从一个标准的分布，卡方或者t或者F或者别的什么，这里不同的模型导出p-value的方法不一样
Wald’s test，p-value算的是β^/SE\hat{\beta}/SE的p-value，以它为例，其它也一样，毕竟都是f(β^)f(\hat{\beta})

p-value的意义是？

p-value就是指，假如β是0，你抽出的数据集算出不小于现在你已经算出来的β^\hat{\beta}（的绝对值）的β^\hat{β}的概率有多大

这里p-value的概率分布就是：
Pr(P≤x)=Pr(test statistic falls within the most extreme region of prob. mass x∣NULL)=xPr(P≤x)=Pr(\text{test statistic falls within the most extreme region of prob. mass x}|NULL)=x
这就很显然是个(0, 1)上的均匀分布

零假设下的−log10(p-value)-log10(\text{p-value})是指数分布

这就直接R证明：

1
2
3
4
5
6
7
8
9
10
11
# 生成遵循均匀分布的随机数据
min_val <- 0  # 均匀分布的最小值
max_val <- 1  # 均匀分布的最大值
n <- 1000     # 生成的数据点数量
data <- runif(n, min = min_val, max = max_val) # 均匀分布
mydata <- -log(data)
theory_quantiles <- qexp(seq(0.01, 0.99, by = 0.01), rate = 1)
# 创建QQ图比较正偏态和标准指数分布
par(mfrow = c(1, 2))  # 设置图形布局为一行两列
qqplot(theory_quantiles, quantiles_standard, main="QQ图", ylab="数据的分位数", xlab="理论分位数", pch=19)
abline(0, 1, col = "red", lty = 2)

GWASQQ-plot_image_1.png
无证明，rate大概是2.3上下的样子。

几种情况下的QQ-plot曲线

为方便写，取对数底数为2，不证明底数不影响结果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
set.seed(123)

mydata <- -log(data)
data_standard <- rexp(1000)  # 标准指数分布数据
# gamma分布
shape_param_positive <- 2  # 形状参数，大于1产生正偏态
data_positive <- rgamma(1000, shape = shape_param_positive)  # 正偏态数据

# 通过使用较小的rate值的指数分布
rate_negative <- 0.5 
data_negative <- rexp(1000, rate = rate_negative)

quantiles_standard <- quantile(data_standard, probs = seq(0.01, 0.99, by = 0.01))
quantiles_negative <- quantile(data_negative, probs = seq(0.01, 0.99, by = 0.01))
quantiles_positive <- quantile(data_positive, probs = seq(0.01, 0.99, by = 0.01))
quantiles_mydata <- quantile(mydata, probs = seq(0.01, 0.99, by = 0.01)) 
# 定义理论指数分布的分位数
theory_quantiles <- qexp(seq(0.01, 0.99, by = 0.01), rate = 1)

# 创建QQ图比较正偏态和标准指数分布
par(mfrow = c(2, 2))  # 设置图形布局为一行两列
qqplot(theory_quantiles, quantiles_standard, main="标准指数分布的QQ图", ylab="数据的分位数", xlab="理论分位数", pch=19)
abline(0, 1, col = "red", lty = 2)
qqplot(theory_quantiles, quantiles_positive, main="Gamma分布的QQ图", ylab="数据的分位数", xlab="理论分位数", pch=19)
abline(0, 1, col = "red", lty = 2)
qqplot(theory_quantiles, quantiles_negative, main="小rate指数分布的QQ图", ylab="数据的分位数", xlab="理论分位数", pch=19)
abline(0, 1, col = "red", lty = 2)
qqplot(theory_quantiles, quantiles_mydata, main="我的数据和指数分布的QQ图", ylab="数据的分位数", xlab="理论分位数", pch=19)
abline(0, 1, col = "red", lty = 2)

hist(quantiles_standard, breaks="Sturges", col="lightblue", border="black", main="标准指数分布", xlab="Value", ylab="Frequency",xlim=c(0,4),ylim=c(0,40))
hist(quantiles_positive, breaks=10, col="lightblue", border="black", main="Gamma分布", xlab="Value", ylab="Frequency", xlim=c(0,4),ylim=c(0,40))
hist(quantiles_negative, breaks=20, col="lightblue", border="black", main="小rate指数分布", xlab="Value", ylab="Frequency", xlim=c(0,4),ylim=c(0,40))
hist(quantiles_mydata, breaks="Sturges", col="lightblue", border="black", main="我的数据", xlab="Value", ylab="Frequency", xlim=c(0,4),ylim=c(0,40))

GWASQQ-plot_image_2.pngGWASQQ-plot_image_3.png

什么情况下QQ-plot曲线会掉下去？

GWASQQ-plot_image_4.png
为什么这种情况不会发生？
因为外面套的那个指数分布的分布函数的λ是来自于-log10的，那玩意就是0.1，不会变动。

前面一段符合理论值，后一段翘起是怎么回事？

大部分的SNP，pvalue都和理论的差不多，有几个pvalue实际值非常小，也就是它的表型和基因型存在较大相关性，这种情况在零假设下，发生的可能性非常小。在图上就体现为尾部高高翘起。

前面一段符合理论值，后一段比理论值还低怎么回事？

说明算出的pvalue比零假设还高，这有可能是模型的问题。

模型认为存在很多本身和品种关联的位点（通过亲缘关系和群体结构判断的），倘若在亚群中表型和基因型关联，也没什么大不了的，pvalue不该这么低，因此对它们作了校正，人工增大pvalue，从而控制假阳性。掉下去了可能是把真阳性都控制掉了。