QTL-seq系列 | 如何计算置信区间

书接上文,在计算出ΔSNP-index之后如何判断该位点是否存在QTL?这里通过随机模拟实验来确定ΔSNP-index的置信区间(H0: 没有QTL)。

计算原理

计算原理

模拟取样过程

双亲杂交(AA × BB)、自交后产生包含大量子代个体的分离群体(RIL或F2),随机抽取两组若干个体构成两个模拟bulk pool。对于构成bulk pool的每一个个体可能的基因型以及概率是固定的,RIL群体为AA:BB=0.5:0.5,而F2群体为AA:AB:BB=0.25:0.5:0.25,由此可以计算出给定子代数目的bulk pool中的等位基因频率P(A)和P(B)。

模拟计算ΔSNP-index

前面计算出了bulk pool中的等位基因频率,那么在特定测序深度(N)下等位基因A出现的次数就相当于进行N次二项分布取样 n(A)=XB(N,P(A))n(A) = X \sim B(N,P(A)),进而可以得到 SNP index=n(A)/NSNP\ index = n(A)/{N},两个bulk pool的SNP index相减得到ΔSNP-index。

计算零假设下的ΔSNP-index分布

将以上过程重复若干次(10000次)可以得到在零假设下的ΔSNP-index的分布,取其上分位数和下分位数可以得到对应置信区间(如95%置信区间取上0.025分位数和下0.025分位数)。

代码实现

模拟计算ΔSNP-index置信区间

这里使用R语言实现上述过程。

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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
# 加载R包
library(tidyverse)
library(cowplot)
# 定义基本信息
popType <- "RIL" # 群体类型,RIL或F2
bulkSizeH <- 30 # 高值表型bulk pool个体数
bulkSizeL <- 30 # 低值表型bulk pool个体数
minDepth <- 5 # 最低覆盖深度
maxDepth <- 150 # 最高覆盖深度
repN <- 10000 # 模拟实验重复次数
# 定义数据框,每一行代表两个bulk pool的不同覆盖深度
dltIndex_CI <- tibble(HB.DP = rep(minDepth:maxDepth, times = maxDepth-minDepth+1),
LB.DP = rep(minDepth:maxDepth, each = maxDepth-minDepth+1),
CI95upper = NA, CI95lower = NA, CI99upper = NA, CI99lower = NA)
# 定义一个进度条,用于显示计算进度
width <- options()$width
pb <- progress::progress_bar$new(
format = 'Progress [:bar] :percent eta: :eta',
total = nrow(dltIndex_CI), clear = FALSE, width = width
)
# 对不同覆盖深度循环计算置信区间,即循环dltIndex_CI的每一行
for (i in 1:nrow(dltIndex_CI)) {
depthH <- dltIndex_CI[i, 1][[1]]
depthL <- dltIndex_CI[i, 2][[1]]
# 模拟计算high bulk pool和low bulk pool的等位基因频率,重复repN次
if (popType == "RIL") {
PH <- rbinom(repN, bulkSizeH, 0.5) / bulkSizeH
PL <- rbinom(repN, bulkSizeL, 0.5) / bulkSizeL
} else if (popType == "F2") {
PH <- apply(rmultinom(repN, bulkSizeH, c(1, 2, 1)) * c(1, 0.5, 0) / bulkSizeH, 2, sum)
PL <- apply(rmultinom(repN, bulkSizeL, c(1, 2, 1)) * c(1, 0.5, 0) / bulkSizeL, 2, sum)
}
# 模拟计算两个bulk pool的ΔSNP-index,重复repN次
indexH <- rbinom(repN, depthH, PH) / depthH
indexL <- rbinom(repN, depthL, PL) / depthL
dltIndex <- indexH - indexL
# 取上分位数和下分位数作为置信区间
dltIndex_CI[i, c("CI95upper", "CI95lower", "CI99upper", "CI99lower")] <- t(
c(quantile(dltIndex, 0.975), quantile(dltIndex, 0.025),
quantile(dltIndex, 0.995), quantile(dltIndex, 0.005))
)
# 打印进度条
pb$tick()
}
# 保存结果,导出指定群体类型、bulk pool大小下不同覆盖深度的置信区间
save(dltIndex_CI,
file = paste("QTLseqCI",
paste(paste("indH", bulkSizeH, sep = ""),
paste("indL", bulkSizeL, sep = ""),
popType, sep = "_"),
paste("Depth", minDepth, maxDepth, sep = "_"),
paste("Rep", repN, sep = "_"),
"RData", sep = "."))

ΔSNP-index置信区间可视化

由此我们得到了特定群体类型、特定样本量大小下,不同测序深度位点的ΔSNP-index的95%和99%的置信区间,把结果导出保存为RData文件,方便下次使用而不用重复计算。ΔSNP-index的置信区间和覆盖深度的关系是怎样的呢?我们通过将上面结果可视化,将横坐标、纵坐标分别表示为high bulk pool、low bulk pool的位点覆盖深度,将ΔSNP-index的置信区间上限和下限值表示为不同颜色。由下面绘图结果可以看出,当提高位点覆盖深度时,ΔSNP-index的置信区间范围将会缩小,相反置信区间范围会增大。按照同样的思路我们也可以探究ΔSNP-index置信区间和群体类型、每个pool取样多少的关系。

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
35
36
37
# 将置信区间和覆盖深度的关系可视化
df <- gather(dltIndex_CI, type, CI, -HB.DP, -LB.DP) %>%
mutate(level = if_else(str_detect(type, "95"), "95 CI", "99 CI"),
direction = if_else(str_detect(type, "upper"), "Upper CI", "Lower CI"))
df$direction <- factor(df$direction, levels = c("Upper CI", "Lower CI"))
P_CI <- ggplot(df, aes(x = HB.DP, y = LB.DP, fill = CI)) +
geom_tile() +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
scale_fill_distiller(palette = "RdBu") +
labs(title = paste("CI",
paste(paste("indH", bulkSizeH, sep = ""),
paste("indL", bulkSizeL, sep = ""),
popType, sep = "_"),
paste("Depth", minDepth, maxDepth, sep = "_"),
paste("Rep", repN, sep = "_"), sep = "."),
fill = NULL) +
facet_grid(direction ~ level) +
theme_half_open() +
theme(strip.background = element_rect(fill = "#90EE90"),
plot.title = element_text(hjust = 0.5))
ggsave(P_CI, filename = paste("QTLseqCI",
paste(paste("indH", bulkSizeH, sep = ""),
paste("indL", bulkSizeL, sep = ""),
popType, sep = "_"),
paste("Depth", minDepth, maxDepth, sep = "_"),
paste("Rep", repN, sep = "_"),
"pdf", sep = "."),
width = 6.5, height = 6)
ggsave(P_CI, filename = paste("QTLseqCI",
paste(paste("indH", bulkSizeH, sep = ""),
paste("indL", bulkSizeL, sep = ""),
popType, sep = "_"),
paste("Depth", minDepth, maxDepth, sep = "_"),
paste("Rep", repN, sep = "_"),
"png", sep = "."),
width = 6.5, height = 6, units = "in", dpi = 500)

置信区间和覆盖深度的关系


参考文献

Takagi H, Abe A, Yoshida K, et al. QTL-seq: rapid mapping of quantitative trait loci in rice by whole genome resequencing of DNA from two bulked populations. Plant J. 2013;74(1):174-183. doi:10.1111/tpj.12105

图片与主题无关


QTL-seq系列 | 如何计算置信区间
https://laowang2023.cn/2023/04/22/20230422-getQTLseqCI/
作者
老王
发布于
2023年4月22日
许可协议