GWAS结果整理丨利用R语言tidyverse自动统计显著位点

GWAS结果文件分析与处理方法

引言

在使用GAPIT进行GWAS分析后，会自动在工作目录下生成若干结果文件，其中相对比较重要的是result.csv文件，该文件中展示了得到的显著位点详细信息，比如染色体、物理位置、p值等，接下来介绍一种算法，对其进行整理计算为绘图所需格式。

---

主要步骤与思路

• 读取数据文件GWAS.Results.csv
• 替换染色体格式
• 计算上下游区域
• 计算region信息
• 生成结果文件

项目运行环境

• centos7 linux
• R4.2.3

具体操作步骤

加载环境和数据

rm(list = ls())
library(tidyverse)
 
ARGS <- commandArgs(T)
print(paste0("Results Working Gene ID:",ARGS[1]))
job <- ARGS[1]
dir_MLM <- paste0("MLM_",job)
phe <- ARGS[2]
file_name <- paste0("/GAPIT.MLM.",phe,".GWAS.Results.csv")
df <- read.csv(paste0("./08_out_GWAS/",dir_MLM,file_name),header = T)

主要实用tidyverse包进行数据处理，ARGS是脚本的参数设置，如果单个任务可以直接读入文件，不用脚本传参，只需要设置好文件名进行读取。

染色体格式转换

###  替换染色体展示方式，1A_to_1 ===========================================================
chr_ref <- read.table("01_scripts/chr_num2str.txt",header = T)
# 读取染色体转换参考信息，可以进行自定义修改
chr_id_translate <- function(data,type){
  # 输入俩参，一为原始数据，二为类型
  if (type == "1_to_chr1A"){
    # 数字转字符型
    old_id <- as.character(data)
    for (k in 1:nrow(chr_ref)){
      if (as.character(chr_ref$chr_num[k]) == old_id){
        return(chr_ref$chr_str[k])
      }
    }
  }else{
    if (type == "chr1A_to_1"){
      # 字符转数字型
      old_id <- as.character(data)
      for (k in 1:nrow(chr_ref)){
        if (as.character(chr_ref$chr_str[k]) == old_id){
          return(chr_ref$chr_num[k])
        }
      }
    }else{
      if (type == "1_to_1A"){
        old_id <- as.character(data)
        for (k in 1:nrow(chr_ref)){
          if (as.character(chr_ref$chr_num[k]) == old_id){
            new <- paste0(chr_ref$atom7[k],chr_ref$atom3[k],sep="")
            return(new)
          }
        }
      }else{
        if (type == "1A_to_1"){
          old_id <- as.character(data)
          for (k in 1:nrow(chr_ref)){
            temp <- paste0(chr_ref$atom7[k],chr_ref$atom3[k],sep="")
            if (as.character(temp) == old_id){
              return(chr_ref$chr_num[k])
            }
          }
        }else{
        print("Please input again! type inaviably")
        }
      }
    }
  }
}

刚刚定义了一个函数chr_id_translate能够对染色体文件进行自定义转换，接下来将其依次应用到数据的染色体列。

for (i in 1:nrow(df)){
  df$Chromosome[i] <- chr_id_translate(df$Chromosome[i],"1A_to_1")
}

物理位置区间计算

根据Postion信息计算最大值和最小值，分别向上下游扩展500bp就能得到想要的区间，将其保存为region，用于后续绘图使用

s_1 <- min(df$Position)
s_2 <- max(df$Position)
s_1 <- s_1 - 500
s_2 <- s_2 + 500
region <- paste0(df$Chromosome[1],":",s_1,":",s_2)

结果保存

绘图需要三列信息，分别是染色体、物理位置、p值，因此将这部分数据单独存放到df_new，然后保存为新文件。

###  生成新文件，染色体-位置-P值 =============================================================
df_new <- df[,2:4]
file_new <- paste0("./09_out_MLM/",job,"_MLM.",phe,".GWAS.Results.csv",sep="")
write_csv(df_new,file_new,col_names=F)

至此，这个方法的原理已分享完毕，如果您在使用过程中有问题或者建议均可提交issues到Github，欢迎转发支持~

本文由mdnice多平台发布