显著SNP位点提取与转化
根据GWAS得到的Rresult文件信息，能够找出每个snp位点对应的显著性情况和基因变异信息，接下来，需要根据表格中的信息进行归纳总结，对不同显著性层次进行区分，找出可能性最大的点，过程比较繁琐。
这里笔者分享一个算法，使统计SNP和变异类型变的更加简便快捷，主要基于R语言的tidyverse完成。
主要步骤与思路解析
- 加载R包与环境，表型和基因列表文件
- 定义变异信息转换函数
- 创建输出数据框，包括基因和注释信息
- 迭代筛选符合要求的SNP
- 按照多个层次依次统计显著情况
- 结果合并与注释
操作步骤
加载R包
library(tidyverse) library(writexl) library(xlsx)
读取输入文件
list_phe <- read.table("./01_scripts/list_phe.txt",header = F) # list_gene <- read.table("./01_scripts/list_gene.txt",header = F) list_gene <- read.table("./17_GWAS_SNP_varient_find/gene.id",header = F) varient_db <- read.table("./01_scripts/function/varient_name.txt",sep = "\t",header = F)
主要依赖三个文件，phe为变形列表，需要与GWAS结果的phe一致，gene为基因ID列表，varient_db是变异类型注释库，包含一一对应的变异信息。
变异信息转换
# 定义一个转换变异的函数 varient_name <- function(x){ if (x %in% varient_db$V1){ for (i in 1:nrow(varient_db)){ if (varient_db$V1[i]==x){ return(varient_db$V2[i]) } } }else{ return(x) } }
这里定义一个函数，对输入的变异类型自动查找匹配的注释信息，若出现不存在于已有的变异类型，则返回原始值，后续结果中方便检查和校正。
创建输出数据框
out <- list_gene colnames(out) <- "gene" out$additon <- NA
在计算开始前，创建一个空数据框，用于迭代过程中添加信息，提前分配储存空间，其中第一列为基因ID，第二列为注释。
迭代筛选算法
下面我提供了两种思路，方法一是先对每个表型下的所有snp进行判断，如果存在大于阈值的显著位点则备注，反之舍弃。
方法二是先找出单个SNP，然后再判断该位点处有多少个表型符合要求，如果存在多个表型均显著，则将其归纳统计到一起。
for (job in list_gene$V1){ print(job) df <- read.xlsx(paste0("./16_out_GWAS_and_T/",job,"_all.xlsx"),sheetIndex = 1) # 法一：寻找每个表型下的SNP # 7 9 11 13 15 17 19 21 23 25 27 29 为待提取的值 # for (i in seq(7,29,2)){ # phe <- colnames(df)[i] # df_p7_snp <- df %>% arrange(!!sym(phe)) %>% filter(!!sym(phe)>7) # df_p3_snp <- df %>% arrange(!!sym(phe)) %>% filter(!!sym(phe)>3) %>% filter(!!sym(phe)<7) # # P值大于7 # var_en <- df_p7_snp$T_eff[1] %>% str_split("[,]") %>% str_split("[|]") # var_en <- var_en[[1]][2] # var_cn <- varient_name(var_en) # } # 法二：寻找每个snp下符合的表型 find <- matrix(ncol = 4,nrow = 0) colnames(find) <- c("snp","var","p","phe") for (i in 1:nrow(df)){ snp_name <- df$SNP[i] if (is.na(df$T_eff[i])){next} snp_var_en <- df$T_eff[i] %>% str_split("[,]") snp_var_en <- snp_var_en[[1]][1] %>% str_split("[|]") if (substr(snp_var_en,4,22)!=job){next} snp_var_en <- snp_var_en[[1]][2] snp_var_en <- varient_name(snp_var_en) snp_phe_p <- df[i,c(seq(7,29,2))] find_phe <- c() for (i in 1:ncol(snp_phe_p)){ if (snp_phe_p[1,i]>7){ find_phe <- c(find_phe,colnames(snp_phe_p)[i]) } } find_snp <- c(snp_name,snp_var_en,"[P>7]",paste0(find_phe,collapse = "+")) if (find_snp[4]!=""){ find <- rbind(find,find_snp) } } if (nrow(find) == 0){ find <- matrix(ncol = 4,nrow = 0) colnames(find) <- c("snp","var","p","phe") for (i in 1:nrow(df)){ snp_name <- df$SNP[i] if (is.na(df$T_eff[i])){next} snp_var_en <- df$T_eff[i] %>% str_split("[,]") snp_var_en <- snp_var_en[[1]][1] %>% str_split("[|]") if (substr(snp_var_en,4,22)!=job){next} snp_var_en <- snp_var_en[[1]][2] snp_var_en <- varient_name(snp_var_en) snp_phe_p <- df[i,c(seq(7,29,2))] find_phe <- c() for (i in 1:ncol(snp_phe_p)){ if (snp_phe_p[1,i]>5){ find_phe <- c(find_phe,colnames(snp_phe_p)[i]) } } find_snp <- c(snp_name,snp_var_en,"[P>5]",paste0(find_phe,collapse = "+")) if (find_snp[4]!=""){ find <- rbind(find,find_snp) } } } if (nrow(find) == 0){ find <- matrix(ncol = 4,nrow = 0) colnames(find) <- c("snp","var","p","phe") for (i in 1:nrow(df)){ snp_name <- df$SNP[i] if (is.na(df$T_eff[i])){next} snp_var_en <- df$T_eff[i] %>% str_split("[,]") snp_var_en <- snp_var_en[[1]][1] %>% str_split("[|]") if (substr(snp_var_en,4,22)!=job){next} snp_var_en <- snp_var_en[[1]][2] snp_var_en <- varient_name(snp_var_en) snp_phe_p <- df[i,c(seq(7,29,2))] find_phe <- c() for (i in 1:ncol(snp_phe_p)){ if (snp_phe_p[1,i]>3){ find_phe <- c(find_phe,colnames(snp_phe_p)[i]) } } find_snp <- c(snp_name,snp_var_en,"[P>3]",paste0(find_phe,collapse = "+")) if (find_snp[4]!=""){ find <- rbind(find,find_snp) } } } var_info <- c() out_info <- c() if (nrow(find)==0){ out_info <- "GAPIT：log10.P < 3" }else{ for (i in 1:nrow(find)){ var_info <- c(var_info,find[i,2],find[i,1],find[i,3],paste0("(",find[i,4],"),")) } out_info <- paste0(nrow(find),"个-GAPIT分析",paste0(var_info,collapse ="")) out_info <- substr(out_info,1,nchar(out_info)-1) } for (i in 1:nrow(out)){ if (identical(out$gene[i],job)){ out$additon[i] <- out_info break } } }
上述算法的核心是先从基因列表中取一个基因，然后找这个基因对应的snp和表型，如果找到某些snp在多个表型中显著性都大于7，则将其添加到注释信息，但是如果没有大于7的位点，则开始继续寻找是否存在大于5的位点，以此类推，若也没有大于5的点，则寻找大于3的位点。
该过程将显著区间分为三层，只有上层个数为零时，才会启动下一层的搜索，因此保证了每次结果的显著性差异保持在相对较平均的范围中，防止过大过小的位点同时选中。
结果保存
write.xlsx(out, "./17_GWAS_SNP_varient_find/gene_infomation.xlsx", sheetName = "varient", row.names = F,col.names = T)
结果文件保存在out变量中，将其输出为excel即可，如有其它想法可以根据out再进行深入分析，本文不做延伸。
本项目测试运行环境
- centos7 linux
- R4.2.3
参考资料：
Plink、Tassel、LDBlockshow、GAPIT、Tidyverse、vcfR、ape、do、multtest、LDheatmap、genetics、scatterplot3d、EMMREML等
声明：
SGAT遵循国际GNU General Public License v3.0，核心算法和代码均开源公布，进行科学研究学习交流，不涉及商业使用，如果有任何问题欢迎联系。
软件公开发布链接：
doi.org/10.5281/zenodo.7783891
感谢您能看到这里，觉得有趣欢迎转发~
END
© 素材来源于网络，侵权请联系后台删除
笔记合集，点击直达 
