GWAS结果显著SNP位点归类提取与变异类型转化
根据GWAS得到的Rresult文件信息，能够找出每个snp位点对应的显著性情况和基因变异信息，接下来，需要根据表格中的信息进行归纳总结，对不同显著性层次进行区分，找出可能性最大的点，过程比较繁琐。
这里笔者分享一个算法，使统计SNP和变异类型变的更加简便快捷，主要基于R语言的tidyverse完成。
主要步骤与思路解析
加载R包与环境，表型和基因列表文件定义变异信息转换函数创建输出数据框，包括基因和注释信息迭代筛选符合要求的SNP按照多个层次依次统计显著情况结果合并与注释项目运行环境
centos7 linuxR4.2.3操作步骤
加载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_genecolnames(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_infobreak}}}
上述算法的核心是先从基因列表中取一个基因，然后找这个基因对应的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再进行深入分析，本文不做延伸。
本文由mdnice多平台发布
如果觉得《R语言实现GWAS结果显著SNP位点归类提取与变异类型转化》对你有帮助，请点赞、收藏，并留下你的观点哦！
