具体操作步骤
加载R包与数据
library(tidyverse) chr_ref <- read.table("01_scripts/chr_num2str.txt",header = T) df <- read_table(paste0("04_hmp/gene_",job,".hmp.txt"),show_col_types = F) trait <- read_table(paste0("05_trait/","trait.txt"),show_col_types = F)
读取三个数据文件，其中第一个是染色体ID个不同格式对应信息，第二个是基因型hmp.txt文件，第三个是表型数据文件。
染色体格式转换
- chr_id_translate 函数
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{ print("Please input again! type inaviably") } } } }
该函数提供了一种对染色体格式的快速转换方法，可以对数字型、字符型、全称之间进行快速转换，第一个参数是原始的编号，第二个参数选择转换方式，返回值是一个新的染色体编码值。
- 批量替换
for (i in 1:nrow(df)){ df$chrom[i] <- chr_id_translate( df$chrom[i],type = "1_to_1A") }
通过迭代将所有的数值型染色体编号换成数字加字母型。
基因型和表型匹配筛选
- 数据转换与处理
df2 <- rbind(colnames(df),df) df_gene <- t(df2) df_add_gene <- matrix(ncol = ncol(df_gene)) df_add_gene <- df_add_gene[-1,] df_add_trait <- matrix(ncol = ncol(trait)) df_add_trait <- df_add_trait[-1,] df_gene <- as.data.frame(df_gene)
对原始数据进行转置，目的是为了让基因型中样品ID按行排布，方便后续筛选，定义一个新的数据框用于储存迭代输出信息。
- 迭代提取匹配观测值
for (i in 1:nrow(df_gene)){ id_gene <- df_gene$V1[i] for (k in 1:nrow(trait)){ id_trait <- trait$ID[k] if (id_gene == id_trait){ my_gene <- df_gene[i,] my_trait <- trait[k,] df_add_gene <- rbind(df_add_gene,my_gene) df_add_trait <- rbind(df_add_trait,my_trait) }else{ next } } }
通过上述方法可以找出两个表格中完全匹配的样品，生成的df_add_gene是所有匹配到的基因型文件，df_add_trait是所有对应的表型文件。后续可以直接拿来做GAPIT分析。
结果输出与保存
out_gene <- rbind(df_gene[1:11,],df_add_gene) out_genet <- t(out_gene) gene_final <- as.data.frame(out_genet) write.table(gene_final,paste0("./06_out_gene/",job,".gene.hmp.txt"), quote = F,sep = "\t",col.names = F,row.names = F) trait_final <- as.data.frame(df_add_trait) write.table(trait_final,paste0("./07_out_trait/",job,".trait.txt"), quote = F,sep = "\t",col.names = T,row.names = F) print(paste0(job," hmp and trait formate finished!"))
重新合并头文件并转置，恢复原有结构，然后分别将两个结果保存到对应文件夹中。
方法：GAPIT进行GWAS分析
GAPIT是张志武老师开发的基于R语言的GWAS分析工具，能够根据表型和基因型数据自动进行不同模型的全基因组关联分析，网上有很多公开的教程。本文分析一种方法，进行单基因GWAS分析。
主要步骤
- 加载分析环境
- 导入数据
- 选择模型并开始分析
- 结果提取
具体操作步骤
加载R包与环境
library(MASS) # required for ginv library(multtest) library(gplots) library(compiler) #required for cmpfun library(scatterplot3d) library(bigmemory) library(ape) library(EMMREML) source("./01_scripts/GAPIT1.txt") source("./01_scripts/GAPIT2.txt")
导入数据
myG <- read.delim(paste0("./06_out_gene/",job,".gene.hmp.txt"), header = F) myY <- read.table(paste0("./07_out_trait/",job,".trait.txt"), header = T,sep = "\t")
这里需要的数据有两个，myG是基因型文件，需要hmp格式，myY是表型文件，需要制表符分隔的txt文件。
设置项目路径
now_dir <- getwd() dir.create(paste0(now_dir,"/08_out_GWAS/MLM_",job)) setwd(paste0(now_dir,"/08_out_GWAS/MLM_",job))
由于GAPIT运行后会自动在当前目录下生成若干结果文件，为了避免紊乱，因此对每次结果设置独立路径。这里会读取当前文件夹，然后创建新文件夹并设为临时工作目录。
GAPIT分析
myGAPIT <- GAPIT( Y=myY, G=myG, PCA.total=3, model="MLM", Random.model = TRUE )
该步骤是GWAS的核心步骤，根据样本数据量的不同，这一步耗费的时间也不同，完成后会看到很多自动生成的图片和表格文件，该步骤可以选择不同的模型，比如MLM等。
setwd(now_dir) print(paste0(job," GWAS finished!"))
完成后重新回到之前的工作目录
方法：GWAS结果整理
在使用GAPIT进行GWAS分析后，会自动在工作目录下生成若干结果文件，其中相对比较重要的是result.csv文件，该文件中展示了得到的显著位点详细信息，比如染色体、物理位置、p值等，接下来介绍一种算法，对其进行整理计算为绘图所需格式。
主要步骤与思路
- 读取数据文件
GWAS.Results.csv - 替换染色体格式
- 计算上下游区域
- 计算region信息
- 生成结果文件
具体操作步骤
加载环境和数据
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)
