以水稻为例教你如何使用BSA方法进行遗传定位（下篇）

往期回顾：

以水稻为例教你如何使用BSA方法进行遗传定位（上篇）

以水稻为例教你如何使用BSA方法进行遗传定位（中篇）

接着上一篇留下的问题，如何按照以1M为窗口，每次移动10Kb计算均值。我们以KY131是"0/0", 而 "DN422" 是 "1/1"的位点结果为例，我们来增加这条线

原本我的想法是使用split根据每个位点的位置信息进行拆分，但是这里有一个问题是，如果只是每个窗口没有重叠的话，那么每个位点应该都有一个唯一的因子，但这里是以10kb将长度为1Mb的窗口从左往右移动，也就意味着有很多位点 是 属于多个窗口的，也就是无法直接用split-apply-combine的方法进行处理。

于是乎，我写了一个函数用于根据窗口来计算均值，它的逻辑就是先创建好窗口，然后根据窗口的起始位置和终止位置去筛选符合要求的位点，最后进行求和。同时考虑有些位置可能没有SNP，那么最后还进行了一波过滤。

函数的输入是之前snp-index的位置和对应的值，以及确定窗口的大小和步长，

calcValueByWindow 
                              window_size = 1000000,
                              step_size = 100000){
  # get the max position in the postion
  max_pos 

  # construct the window
  window_start 
  window_end 
  mean_value 

  # select the value inside the window
  for (j in seq_along(window_start)){

    pos_in_window  window_start[j] &
                             pos < window_end[j])
    value_in_window 

    mean_value[j] 

  }
  # remove the Not A Number position
  nan_pos 
  mean_value 
  window_pos 
  df 
                   value = mean_value)
  return(df)
}

得到结果就可以用 lines 在之前结果上加上均值线

par(mfrow = c(3,4))

for (i in paste0("chr", formatC(1:12, width = 2, flag=0)) ){

  freq_flt   pos 7))
  snp_index 1] - freq_flt[,2]

  # bin
  df pos = pos, value = snp_index)

  plot(x = pos, y =snp_index,
       ylim = c(-1,1),
       pch = 20, cex = 0.2,
       xlab = i,
       ylab = expression(paste(Delta, " " ,"SNP index")))
  lines(x = df$pos, y = df$value, col = "red")
}

最后再对文章总结一下。文章并不是只用了BSA的方法进行定位，他们花了几年的时间用SSR分子标记确定了候选基因可能区间，用BSA的方法在原有基础上缩小了定位区间。当然即便如此，候选基因也有上百个，作者通过BLAST的方式，对这些基因进行了注释。尽管中间还加了一些GO富集分析的内容，说这些基因富集在某个词条里，有一个是DNA metabolic processes(GO:0006259)，但我觉得如果作者用clusterProfiler做富集分析，它肯定无法得到任何富集结果。他做富集分析的目的是其实下面这个描述，也就是找到和抗冻相关的基因

LOC_Os06g39740 and LOC_Os06g39750,were annotated as the function of “response to cold (GO: 0009409)”, suggesting their key roles in regulating cold tolerance in rice. "

当然他还做了qRT-PCR进行了验证，最后推测LOC_Os06g39750应该是目标基因，这个基因里还有8个SNP位点。

关于水稻的BSA分析的全文就这些了，完整版阅读原文即可直达。明天还会更新一篇如何用qtlseqr来完成下游分析。