GWAS結(jié)果顯著SNP位點歸類提取與變異類型轉(zhuǎn)化

根據(jù)GWAS得到的Rresult文件信息，能夠找出每個snp位點對應(yīng)的顯著性情況和基因變異信息，接下來，需要根據(jù)表格中的信息進行歸納總結(jié)，對不同顯著性層次進行區(qū)分，找出可能性最大的點，過程比較繁瑣。

這里筆者分享一個算法，使統(tǒng)計SNP和變異類型變的更加簡便快捷，主要基于R語言的tidyverse完成。

主要步驟與思路解析

• 加載R包與環(huán)境，表型和基因列表文件
• 定義變異信息轉(zhuǎn)換函數(shù)
• 創(chuàng)建輸出數(shù)據(jù)框，包括基因和注釋信息
• 迭代篩選符合要求的SNP
• 按照多個層次依次統(tǒng)計顯著情況
• 結(jié)果合并與注釋

項目運行環(huán)境

• 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結(jié)果的phe一致，gene為基因ID列表，varient_db是變異類型注釋庫，包含一一對應(yīng)的變異信息。

變異信息轉(zhuǎn)換

# 定義一個轉(zhuǎn)換變異的函數(shù)
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)
      }
}

這里定義一個函數(shù)，對輸入的變異類型自動查找匹配的注釋信息，若出現(xiàn)不存在于已有的變異類型，則返回原始值，后續(xù)結(jié)果中方便檢查和校正。

創(chuàng)建輸出數(shù)據(jù)框

out <- list_gene
colnames(out) <- "gene"
out$additon <- NA

在計算開始前，創(chuàng)建一個空數(shù)據(jù)框，用于迭代過程中添加信息，提前分配儲存空間，其中第一列為基因ID，第二列為注釋。

迭代篩選算法

下面我提供了兩種思路，方法一是先對每個表型下的所有snp進行判斷，如果存在大于閾值的顯著位點則備注，反之舍棄。方法二是先找出單個SNP，然后再判斷該位點處有多少個表型符合要求，如果存在多個表型均顯著，則將其歸納統(tǒng)計到一起。

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
            }
      }
}

上述算法的核心是先從基因列表中取一個基因，然后找這個基因?qū)?yīng)的snp和表型，如果找到某些snp在多個表型中顯著性都大于7，則將其添加到注釋信息，但是如果沒有大于7的位點，則開始繼續(xù)尋找是否存在大于5的位點，以此類推，若也沒有大于5的點，則尋找大于3的位點。

該過程將顯著區(qū)間分為三層，只有上層個數(shù)為零時，才會啟動下一層的搜索，因此保證了每次結(jié)果的顯著性差異保持在相對較平均的范圍中，防止過大過小的位點同時選中。

結(jié)果保存

write.xlsx(out,
    "./17_GWAS_SNP_varient_find/gene_infomation.xlsx",
    sheetName = "varient",
    row.names = F,col.names = T)

結(jié)果文件保存在out變量中，將其輸出為excel即可，如有其它想法可以根據(jù)out再進行深入分析，本文不做延伸。

本文由mdnice多平臺發(fā)布