批量查找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 in1: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 in1: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)
END
© 素材来源于网络，侵权请联系后台删除
