干货丨 一文详解SGAT单基因关联分析工具(二)

简介: 干货丨 一文详解SGAT单基因关联分析工具(二)

具体操作步骤

加载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)
相关文章
|
算法 Linux Shell
干货丨 一文详解SGAT单基因关联分析工具
干货丨 一文详解SGAT单基因关联分析工具
|
算法 Linux Python
干货丨 一文详解SGAT单基因关联分析工具(三)
干货丨 一文详解SGAT单基因关联分析工具(三)
|
新零售 算法 数据挖掘
小白学数据分析之关联分析理论篇
关联分析的学习 在说关联分析之前,先说说自己这段时间的一些感受吧,这段时间相对轻松一些,有一些时间自己自己来学习一些新东西和知识,然而却发现捧着一本数据挖掘理论的书籍在一点一点的研读实在是很漫长,而且看过了没有什么感觉。
1229 0
|
存储 算法 数据挖掘
小白学数据分析之关联分析算法篇Apriori
早些时候写过关于购物篮分析的文章,其中提到了C5.0和Apriori算法,没有仔细说说这算法的含义,昨天写了一下关联分析的理论部分,今天说说关联分析算法之一的Apriori算法,很多时候大家都说,数据分析师更多的是会用就可以了,不必纠结于那些长篇累牍的理论,其实我觉得还是有点必要的,你未必要去设计算法,但是如果你掌握和熟知一个算法,这对于你如何驾驭和使用这个算法是很有帮助的,此外每个算法都有使用的局限性,比如空间和时间复杂度,使用条件约束。
1568 0
|
9月前
|
算法 数据挖掘 C++
九大数据分析方法-综合型分析方法以及如何使用这九大分析方法
九大数据分析方法-综合型分析方法以及如何使用这九大分析方法
《大数据分析原理与实践》——3.3 相关分析
本节书摘来自华章计算机《大数据分析原理与实践》一书中的第3章,第3.3节,作者 王宏志,更多章节内容可以访问云栖社区“华章计算机”公众号查看。
1394 0
《大数据分析原理与实践》一一3.3 相关分析
本节书摘来自华章出版社《大数据分析原理与实践》一 书中的第3章,第3.3节,作者:王宏志 ,更多章节内容可以访问云栖社区“华章计算机”公众号查看。
1342 0
|
算法 数据挖掘 数据库
|
存储 算法 数据可视化
利用TCseq包进行基因表达趋势分析
TCseq包提供了一个统一的套件去处理不同时序类型的数据分析,可以应用于转录组或者像ATAC-seq,Chip-seq的表观基因组时序型数据分析。该包主要的集中于不同时间点的差异分析,时间趋势分析及可视化作图。
525 0

热门文章

最新文章