批量查找SNP位点连锁区内对应的QTL以及基因

如果已知SNP位点的物理位置和其LDblock区间的端点，想要快速找到该区间内的QTL，之后根据参考基因组找到与连锁区域存在交集的基因，最终得到与SNP和QTL相匹配的基因集

通常的做法是在Excel中先对每个SNP计算出相应区间，然后找到对应的QTL，然后打开全部基因的参考信息，寻找有关基因。

上述方法较复杂，如果有成千上万条SNP或者基因将会非常耗费时间和精力，而且操作过程中容易出错，以下介绍一种快速实现的算法。

运行方法

• 安装R4.2.3，并加载tidyverse
• 将SNP和基因数据保存在data目录
• 运行run.R，结果将输出在out目录

算法原理

总体的思路是先对SNP进行区间进行合并，由于一个QTL可能对应多个SNP，因此要将这些位点扩展为一个区间，保证每个QTL都独立映射一个区间范围。然后从参考基因组找到最优的匹配项，与QTL合并保存。

数据读取

### 本算法用于批量查找snp位点对应的LDblock区间内QTL和基因 ============================================
library(tidyverse)
## QTL和SNP信息文件
f <- "./data/xxx.csv"
df <- read_csv(f)
## gene参考文件
f_ref <- "./data/ref.txt"
df_ref <- read_tsv(f_ref)

SNP和QTL数据整理

df <- df %>% arrange(QTL)
add_list <- list()
for (i in 1:nrow(df)){
  #取出一行进行添加信息
  snp <- df$target_SNP[i]
  atom <- str_split(snp,"_")
  snp_chr <- str_sub(atom[[1]][1],-2,-1)
  snp_loc <- as.numeric(atom[[1]][2])
  loc_L <- df$L_region[i]
  loc_R <- df$R_region[i]
  name_QTL <- df$QTL[i]
   
  # i=1单独处理
  if (i==1){
    add_list[[`name_QTL`]][1] <- name_QTL
    add_list[[`name_QTL`]][4] <- loc_L
    add_list[[`name_QTL`]][3] <- loc_R
    add_list[[`name_QTL`]][2] <- snp_chr
    next
  }
   
  # 存在多个连续QTL时，扩展右侧端点进行合并
  if (df$QTL[i]!=df$QTL[i-1]){
    add_list[[`name_QTL`]][1] <- name_QTL
    add_list[[`name_QTL`]][4] <- loc_L
    add_list[[`name_QTL`]][3] <- loc_R
    add_list[[`name_QTL`]][2] <- snp_chr
  }else{
    add_list[[`name_QTL`]][3] <- loc_R
  }
}

上面的算法对每个SNP进行迭代，将SNP切分成染色体和物理位置元素，然后记录左右区间端点，对于连续存在的QTL将右侧端点处的值进行扩展，实现了每个QTL对应位置区间的整理。

筛选QTL区间内的基因

首先创建一个空数据框，第一列是QTL名称，第二列是染色体位置，第三列是提取到的基因，并对参考基因组进行过滤，剔除NA值，用于后续迭代搜索。

out <- matrix(ncol = 3)
colnames(out) <- c("QTL","chr","gene")
out <- as.data.frame(out)
out <- out[-1,]
df_ref <- df_ref %>% drop_na()

然后对每个QTL进行搜索，在保证染色体一致的情况下，假如某基因的起始和终止位点都小于QTL区间左端点，代表该基因在QTL上游，如果都大于QTL区间右端点，则在下游。否则出现三种情况：基因在QTL区域内部、QTL在基因区域内部、QTL区间与基因区间有交集，这三种情况的基因就是目标基因集。

### 筛选处于QTL区间内的基因 ================================================================
for (QTL in add_list){
  print(QTL) 
  for (i in 1:nrow(df_ref)){
    if (df_ref$chr[i]==QTL[2]){
      if (df_ref$start[i]<QTL[3]|df_ref$end[i]<QTL[3]){
        next
      }else{
        if (df_ref$start[i]>QTL[4]|df_ref$end[i]>QTL[4]){
          next
        }else{
          tem <- c(QTL[1],df_ref$chr[i],df_ref$gene[i])
          out <- rbind(out,tem)
        }
      }
    }
  }
}

结果保存

write.csv(out,"./out/QTL_chr_gene.csv",quote = F,row.names = F)