总结丨SGAT单基因关联分析工具，一文上手使用

SGAT是一个免费开源的单基因分析工具，基于Linux系统实现自动化批量处理，能够快速准确的完成单基因和表型的关联分析，只需要输入基因型和表型原始数据，即可计算出显著关联的SNP位点，并自动生成结果报告。

前段时间陆续的分享了SGAT（Single Gene Analysis Tool）的相关介绍，今天做一个总结整理，该工具是一个基于R语言tidyverse开发的快速分析流程化小工具，还存在很多的不足之处，欢迎大家多多指导。

接下来，将用5000字长文详解SGAT的使用方法和算法原理，既是一个分享的过程，也是一个学习的过程。

背景信息

什么是单基因关联分析？

单基因关联分析是一种遗传学和生物统计学方法，用于研究基因与特定表型之间的关系。在单基因关联分析中，通常比较来自不同群体的不同等位基因频率。

如果某个等位基因在处理组中出现的频率显著高于对照组，则可以认为该等位基因与特定表型相关联。

单基因关联分析具有广泛应用，在医学、农业、动植物遗传学等领域都得到了广泛的应用！

待解决的问题

传统方式人工进行单基因关联分析需要从VCF文件开始，修改基因型文件，经过plink和taseel等软件转换文件格式，并手动修改变异信息，整理表型和基因型并互相匹配，逐步进行GWAS分析并根据结果作图，整个过程费时费力，而且极易出错。

因此，基于以上问题，开发了SGAT自动化单基因关联分析工具，快速完成多个基因多个表型多个模型的关联分析。

核心功能

变异信息自动识别与替换
染色体编号转换
基因型文件转换
表型与基因型匹配筛选
批量进行多模型GWAS分析
连锁不平衡作图
GWAS结果汇总整理
自动筛选显著性位点并提取变异信息
基因变异注释转换

定制化开发

GWAS分析模型自由选择
区间长度自由选择
筛选阈值自由选择
结果图片类型自由选择

源码开放性

 Mar 29 22:55 0_README.mdMar 22 20:25 1_check.RMar 19 21:40 2_gene_vcf2txt.RMar 22 20:12 3_hmp_trait_formate.RMar 20 11:05 4_GWAS_gapit.RMar 23 20:29 5_GWAS_results_translate.RMar 29 22:43 6_GWAS_Ttest_Result.RMar 22 20:14 clearn.shMar 31 11:53 start.sh

上述所有源码均在Github存放，其中bash脚本clearn.sh的功能是初始化工作目录并清除临时数据，start.sh的功能是启动自动化进程。

安装与部署运行环境

官网渠道（推荐）

curl https://www.jewin.love/install.sh |sh

Github仓库

git clone https://github.com/JewinZao/SGAT.git

本地安装

wget https://www.jewin.love/SGAT-V1.1.0.zip
unzip SGAT-V1.1.0.zip

通过上述方式安装SGAT工具，安装完成后可以在当前目录下看到脚本文件即成功！

$ curl https://www.jewin.love/install.sh |shArchive:  SGAT-V1.1.0.zip
1090a66274055c0b2cc578a43f0a4bce083ede4bGood finished！

依赖软件检查与安装

运行$ Rscript 1_check.R进行检查，根据提示安装相应软件和R包，直到所有依赖软件安装完成后提示finished，该过程也会自动检查基因型文件和表型文件，并对其进行提取，输出为列表，用于后续迭代计算。

###################### 单基因关联分析 ###########################Design by Jewel                           使用方法：                                                   1.将所有的基因型文件放在02文件夹中                           例如"TraesCS1A01G0123456.filter.vcf.gz"                    2.将表型文件放在05文件夹中，命名为trait.txt                  第一列名称为ID，后面每一列代表一个表型，例如"HT32L"        3.软件自动识别基因与表型信息                                 4.在当前文件夹下执行". ./start.sh"                           5.结果将在后续生成                                           6.初始化与清除工作空间请执行". ./clearn.sh"                  【 版本：V1.3.0 】                         #################################################################

方法：vcf转txt并自动规范化

vcf文件是存放基因变异信息的一种方式，本文提供一种算法，用于读取vcf文件并转换等位基因展示方法、替换染色体展示格式、以及自动识别非唯一变异并进行修改，用于对变异信息进行整理。

主要步骤与设计思路

读取VCF文件并分为三部分储存
提取变异信息并批量替换
修改染色体格式
SNP位点的判断与校正
单点碱基差异唯一化

具体操作步骤

加载R包与数据

