SGAT是一个免费开源的单基因分析工具，基于Linux系统实现自动化批量处理，能够快速准确的完成单基因和表型的关联分析，只需要输入基因型和表型原始数据，即可计算出显著关联的SNP位点，并自动生成结果报告。
前段时间陆续的分享了SGAT（Single Gene Analysis Tool）的相关介绍，今天做一个总结整理，该工具是一个基于R语言tidyverse开发的快速分析流程化小工具，还存在很多的不足之处，欢迎大家多多指导。
接下来，将用8000字长文详解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 |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 |sh Archive: SGAT-V1.1.0.zip 1090a66274055c0b2cc578a43f0a4bce083ede4b Good finished！
依赖软件检查与安装
运行$ Rscript 1_check.R进行检查，根据提示安装相应软件和R包，直到所有依赖软件安装完成后提示finished，该过程也会自动检查基因型文件和表型文件，并对其进行提取，输出为列表，用于后续迭代计算。
###################### 单基因关联分析 ########################### Design by Jewel 使用方法： 1.将所有的基因型文件放在02文件夹中 例如"GeneID1A01G0123456.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文件中的染色体编号格式
- 两表关联后迭代提取匹配的观测值
- 基因型和表型文件整理
