GWAS结果显著SNP位点归类提取与变异类型转化

根据GWAS得到的Rresult文件信息，能够找出每个snp位点对应的显著性情况和基因变异信息，接下来，需要根据表格中的信息进行归纳总结，对不同显著性层次进行区分，找出可能性最大的点，过程比较繁琐。

这里笔者分享一个算法，使统计SNP和变异类型变的更加简便快捷，主要基于R语言的tidyverse完成。

主要步骤与思路解析

• 加载R包与环境，表型和基因列表文件
• 定义变异信息转换函数
• 创建输出数据框，包括基因和注释信息
• 迭代筛选符合要求的SNP
• 按照多个层次依次统计显著情况
• 结果合并与注释

项目运行环境

• centos7 linux
• R4.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_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再进行深入分析，本文不做延伸。