失眠网 > R语言算法丨批量查找SNP位点连锁区内对应的QTL以及基因

R语言算法丨批量查找SNP位点连锁区内对应的QTL以及基因

时间：2022-07-13 23:44:59

批量查找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_QTLadd_list[[`name_QTL`]][4] <- loc_Ladd_list[[`name_QTL`]][3] <- loc_Radd_list[[`name_QTL`]][2] <- snp_chrnext}# 存在多个连续QTL时，扩展右侧端点进行合并if (df$QTL[i]!=df$QTL[i-1]){add_list[[`name_QTL`]][1] <- name_QTLadd_list[[`name_QTL`]][4] <- loc_Ladd_list[[`name_QTL`]][3] <- loc_Radd_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)}}}}}