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 |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文件夹中
例如"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 读取VCF文件信息
fix 利用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) 该步骤的目的是为了将|修改为/,这是后面转hmp格式所需的条件。
替换染色体编号
### 替换染色体 =====================================================================
for (i in 1:nrow(df)){
old_chr 利用for循环查找逐一取出染色体元素值，然后从参考信息中查找对应的正确格式，然后赋值给染色体信息，这一步中使用的chr_ref是染色体不同格式的对应信息。
参数识别与矫正
因为有插入缺失的存在，所以参考位置和实际位置的碱基并非完全唯一且差异，这将导致后面运行出错。这里提供一个算法，批量实现对SNP位点的检测与矫正。
- snp_reverse函数
snp_reverse 该函数输入两个参数，如“A,CATG”，首先将第二个参数分割成单个字母，然后迭代判断第一个字母是否与第二个一致，一旦出现与第一个参数不相同的值则返回该值。目的是为了让两个值长度为1且不相同。
批量处理ALT和REF位点
# 对每行的REF和ALT进行处理，将其变成不同值
for (i in 1:nrow(df)){
ref 结果保存与输出
colnames(df)[1] 通过该算法能够对vcf文件进行转换，并得到规范化的txt文件，用于后续的分析。
方法：hmp文件与表型匹配
分析过程中，如果已经得到了hmp文件，下一步是将表型数据与hmp中的基因型数据一一对应，保证两者的样品ID信息一致，还需要对数据的格式进行规范化处理，用于后续的GWAS分析。
在此提供一种算法，能够实现对hmp文件和表型数据的关联筛选与校正。
主要步骤与设计思路
- 读取hmp文件和表型数据
- 替换hmp文件中的染色体编号格式
- 两表关联后迭代提取匹配的观测值
- 基因型和表型文件整理
具体操作步骤
加载R包与数据
library(tidyverse)
chr_ref 读取三个数据文件，其中第一个是染色体ID个不同格式对应信息，第二个是基因型hmp.txt文件，第三个是表型数据文件。
染色体格式转换
- chr_id_translate 函数
chr_id_translate 该函数提供了一种对染色体格式的快速转换方法，可以对数字型、字符型、全称之间进行快速转换，第一个参数是原始的编号，第二个参数选择转换方式，返回值是一个新的染色体编码值。
- 批量替换
for (i in 1:nrow(df)){
df$chrom[i] 通过迭代将所有的数值型染色体编号换成数字加字母型。
基因型和表型匹配筛选
-
数据转换与处理
df2对原始数据进行转置，目的是为了让基因型中样品ID按行排布，方便后续筛选，定义一个新的数据框用于储存迭代输出信息。
-
迭代提取匹配观测值
for (i in 1:nrow(df_gene)){ id_gene通过上述方法可以找出两个表格中完全匹配的样品，生成的
df_add_gene是所有匹配到的基因型文件，df_add_trait是所有对应的表型文件。后续可以直接拿来做GAPIT分析。
结果输出与保存
out_gene 重新合并头文件并转置，恢复原有结构，然后分别将两个结果保存到对应文件夹中。
方法：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 这里需要的数据有两个，myG是基因型文件，需要hmp格式，myY是表型文件，需要制表符分隔的txt文件。
设置项目路径
now_dir 由于GAPIT运行后会自动在当前目录下生成若干结果文件，为了避免紊乱，因此对每次结果设置独立路径。这里会读取当前文件夹，然后创建新文件夹并设为临时工作目录。
GAPIT分析
myGAPIT 该步骤是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 主要实用tidyverse包进行数据处理，ARGS是脚本的参数设置，如果单个任务可以直接读入文件，不用脚本传参，只需要设置好文件名进行读取。
染色体格式转换
### 替换染色体展示方式，1A_to_1 ===========================================================
chr_ref 刚刚定义了一个函数chr_id_translate能够对染色体文件进行自定义转换，接下来将其依次应用到数据的染色体列。
for (i in 1:nrow(df)){
df$Chromosome[i] 物理位置区间计算
根据Postion信息计算最大值和最小值，分别向上下游扩展500bp就能得到想要的区间，将其保存为region，用于后续绘图使用
s_1 结果保存
绘图需要三列信息，分别是染色体、物理位置、p值，因此将这部分数据单独存放到df_new，然后保存为新文件。
### 生成新文件，染色体-位置-P值 =============================================================
df_new 显著SNP位点提取与转化
根据GWAS得到的Rresult文件信息，能够找出每个snp位点对应的显著性情况和基因变异信息，接下来，需要根据表格中的信息进行归纳总结，对不同显著性层次进行区分，找出可能性最大的点，过程比较繁琐。
这里笔者分享一个算法，使统计SNP和变异类型变的更加简便快捷，主要基于R语言的tidyverse完成。
主要步骤与思路解析
- 加载R包与环境，表型和基因列表文件
- 定义变异信息转换函数
- 创建输出数据框，包括基因和注释信息
- 迭代筛选符合要求的SNP
- 按照多个层次依次统计显著情况
- 结果合并与注释
操作步骤
加载R包
library(tidyverse)
library(writexl)
library(xlsx)读取输入文件
list_phe 主要依赖三个文件，phe为变形列表，需要与GWAS结果的phe一致，gene为基因ID列表，varient_db是变异类型注释库，包含一一对应的变异信息。
变异信息转换
# 定义一个转换变异的函数
varient_name 这里定义一个函数，对输入的变异类型自动查找匹配的注释信息，若出现不存在于已有的变异类型，则返回原始值，后续结果中方便检查和校正。
创建输出数据框
out 在计算开始前，创建一个空数据框，用于迭代过程中添加信息，提前分配储存空间，其中第一列为基因ID，第二列为注释。
迭代筛选算法
下面我提供了两种思路，方法一是先对每个表型下的所有snp进行判断，如果存在大于阈值的显著位点则备注，反之舍弃。
方法二是先找出单个SNP，然后再判断该位点处有多少个表型符合要求，如果存在多个表型均显著，则将其归纳统计到一起。
for (job in list_gene$V1){
print(job)
df % arrange(!!sym(phe)) %>% filter(!!sym(phe)>7)
# df_p3_snp % arrange(!!sym(phe)) %>% filter(!!sym(phe)>3) %>% filter(!!sym(phe)% str_split("[,]") %>% str_split("[|]")
# var_en % str_split("[,]")
snp_var_en % str_split("[|]")
if (substr(snp_var_en,4,22)!=job){next}
snp_var_en 7){
find_phe 7]",paste0(find_phe,collapse = "+"))
if (find_snp[4]!=""){
find % str_split("[,]")
snp_var_en % str_split("[|]")
if (substr(snp_var_en,4,22)!=job){next}
snp_var_en 5){
find_phe 5]",paste0(find_phe,collapse = "+"))
if (find_snp[4]!=""){
find % str_split("[,]")
snp_var_en % str_split("[|]")
if (substr(snp_var_en,4,22)!=job){next}
snp_var_en 3){
find_phe 3]",paste0(find_phe,collapse = "+"))
if (find_snp[4]!=""){
find 上述算法的核心是先从基因列表中取一个基因，然后找这个基因对应的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多平台发布
服务器托管，北京服务器托管，服务器租用 http://www.fwqtg.net
