共计 13644 个字符，预计需要花费 35 分钟才能阅读完成。
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 <- 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_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 多平台公布