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.md Mar 22 20:25 1_check.R Mar 19 21:40 2_gene_vcf2txt.R Mar 22 20:12 3_hmp_trait_formate.R Mar 20 11:05 4_GWAS_gapit.R Mar 23 20:29 5_GWAS_results_translate.R Mar 29 22:43 6_GWAS_Ttest_Result.R Mar 22 20:14 clearn.sh Mar 31 11:53 start.sh上述所有源码均在Github寄存，其中bash脚本clearn.sh的性能是初始化工作目录并革除长期数据，start.sh的性能是启动自动化过程。
装置与部署运行环境
官网渠道（举荐）
curl https://www.jewin.love/install.sh |shGithub仓库
git clone https://github.com/JewinZao/SGAT.git本地装置
wget https://www.jewin.love/SGAT-V1.1.0.zipunzip SGAT-V1.1.0.zip通过上述形式装置SGAT工具，装置实现后能够在当前目录下看到脚本文件即胜利！
$ curl https://www.jewin.love/install.sh |shArchive: SGAT-V1.1.0.zip1090a66274055c0b2cc578a43f0a4bce083ede4bGood 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@fixgt <- vcf@gtmeta <- 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 ginvlibrary(multtest)library(gplots)library(compiler) #required for cmpfunlibrary(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 - 500s_2 <- s_2 + 500region <- 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_genecolnames(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_info break } }}上述算法的外围是先从基因列表中取一个基因，而后找这个基因对应的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多平台公布