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多平台发布