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

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

计算原理

模拟取样过程

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

模拟计算 ΔSNP-index

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

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

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

代码实现

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

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

# 加载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 取样多少的关系。

# 将置信区间和覆盖深度的关系可视化
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)