R语言实现GWAS结果显著SNP位点归类提取与变异类型转化
GWAS结果显著SNP位点归类提取与变异类型转化
根据GWAS得到的Rresult文件信息，能够找出每个snp位点对应的显著性情况和基因变异信息，接下来，需要根据表格中的信息进行归纳总结，对不同显著性层次进行区分，找出可能性最大的点，过程比较繁琐。
这里笔者分享一个算法，使统计SNP和变异类型变的更加简便快捷，主要基于R语言的tidyverse完成。
主要步骤与思路解析
- 加载R包与环境，表型和基因列表文件
- 定义变异信息转换函数
- 创建输出数据框，包括基因和注释信息
- 迭代筛选符合要求的SNP
- 按照多个层次依次统计显著情况
- 结果合并与注释
项目运行环境
- centos7 linux
- R4.2.3
操作步骤
加载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再进行深入分析，本文不做延伸。
本文由mdnice多平台发布