library(tidyverse)
library(vcfR)
library(do)
library(R.utils)
df <- read.table(paste0("02_ordata/",job,".filter.vcf"),header = F)
vcf <- read.vcfR(paste0("02_ordata/",job,".filter.vcf.gz"))
chr_ref <- read.table("01_scripts/chr_num2str.txt",header = T)

读取VCF文件信息

fix <- vcf@fix
gt <- vcf@gt
meta <- vcf@meta

利用vcfR包读取入VCF文件后，分别提取出不同部分存放于临时变量中，以供后续使用。

批量替换变异信息

### 批量替换“|”为“/” ==================================================================
df[df == "0|0"] = "0/0"
df[df == "1|0"] = "1/0"
df[df == "0|1"] = "0/1"
df[df == "1|1"] = "1/1"
colnames(df) <- c(colnames(fix),colnames(gt))

该步骤的目的是为了将|修改为/,这是后面转hmp格式所需的条件。

替换染色体编号

###  替换染色体 =====================================================================
for (i in 1:nrow(df)){old_chr <- df$CHROM[i]for (k in 1:nrow(chr_ref)){if (chr_ref$chr_str[k] == old_chr){new_chr <- chr_ref$chr_num[k]df$CHROM[i] <- new_chr}}
}

利用for循环查找逐一取出染色体元素值，然后从参考信息中查找对应的正确格式，然后赋值给染色体信息，这一步中使用的chr_ref是染色体不同格式的对应信息。

参数识别与矫正

因为有插入缺失的存在，所以参考位置和实际位置的碱基并非完全唯一且差异，这将导致后面运行出错。这里提供一个算法，批量实现对SNP位点的检测与矫正。

snp_reverse函数

snp_reverse <- function(one,more){# 输入俩参，一为单二为多，返回存在于多但不与单同之值list_snp <- str_split(more,"")for (i in 1:str_length(more)){snp_now <- list_snp[[1]][i]ifelse(one==snp_now,next,return(snp_now))}
}

该函数输入两个参数，如“A,CATG”，首先将第二个参数分割成单个字母，然后迭代判断第一个字母是否与第二个一致，一旦出现与第一个参数不相同的值则返回该值。目的是为了让两个值长度为1且不相同。

批量处理ALT和REF位点

# 对每行的REF和ALT进行处理，将其变成不同值
for (i in 1:nrow(df)){ref <- df$REF[i]alt <- df$ALT[i]# 情况有三，均为单或其一为多if (str_length(ref) == 1){if (str_length(alt) == 1){}else{df$ALT[i] <- snp_reverse(ref,alt)}}else{if (str_length(alt) == 1){df$REF[i] <- snp_reverse(alt,ref)}else{print(paste0("ERROR：",df$ID[i]," this snp has more REF、ALT !"))}}
}

结果保存与输出

colnames(df)[1] <- "#CHROM"
write.table(df,paste0("03_vcf2txt/","gene_",job,".txt"),sep = "\t",row.names = F,col.names = T,quote = F)
print(paste0(job," Step ordata gene vcf to txt finished!"))

通过该算法能够对vcf文件进行转换，并得到规范化的txt文件，用于后续的分析。

方法：hmp文件与表型匹配

分析过程中，如果已经得到了hmp文件，下一步是将表型数据与hmp中的基因型数据一一对应，保证两者的样品ID信息一致，还需要对数据的格式进行规范化处理，用于后续的GWAS分析。

在此提供一种算法，能够实现对hmp文件和表型数据的关联筛选与校正。

主要步骤与设计思路

读取hmp文件和表型数据
替换hmp文件中的染色体编号格式
两表关联后迭代提取匹配的观测值
基因型和表型文件整理

具体操作步骤

加载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)

显著SNP位点提取与转化

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

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

主要步骤与思路解析

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

操作步骤

加载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_infobreak}}
}

上述算法的核心是先从基因列表中取一个基因，然后找这个基因对应的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再进行深入分析，本文不做延伸。

本项目测试运行环境

centos7 linux
R4.2.3

参考资料：

Plink、Tassel、LDBlockshow、GAPIT、Tidyverse、vcfR、ape、do、multtest、LDheatmap、genetics、scatterplot3d、EMMREML等

声明：

SGAT遵循国际GNU General Public License v3.0，核心算法和代码均开源公布，进行科学研究学习交流，不涉及商业使用，如果有任何问题欢迎联系。

软件公开发布链接：

doi.org/10.5281/zenodo.7783891

感谢您能看到这里，觉得有趣欢迎转发~

本文由mdnice多平台发布