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，欢迎转发支持~
END
© 素材来源于网络，侵权请联系后台删除
笔记合集，点击直达 
