TidyGWAS
GWAS分析關(guān)鍵結(jié)果之一是顯著性SNP位點(diǎn)的P值粉楚，通常多年份多地點(diǎn)多模型的GWAS分析將會(huì)產(chǎn)生很多結(jié)果文件怀喉，如何對(duì)這些數(shù)據(jù)進(jìn)行整理？
匯總這些結(jié)果啼县，并將顯著性的位點(diǎn)或區(qū)域找出來率翅，更加清晰的展示關(guān)鍵信息练俐。
今天介紹TidyGWAS結(jié)果整理新方法，前段時(shí)間曾發(fā)過一篇筆記（GWAS結(jié)果整理算法）冕臭，但是有一些地方比較繁瑣腺晾，仍有優(yōu)化空間。
本篇筆記將分享一種基于R語言自動(dòng)實(shí)現(xiàn)GWAS結(jié)果整理的升級(jí)版方法浴韭，通過優(yōu)化關(guān)鍵步驟和算法丘喻，將代碼量從2000行縮減到了400行，速度提高10倍以上念颈。
關(guān)鍵詞：GWAS泉粉、R語言、Tidyverse
前期準(zhǔn)備工作
軟件安裝
本次使用的R語言版本是R4.3.0榴芳，需要以下R包嗡靡，如果沒有安裝需要提前安裝。
# Install packages
install.packages("tidyverse")
install.packages("data.table")
install.packages("foreach")
install.packages("doParallel")
install.packages("stringr")
# Load libraries
library(tidyverse)
library(data.table)
library(foreach)
library(doParallel)
library(stringr)
項(xiàng)目文件
通常建議每個(gè)任務(wù)建立項(xiàng)目文件夾窟感，請(qǐng)新建一個(gè)文件夾并設(shè)置為工作目錄讨彼，然后創(chuàng)建如下文件結(jié)構(gòu)。
其中Ref是參考基因組信息柿祈，data子文件夾存放原始的數(shù)據(jù)哈误，可以包含很多個(gè)以txt結(jié)尾的文件，該文件內(nèi)容如下所示：
主要信息是SNP躏嚎、染色體蜜自、物理位置、顯著性P值卢佣。每個(gè)文件的命名方式是“類型.年份環(huán)境.模型.P閾值.txt”
核心操作步驟
人工設(shè)置參數(shù)
# 參數(shù)設(shè)置-----
prefix <- "xxx"
windows_near <- 300*1000
#默認(rèn)300kb內(nèi)為連續(xù)顯著區(qū)段
id_list <- list.files("./data/",pattern = "*.txt") # 待整理的所有文件
Ref <- fread("./Ref.csv") # 參考基因組信息
其實(shí)重荠，使用起來非常簡單，只用設(shè)置一個(gè)項(xiàng)目輸出結(jié)果前綴和窗口大小即可虚茶，這個(gè)窗口大小指的是在判定連續(xù)型顯著區(qū)域時(shí)最大的閾值戈鲁。
比如默認(rèn)300Kb仇参，假如某幾個(gè)顯著的SNP位點(diǎn)之間的物理距離都在300Kb之內(nèi)，則把它們當(dāng)做一個(gè)連續(xù)的顯著區(qū)段婆殿。
之后的步驟都是全自動(dòng)執(zhí)行诈乒，不用再進(jìn)行任何修改，如果您想使用此方法只需要直接運(yùn)行一遍全部代碼即可鸣皂。
自動(dòng)生成參數(shù)
id_df <- str_split(id_list,"[.]") %>% as.data.frame() %>% t() %>% as.data.frame()
type_list <- id_df$V1[!duplicated(id_df$V1)]
phe_list <- id_df$V2[!duplicated(id_df$V2)]
model_list <- id_df$V3[!duplicated(id_df$V3)]
p_value <- 1e-5 #顯著性
suffix <- ".txt"
write.csv(id_df,str_c(prefix,"_parment.csv"),row.names = F)
這段代碼的目的是將一個(gè)包含點(diǎn)號(hào)分隔字符串的列表（文件名稱列表）分割成多個(gè)部分抓谴，轉(zhuǎn)換為數(shù)據(jù)框暮蹂，然后從每一列中提取出不重復(fù)的元素寞缝，分別存儲(chǔ)在三個(gè)不同的列表中，這樣就得到了所有待整理的信息清單仰泻。
這份文件將會(huì)自動(dòng)輸出保存荆陆，記錄了所有待整理的文件信息和參數(shù)。
創(chuàng)建輸出文件夾
out_dirs <- c("1_GWAS_Result_txt2csv",
"2_SNP_Infomation",
"3_Gene_Maping_Result",
"4_Rebind_All_Output")
for (mydir in out_dirs){
if (dir.exists(paste0("./out/",mydir))){
cat(paste0("[check] ./out/",mydir," is exist !\n"))
}else{
dir.create(paste0("./out/",mydir))
cat(paste0("[check] ./out/",mydir," create finished !\n"))
}
}
這一步自動(dòng)檢測是否存在目標(biāo)文件夾集侯，如果不存在的話創(chuàng)建一個(gè)被啼，后續(xù)的中間結(jié)果和文件將自動(dòng)寫入這些文件夾。
Step1：文件整理
for (id in id_list){
file_name <- paste0("./data/",id)
atom <- str_split(id,"[.]")
type <- atom[[1]][1]
phe <- atom[[1]][2]
way <- atom[[1]][3] %>% str_replace("Farm","") # 將模型替換為CPU
plast <- atom[[1]][4]
# 特異性標(biāo)注P值并將其裝換為數(shù)字型
if (plast == "1e-5"){plast <- 6}else{plast <- as.numeric(plast)}
print(file_name)
# 計(jì)算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
# 轉(zhuǎn)換染色體編號(hào)
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")
# 替換染色體編號(hào)
data %>%
left_join(new,by = "CHROM") ->data2
data2$loc <- phe
# 待標(biāo)注的log值篩選
data2$logwt <- ifelse(data2$log > 10,paste0('log=',data2$log,sep=""),NA)
data2$MB <- data2$POS/1000000
# 寫出為中間結(jié)果
write_csv(data2,paste0("./out/1_GWAS_Result_txt2csv/",type,".",phe,".",
way,".csv"))
}
這段代碼主要用于處理基因組數(shù)據(jù)棠枉，涉及文件讀取浓体、數(shù)據(jù)分割、條件篩選辈讶、數(shù)據(jù)轉(zhuǎn)換和導(dǎo)出等步驟命浴。
Step2：統(tǒng)計(jì)顯著位點(diǎn)
Ref_chr <- Ref
all_single <- list() #匯總單標(biāo)記
all_near <- list() #匯總連續(xù)標(biāo)記
id_list_step2 <- list.files("./out/1_GWAS_Result_txt2csv/")
for (id in id_list_step2){
# 創(chuàng)建染色體
chr_list <- list()
for (tmp_chr in new$chr){
chr_list[[tmp_chr]] <- filter(Ref_chr,chr == tmp_chr)
}
# 開始計(jì)算----
file_name <- paste0("./out/1_GWAS_Result_txt2csv/",id)
atom <- str_split(id,"[.]")
# print(file_name)
data <- read_csv(file_name,show_col_types = FALSE)
loc <- data$loc[1] #
job <- paste0(atom[[1]][1],"_",atom[[1]][2],"_",atom[[1]][3])
way <- atom[[1]][3]
### 單標(biāo)記篩選 ========================================================================
# 計(jì)算基因位置間距
data$longH <- NA
data$longQ <- NA
data$class <- NA
# 顯著位點(diǎn)小于3個(gè)的情況下跳過
if (nrow(data) < 3){
next
}
for (i in 2:nrow(data)){
a <- data$POS[i]
i <- i+1
b <- data$POS[i]
i <- i-2
c <- data$POS[i]
i <- i+1
if(i == nrow(data)){
data$class[i] <- "wei"
break
}
if(a-c<0){
data$class[i] <- "shou"
next
}
if(b-a<0){
data$class[i] <- "wei"
next
}
data$longH[i] <- (b-a)
data$longQ[i] <- (a-c)
}
data$class[1] <- "shou"
# 對(duì)距離進(jìn)行區(qū)分,按照windows_near為區(qū)分閾值
for (i in 1:nrow(data)){
if (is.na(data$longH[i]) | is.na(data$longQ[i])){
next
}
if (data$longH[i]>windows_near & data$longQ[i]>windows_near){
data$class[i] <- "single"
}
}
# 單標(biāo)記位點(diǎn)處理
data$ws <- ifelse(is.na(data$logwt),
paste0(data$SNP,",Find in ",str_replace(id,".csv",""),sep=""),
paste0(data$SNP,",Find in ",str_replace(id,".csv",""),",",data$logwt,sep=""))
### 單標(biāo)記信息位置注釋 ===================================================================
# 單標(biāo)記位置信息寫入single
single <- data.frame(matrix(ncol = 4))
single <- single[-1,]
colnames(single) <- c("positon","info","chr","loc")
for (i in 1:nrow(data)){
tem_class <- data$class[i]
tem_add <- c(data$POS[i],data$ws[i],data$chr[i],data$loc[i])
ifelse(tem_class == 'single',single <- rbind(single,tem_add),"1")
ifelse(tem_class == 'shou',single <- rbind(single,tem_add),"2")
ifelse(tem_class == 'wei',single <- rbind(single,tem_add),"3")
}
colnames(single) <- c("positon","info","chr","loc")
### 連續(xù)區(qū)間篩選 ===================================================================
near <- data.frame(matrix(ncol = 5)) #初始化矩陣
near <- near[-1,]
colnames(near) <- c("p1","p2","info","chr","number")
for (x in c(1:7)){
for (y in c("A","B","D")){
chr_id <- paste0(x,y,sep="")
foot <- 0 # 步長，用于迭代計(jì)算閱讀框
for (i in which(data$chr==chr_id)){
if (sum(data$chr==chr_id)<2){next} #若某個(gè)染色體的位點(diǎn)數(shù)小于3則跳過
if (foot>0){ # 如果foot變量大于0贱除，說明兩個(gè)位點(diǎn)存在跨越關(guān)系生闲，進(jìn)行歸零
foot <- foot - 1 #
next
}
n_pos_1 <- data$POS[i]
for (m in which(data$chr==chr_id)){
n_pos_2 <- data$POS[m]
if (n_pos_2 - n_pos_1 < 0){ # 后一個(gè)值小于前一個(gè)值時(shí)跳過
next
}
else{
if (n_pos_2 - n_pos_1 == 0){ # 兩個(gè)值相等時(shí)跳過
next
}
else{
if (n_pos_2 - n_pos_1 < windows_near){ # 任意兩個(gè)位點(diǎn)距離小于預(yù)設(shè)窗口大小
foot <- foot+1 # 向前進(jìn)行一步，跨越一個(gè)位點(diǎn)
if (is.na(data$class[i])){
data$class[i] <- "near" # 如果此時(shí)位點(diǎn)尚不屬于single月幌、shou碍讯、wei，則為near
}
}
else{
if (foot > 10){
n_wn <- paste0(data$SNP[i],"-",data$SNP[m-1],", [",foot,"],Find in ",job)
}else{
n_wn <- paste0(data$SNP[i],"-",data$SNP[m-1],",Find in ",job)
}
n_add <- c(n_pos_1,data$POS[m-1],n_wn,data$chr[i],foot)
if (is.na(data$class[i])){
data$class[i] <- "near"
}
break
}
}
}
}
if (length(data$POS[m-1])>0){
if (n_pos_1 !=data$POS[m-1]){
near <- rbind(near,n_add)
n_add <- c()
}
}
}
}
}
colnames(near) <- c("p1","p2","info","chr","number")
# 刪除重復(fù)行
near_new <- data.frame(matrix(ncol =5))
near_new <- near_new[-1,]
for (i in 1:nrow(near)){
if (!identical(near$p1[i],near$p2[i])){
new_add <- c(near$p1[i],near$p2[i],near$info[i],near$chr[i],near$number[i])
near_new <- rbind(near_new,new_add)
}
}
colnames(near_new) <- c("p1","p2","info","chr","number")
near <- near_new
### 連續(xù)標(biāo)記回帖參考基因組----
OK <- 0 #成功添加的個(gè)數(shù)
for (chr in names(chr_list)){
for (i in 1:nrow(near)){
if (identical(near$chr[i],chr)){
my_a <- which.min(abs(as.numeric(near$p1[i]) - as.numeric(chr_list[[chr]]$X3G.Start.1)))
my_b <- which.min(abs(as.numeric(near$p2[i]) - as.numeric(chr_list[[chr]]$X3G.Start.1)))
chr_list[[chr]]$out[my_a:my_b] <- near$info[i]
OK <- OK + (my_b - my_a)
}
}
cli::cli_alert_success(str_c("[Chromosomes are currently being processed]：",chr))
}
num_near <- OK
### 迭代添加單標(biāo)記信息----
OK <- 0 #成功添加的個(gè)數(shù)
for (chr in names(chr_list)){
for (i in 1:nrow(single)){
if (identical(single$chr[i],chr)){
index <- which.min(abs(as.numeric(single$positon[i]) - as.numeric(chr_list[[chr]]$X3G.Start.1)))
chr_list[[chr]]$out[index] <- single$info[i]
OK <- OK+1
}
}
cli::cli_alert_success(str_c("[Chromosomes are currently being processed]：",chr))
}
num_single <- OK
### 迭代添加首尾標(biāo)記 =====================================================================
cat(paste0("[",str_sub(data$ws[1],nchar(data$ws[1])-2,nchar(data$ws[1])),
"-",data$loc[1],"] \t total near mark: ",num_near," \t total single mark: ",num_single,"\n"))
all_near[[id]] <- near
all_single[[id]] <- single
write_excel_csv(single,paste0("./out/2_SNP_Infomation/",id,"_single.csv"))
write_excel_csv(near,paste0("./out/2_SNP_Infomation/",id,"_near.csv"))
write_excel_csv(data,paste0("./out/2_SNP_Infomation/",id,"_data.csv"),na = "near")
out <- do.call(rbind,lapply(chr_list,function(x)x))
write_excel_csv(out,paste0("./out/3_Gene_Maping_Result/",job,".snp.csv"),na = "")
}
這段R語言代碼執(zhí)行了一系列復(fù)雜的數(shù)據(jù)處理操作扯躺，主要用于處理基因組關(guān)聯(lián)研究（GWAS）中的SNP（單核苷酸多態(tài)性）數(shù)據(jù)捉兴，包括識(shí)別顯著的單個(gè)位點(diǎn)和連續(xù)區(qū)間，以及將這些信息映射到參考基因組上录语，這對(duì)于理解基因與特定表型之間的關(guān)系非常重要倍啥。
Step3：標(biāo)注
id_list_step3 <- list.files("./out/3_Gene_Maping_Result/")
tem <- Ref[1:nrow(out),1:7] # 需要格外注意：tem和out文件必須一一對(duì)應(yīng)
index <- 8 #從第8列開始標(biāo)注
# 提取并標(biāo)注注釋信息
for (id in id_list_step3){
now <- read.csv(paste0("./out/3_Gene_Maping_Result/",id))
atom <- str_split(id,"[.]")
job <- atom[[1]][1]
tem <- bind_cols(tem,now[,8])
colnames(tem)[index] <- job
index <- index + 1
# print(id)
}
# 將多列信息合并為一列（優(yōu)化算法）
tem <- tem %>% as.data.frame()
tem$all <- NA
tem_rm_na <- tem[,colSums(is.na(tem)) < nrow(tem)]
tem_rm_na_info <- tem_rm_na[,9:ncol(tem_rm_na)-1]
tem_rm_na$all <- apply(tem_rm_na_info, 1, function(x){
x <- na.omit(x) # 刪除NA值
x <- x[nchar(x) > 3] # 保留字符長度大于3的元素
paste(x, collapse = " ; ") # 使用分號(hào)作為分隔符連接字符串
}
)
這段代碼的主要目的是將多個(gè)基因組映射結(jié)果文件中的注釋信息提取出來，并合并到一個(gè)主數(shù)據(jù)框中钦无。每個(gè)文件的特定列（通常是第8列）被提取并添加到tem數(shù)據(jù)框中逗栽，最后將這些信息合并，以便進(jìn)一步的分析和解釋失暂。
結(jié)果保存與輸出
final_out <- cbind(tem_rm_na[,1:7],tem_rm_na$all,tem_rm_na[,8:(ncol(tem_rm_na)-1)])
write_tsv(final_out,str_c("./out/4_Rebind_All_Output/",prefix,"_IT_DS_MLM_CPU.Output.final.tsv"))
這一步是為了將最終結(jié)果整理輸出保存彼宠，自動(dòng)根據(jù)項(xiàng)目名稱建立結(jié)果文件鳄虱。tsv文件可以直接選擇以Excel打開，就是常規(guī)表格格式凭峡。
查看結(jié)果
正常情況下拙已，運(yùn)行完上述流程后，能夠在out文件夾發(fā)現(xiàn)如上信息摧冀。
其中第一個(gè)文件夾儲(chǔ)存了原始文件轉(zhuǎn)換后的結(jié)果倍踪，第二個(gè)文件夾儲(chǔ)存了每個(gè)SNP的詳細(xì)信息，第三個(gè)文件夾是顯著區(qū)域回帖到參考基因組的結(jié)果索昂，第四個(gè)文件夾內(nèi)是最終的一個(gè)結(jié)果文件建车。
其中最后一個(gè)結(jié)果文件很重要，包含了所有的顯著信息椒惨，并對(duì)多環(huán)境同時(shí)共定位到的位點(diǎn)進(jìn)行標(biāo)注缤至，可以用于后續(xù)研究。
補(bǔ)充：優(yōu)化思路與方法
在寫代碼的時(shí)候康谆，最開始并沒有想到向量化編程的思路领斥，因此在早期版本中采用for循環(huán)迭代，速度巨慢沃暗。
for (i in 1:nrow(tem)){
var_add <- c()
for (m in 8:(ncol(tem)-1)){
if (tem[i,m] == ""){
next
}else{
var_add <- c(var_add,tem[i,m])
}
}
add_info <- str_c(var_add,sep = "",collapse = " ; ")
tem$all[i] <- add_info
}
該流程中最耗時(shí)的步驟是對(duì)結(jié)果進(jìn)行合并月洛，也就是Step3中將不同年份、地點(diǎn)孽锥、類型嚼黔、模型的顯著性關(guān)鍵信息進(jìn)行整合，合并為一列信息忱叭。
在實(shí)際計(jì)算中隔崎，這個(gè)數(shù)據(jù)維度大概是幾十萬行，每行進(jìn)行依次迭代的速度很慢韵丑，由于計(jì)算過程中不需要考慮不同行之間的相互影響爵卒，因此考慮改成多線程并行計(jì)算，同時(shí)在CPU上計(jì)算多行數(shù)據(jù)撵彻。
num_cores <- parallel::detectCores() # 設(shè)置線程數(shù)
cl <- makeCluster(num_cores)
registerDoParallel(cl)
# 原始數(shù)據(jù)框?yàn)閠em
n <- nrow(tem)
result <- foreach(i = 1:n, .combine = rbind) %dopar% {
var_add <- c()
for (m in 8:(ncol(tem)-1)){
# 如果出現(xiàn)缺失則跳過
if (is.na(tem[i,m])){next}
# 如果出現(xiàn)空位則跳過
if (tem[i,m] == ""){
next
} else {
# 追加新結(jié)果
var_add <- c(var_add, tem[i,m])
}
}
add_info <- stringr::str_c(var_add, sep = "", collapse = " ; ")
tem$all[i] <- add_info
return(tem[i, ])
}
# 關(guān)閉并行計(jì)算
stopCluster(cl)
registerDoSEQ()
并行計(jì)算的速度能有一定提升钓株，理論上64核心處理器的速度會(huì)比單純for循環(huán)提高幾十倍，但是缺陷也比較明顯陌僵，這個(gè)在計(jì)算的過程中每個(gè)線程都會(huì)復(fù)制一份內(nèi)存空間轴合，導(dǎo)致內(nèi)存占用量攀升。
最佳的方法是采用R語言向量化編程碗短，使用apply函數(shù)直接按行應(yīng)用函數(shù)受葛，這個(gè)速度嘎嘎快，而且還節(jié)省內(nèi)存空間。
tem_rm_na$all <- apply(tem_rm_na_info, 1, function(x){
x <- na.omit(x) # 刪除NA值
x <- x[nchar(x) > 3] # 保留字符長度大于3的元素
paste(x, collapse = " ; ") # 使用分號(hào)作為分隔符連接字符串
}
這回看著比較優(yōu)雅总滩，運(yùn)行速度也相對(duì)提升了一大截纲堵。
另外還有一個(gè)地方需要進(jìn)行優(yōu)化，在不同染色體的分界處需要考慮首尾位置闰渔，每條染色體之間是獨(dú)立的席函，同一條染色體是按物理位置依次排序，因此確定邊界很重要冈涧。
以下是原來計(jì)算首尾SNP的方法茂附，邏輯是根據(jù)當(dāng)前SNP物理位置與上一行SNP物理位置的大小來比較，如果是結(jié)束位置督弓，那么當(dāng)前SNP減去下一行的值為正值营曼，否則為負(fù)值。
if(i == nrow(data)){
data$class[i] <- "wei"
break
}
if(a-c<0){
data$class[i] <- "shou"
next
}
if(b-a<0){
data$class[i] <- "wei"
next
}
上述算法有個(gè)隱藏BUG咽筋，當(dāng)SNP數(shù)量多的時(shí)候能夠正常判斷溶推，但是當(dāng)SNP數(shù)量只有幾個(gè)的時(shí)候，有可能會(huì)出現(xiàn)某條染色體上最后一個(gè)顯著的SNP恰好比下一條染色體的第一條SNP位置大奸攻，此時(shí)算法會(huì)將其認(rèn)為是同一條染色體。
為了解決上述問題虱痕，重新修改了判定SNP首尾位置的算法睹耐，采用染色體信息直接判斷：
# 更新判斷SNP首尾位置的方法
if (data$chr[i-1] != data$chr[i] &
data$chr[i+1] == data$chr[i]){
data$class[i] <- "shou"
next
}else{
if (data$chr[i-1] == data$chr[i] &
data$chr[i+1] != data$chr[i]){
data$class[i] <- "wei"
next
}else{
if (data$chr[i-1] != data$chr[i] &
data$chr[i+1] != data$chr[i]){
data$class[i] <- "single"
next
}
}
}
本文由mdnice多平臺(tái)發(fā)布
