GWAS 结果批量整理
今天分享一个算法，主要是利用R语言对GWAS分析得到的结果进行批量整理，适用于批量筛选关键SNP位点与对应的基因或QTL，尤其是多个表型多个材料的批量分析。
GWAS（全基因组关联分析）是一种常用的遗传学研究方法，用于探索基因与表型之间的关联。它通过对大规模样本集合进行基因组广泛扫描，寻找与特定表型（如疾病、性状等）相关的遗传变异。
GWAS结果示例
通过GWAS分析后会得到一个储存显著SNP的结果文件，假设其命名为“组别.表型.模型.阈值.result.txt”,理论上会有很多个这种格式的文件，文件内容如下：
"INDEX" "SNP" "CHROM" "POS" "REF" "ALT" "Effect" "SE" "trait" "1" "rs151515" 1 25554 "C" "A" -1.630526 0.39931444 6.688587e-06 "2" "rs32151" 1 72857 "T" "TAC" 0.734972 0.17574796 9.271553e-06
第二列是SNP的位点名称，最后一列是对应的P值，这两个信息非常关键。
整理思路与算法解释
首先，library(tidyverse)加载R包，代码流程均使用优雅的tidy框架。
library(tidyverse)
step1：原始数据整理
大致思路是先读取当前目录下的文件列表，然后依次循环执行计算过程，识别表型、模型、P值等参数，然后传递给筛选函数，并对符合条件的值进行标注，最后会将结果写出为一个csv文件，用于下一步的计算。
id_list <- list.files("./data/",pattern = "*.txt") for (id in id_list){ file_name <- paste0("./",id) atom <- str_split(id,"[.]") # Group A ----- if (length(atom[[1]]) == 5){ phe <- atom[[1]][1] way <- atom[[1]][2] %>% str_replace("Farm","") plast <- atom[[1]][3] # 特异性标注P值并将其装换为数字型 if (plast == "1e-5"){plast <- 6}else{plast <- as.numeric(plast)} print(file_name) # 计算P值并筛选 df <- read_delim(file_name,delim = " ", col_types = cols(CHROM = col_character())) colnames(df)[9] <- way df %>% mutate(log = round(-log10(!!sym(way)),1)) %>% filter(log > plast) ->data # 转换染色体编号 i <- 1 new <- data.frame(matrix(ncol = 2)) new <- new[-1,] for (x in c(1:7)){ for (y in c("A","B","D")){ chr <- paste0(x,y) # print(chr) new_add <- c(i,chr) new <- rbind(new,new_add) i <- i + 1 } } colnames(new) <- c("CHROM","chr") # 替换染色体编号 data %>% left_join(new,by = "CHROM") ->data2 data2$loc <- phe # 待标注的log值筛选 data2$logwt <- ifelse(data2$log > 10,paste0('log=',data2$log,sep=""),NA) data2$MB <- data2$POS/1000000 # 写出为中间结果 write_csv(data2,paste0("./out/1_GWAS_Result_txt2csv/",phe,".", way,".csv")) } }
获取文件列表:
使用list.files()函数，查找目录中所有满足文件名模式为"*.txt"的文件，并将它们的文件名存储在id_list列表中
遍历文件列表:
使用一个循环来遍历id_list中的每个文件名（变量名为id）
文件名解析:
使用str_split()函数，将文件名（id）按照"."分割成多个部分，存储在atom中
条件判断 (Group A):
检查atom中的部分数量是否为5，如果是，进入下一步
变量赋值:
从atom中获取不同部分的值，分别赋给phe、way和plast变量。 对way进行处理，将其中的"Farm"替换为空字符串
特殊值处理:
将plast的值与字符串"1e-5"进行比较，如果相等，将plast设置为6，否则将其转换为数值类型（这里如果需要人工规定阈值则修改此参数）
数据读取与处理:
使用read_delim()函数读取文件内容，并指定列类型。数据将被读入名为df的数据框。 修改df的第9列名称为way。 使用mutate()函数，计算并添加新的"log"列，值为way列的负对数的舍入值。 使用filter()函数，筛选出"log"列大于plast的行，结果存储在名为data的数据框中。
染色体编号转换:
使用两个嵌套循环，生成一个新的数据框new，其中存储了染色体编号的对应关系，生成的new数据框将被用于将染色体编号从数字格式转换为字符格式。
数据处理与替换:
使用left_join()函数，将之前筛选得到的data数据框与new数据框按照"CHROM"列进行连接，结果存储在data2中。 将data2中的"CHROM"列替换为之前提取的phe值。 根据条件，计算并添加"logwt"列和"MB"列。
结果输出:
使用write_csv()函数，将经过处理的data2数据框写入指定路径（"./out/1_GWAS_Result_txt2csv/"）下，以"phe.way.csv"的格式命名。
