GWAS analysis for salt-induced changes in root:shoot ratio in Arabidopsis
This notebook describes the selection of the top candidate genes using GWAS output for the data collected over 360 Arabidopsis accessions in HapMap population. The GWAS was run using ASReml script, as described in Awlia et al., 2021, using initially 250 K SNPs and subsequently 4 M SNPs. The raw GWAS output data can be found on Zenodo repository for mapping using 250 k SNP and 4 M SNP.
Loading and analyzing the 250 k SNP data:
my_GWAS_files <- list.files(pattern="rda")
list.files()## [1] "20221102_GWAS_notebook.nb.html"
## [2] "20221102_GWAS_notebook.Rmd"
## [3] "250kManhattan. R_125p0 .pdf"
## [4] "250kManhattan. R_75p0 .pdf"
## [5] "250kManhattan. reGF_0mM .pdf"
## [6] "250kManhattan. reGF_125mM .pdf"
## [7] "250kManhattan. reGF_75mM .pdf"
## [8] "250kManhattan. reSTART_0mM .pdf"
## [9] "250kManhattan. reSTART_125mM .pdf"
## [10] "250kManhattan. reSTART_75mM .pdf"
## [11] "250kManhattan. Root_d2_0mM .pdf"
## [12] "250kManhattan. Root_d2_125mM .pdf"
## [13] "250kManhattan. Root_d2_75mM .pdf"
## [14] "250kManhattan. Root_d4_0mM .pdf"
## [15] "250kManhattan. Root_d4_125mM .pdf"
## [16] "250kManhattan. Root_d4_75mM .pdf"
## [17] "250kManhattan. Root_d6_0mM .pdf"
## [18] "250kManhattan. Root_d6_125mM .pdf"
## [19] "250kManhattan. Root_d6_75mM .pdf"
## [20] "250kManhattan. Root_d8_0mM .pdf"
## [21] "250kManhattan. Root_d8_125mM .pdf"
## [22] "250kManhattan. Root_d8_75mM .pdf"
## [23] "250kManhattan. RpS_125p0 .pdf"
## [24] "250kManhattan. RpS_75p0 .pdf"
## [25] "250kManhattan. RpS_d2_0 .pdf"
## [26] "250kManhattan. RpS_d2_125 .pdf"
## [27] "250kManhattan. RpS_d2_75 .pdf"
## [28] "250kManhattan. RpS_d4_0 .pdf"
## [29] "250kManhattan. RpS_d4_125 .pdf"
## [30] "250kManhattan. RpS_d4_75 .pdf"
## [31] "250kManhattan. RpS_d6_0 .pdf"
## [32] "250kManhattan. RpS_d6_125 .pdf"
## [33] "250kManhattan. RpS_d6_75 .pdf"
## [34] "250kManhattan. RpS_d8_0 .pdf"
## [35] "250kManhattan. RpS_d8_125 .pdf"
## [36] "250kManhattan. RpS_d8_75 .pdf"
## [37] "250kManhattan. RpS_GF_0 .pdf"
## [38] "250kManhattan. RpS_GF_125 .pdf"
## [39] "250kManhattan. RpS_GF_75 .pdf"
## [40] "250kManhattan. S_125p0 .pdf"
## [41] "250kManhattan. S_75p0 .pdf"
## [42] "250kManhattan. seGF_0mM .pdf"
## [43] "250kManhattan. seGF_125mM .pdf"
## [44] "250kManhattan. seGF_75mM .pdf"
## [45] "250kManhattan. seSTART_0mM .pdf"
## [46] "250kManhattan. seSTART_125mM .pdf"
## [47] "250kManhattan. seSTART_75mM .pdf"
## [48] "250kManhattan. Shoot_d2_0mM .pdf"
## [49] "250kManhattan. Shoot_d2_125mM .pdf"
## [50] "250kManhattan. Shoot_d2_75mM .pdf"
## [51] "250kManhattan. Shoot_d4_0mM .pdf"
## [52] "250kManhattan. Shoot_d4_125mM .pdf"
## [53] "250kManhattan. Shoot_d4_75mM .pdf"
## [54] "250kManhattan. Shoot_d6_0mM .pdf"
## [55] "250kManhattan. Shoot_d6_125mM .pdf"
## [56] "250kManhattan. Shoot_d6_75mM .pdf"
## [57] "250kManhattan. Shoot_d8_0mM .pdf"
## [58] "250kManhattan. Shoot_d8_125mM .pdf"
## [59] "250kManhattan. Shoot_d8_75mM .pdf"
## [60] "250kManhattan. SpT_d2_0 .pdf"
## [61] "250kManhattan. SpT_d2_125 .pdf"
## [62] "250kManhattan. SpT_d2_75 .pdf"
## [63] "250kManhattan. SpT_d4_0 .pdf"
## [64] "250kManhattan. SpT_d4_125 .pdf"
## [65] "250kManhattan. SpT_d4_75 .pdf"
## [66] "250kManhattan. SpT_d6_0 .pdf"
## [67] "250kManhattan. SpT_d6_125 .pdf"
## [68] "250kManhattan. SpT_d6_75 .pdf"
## [69] "250kManhattan. SpT_d8_0 .pdf"
## [70] "250kManhattan. SpT_d8_125 .pdf"
## [71] "250kManhattan. SpT_d8_75 .pdf"
## [72] "250kManhattan. SpT_SHIIT_d2_125 .pdf"
## [73] "250kManhattan. SpT_SHIIT_d2_75 .pdf"
## [74] "250kManhattan. SpT_SHIIT_d4_125 .pdf"
## [75] "250kManhattan. SpT_SHIIT_d4_75 .pdf"
## [76] "250kManhattan. SpT_SHIIT_d6_125 .pdf"
## [77] "250kManhattan. SpT_SHIIT_d6_75 .pdf"
## [78] "250kManhattan. SpT_SHIIT_d8_125 .pdf"
## [79] "250kManhattan. SpT_SHIIT_d8_75 .pdf"
## [80] "250kManhattan. Tot_d2_0 .pdf"
## [81] "250kManhattan. Tot_d2_125 .pdf"
## [82] "250kManhattan. Tot_d2_75 .pdf"
## [83] "250kManhattan. Tot_d4_0 .pdf"
## [84] "250kManhattan. Tot_d4_125 .pdf"
## [85] "250kManhattan. Tot_d4_75 .pdf"
## [86] "250kManhattan. Tot_d6_0 .pdf"
## [87] "250kManhattan. Tot_d6_125 .pdf"
## [88] "250kManhattan. Tot_d6_75 .pdf"
## [89] "250kManhattan. Tot_d8_0 .pdf"
## [90] "250kManhattan. Tot_d8_125 .pdf"
## [91] "250kManhattan. Tot_d8_75 .pdf"
## [92] "250kQQplot. R_125p0 .pdf"
## [93] "250kQQplot. R_75p0 .pdf"
## [94] "250kQQplot. reGF_0mM .pdf"
## [95] "250kQQplot. reGF_125mM .pdf"
## [96] "250kQQplot. reGF_75mM .pdf"
## [97] "250kQQplot. reSTART_0mM .pdf"
## [98] "250kQQplot. reSTART_125mM .pdf"
## [99] "250kQQplot. reSTART_75mM .pdf"
## [100] "250kQQplot. Root_d2_0mM .pdf"
## [101] "250kQQplot. Root_d2_125mM .pdf"
## [102] "250kQQplot. Root_d2_75mM .pdf"
## [103] "250kQQplot. Root_d4_0mM .pdf"
## [104] "250kQQplot. Root_d4_125mM .pdf"
## [105] "250kQQplot. Root_d4_75mM .pdf"
## [106] "250kQQplot. Root_d6_0mM .pdf"
## [107] "250kQQplot. Root_d6_125mM .pdf"
## [108] "250kQQplot. Root_d6_75mM .pdf"
## [109] "250kQQplot. Root_d8_0mM .pdf"
## [110] "250kQQplot. Root_d8_125mM .pdf"
## [111] "250kQQplot. Root_d8_75mM .pdf"
## [112] "250kQQplot. RpS_125p0 .pdf"
## [113] "250kQQplot. RpS_75p0 .pdf"
## [114] "250kQQplot. RpS_d2_0 .pdf"
## [115] "250kQQplot. RpS_d2_125 .pdf"
## [116] "250kQQplot. RpS_d2_75 .pdf"
## [117] "250kQQplot. RpS_d4_0 .pdf"
## [118] "250kQQplot. RpS_d4_125 .pdf"
## [119] "250kQQplot. RpS_d4_75 .pdf"
## [120] "250kQQplot. RpS_d6_0 .pdf"
## [121] "250kQQplot. RpS_d6_125 .pdf"
## [122] "250kQQplot. RpS_d6_75 .pdf"
## [123] "250kQQplot. RpS_d8_0 .pdf"
## [124] "250kQQplot. RpS_d8_125 .pdf"
## [125] "250kQQplot. RpS_d8_75 .pdf"
## [126] "250kQQplot. RpS_GF_0 .pdf"
## [127] "250kQQplot. RpS_GF_125 .pdf"
## [128] "250kQQplot. RpS_GF_75 .pdf"
## [129] "250kQQplot. S_125p0 .pdf"
## [130] "250kQQplot. S_75p0 .pdf"
## [131] "250kQQplot. seGF_0mM .pdf"
## [132] "250kQQplot. seGF_125mM .pdf"
## [133] "250kQQplot. seGF_75mM .pdf"
## [134] "250kQQplot. seSTART_0mM .pdf"
## [135] "250kQQplot. seSTART_125mM .pdf"
## [136] "250kQQplot. seSTART_75mM .pdf"
## [137] "250kQQplot. Shoot_d2_0mM .pdf"
## [138] "250kQQplot. Shoot_d2_125mM .pdf"
## [139] "250kQQplot. Shoot_d2_75mM .pdf"
## [140] "250kQQplot. Shoot_d4_0mM .pdf"
## [141] "250kQQplot. Shoot_d4_125mM .pdf"
## [142] "250kQQplot. Shoot_d4_75mM .pdf"
## [143] "250kQQplot. Shoot_d6_0mM .pdf"
## [144] "250kQQplot. Shoot_d6_125mM .pdf"
## [145] "250kQQplot. Shoot_d6_75mM .pdf"
## [146] "250kQQplot. Shoot_d8_0mM .pdf"
## [147] "250kQQplot. Shoot_d8_125mM .pdf"
## [148] "250kQQplot. Shoot_d8_75mM .pdf"
## [149] "250kQQplot. SpT_d2_0 .pdf"
## [150] "250kQQplot. SpT_d2_125 .pdf"
## [151] "250kQQplot. SpT_d2_75 .pdf"
## [152] "250kQQplot. SpT_d4_0 .pdf"
## [153] "250kQQplot. SpT_d4_125 .pdf"
## [154] "250kQQplot. SpT_d4_75 .pdf"
## [155] "250kQQplot. SpT_d6_0 .pdf"
## [156] "250kQQplot. SpT_d6_125 .pdf"
## [157] "250kQQplot. SpT_d6_75 .pdf"
## [158] "250kQQplot. SpT_d8_0 .pdf"
## [159] "250kQQplot. SpT_d8_125 .pdf"
## [160] "250kQQplot. SpT_d8_75 .pdf"
## [161] "250kQQplot. SpT_SHIIT_d2_125 .pdf"
## [162] "250kQQplot. SpT_SHIIT_d2_75 .pdf"
## [163] "250kQQplot. SpT_SHIIT_d4_125 .pdf"
## [164] "250kQQplot. SpT_SHIIT_d4_75 .pdf"
## [165] "250kQQplot. SpT_SHIIT_d6_125 .pdf"
## [166] "250kQQplot. SpT_SHIIT_d6_75 .pdf"
## [167] "250kQQplot. SpT_SHIIT_d8_125 .pdf"
## [168] "250kQQplot. SpT_SHIIT_d8_75 .pdf"
## [169] "250kQQplot. Tot_d2_0 .pdf"
## [170] "250kQQplot. Tot_d2_125 .pdf"
## [171] "250kQQplot. Tot_d2_75 .pdf"
## [172] "250kQQplot. Tot_d4_0 .pdf"
## [173] "250kQQplot. Tot_d4_125 .pdf"
## [174] "250kQQplot. Tot_d4_75 .pdf"
## [175] "250kQQplot. Tot_d6_0 .pdf"
## [176] "250kQQplot. Tot_d6_125 .pdf"
## [177] "250kQQplot. Tot_d6_75 .pdf"
## [178] "250kQQplot. Tot_d8_0 .pdf"
## [179] "250kQQplot. Tot_d8_125 .pdf"
## [180] "250kQQplot. Tot_d8_75 .pdf"
## [181] "250kSNP.Bonf.csv"
## [182] "R_125p0_gwas.rda"
## [183] "R_75p0_gwas.rda"
## [184] "reGF_0mM_gwas.rda"
## [185] "reGF_125mM_gwas.rda"
## [186] "reGF_75mM_gwas.rda"
## [187] "reSTART_0mM_gwas.rda"
## [188] "reSTART_125mM_gwas.rda"
## [189] "reSTART_75mM_gwas.rda"
## [190] "Root_d2_0mM_gwas.rda"
## [191] "Root_d2_125mM_gwas.rda"
## [192] "Root_d2_75mM_gwas.rda"
## [193] "Root_d4_0mM_gwas.rda"
## [194] "Root_d4_125mM_gwas.rda"
## [195] "Root_d4_75mM_gwas.rda"
## [196] "Root_d6_0mM_gwas.rda"
## [197] "Root_d6_125mM_gwas.rda"
## [198] "Root_d6_75mM_gwas.rda"
## [199] "Root_d8_0mM_gwas.rda"
## [200] "Root_d8_125mM_gwas.rda"
## [201] "Root_d8_75mM_gwas.rda"
## [202] "RpS_125p0_gwas.rda"
## [203] "RpS_75p0_gwas.rda"
## [204] "RpS_d2_0_gwas.rda"
## [205] "RpS_d2_125_gwas.rda"
## [206] "RpS_d2_75_gwas.rda"
## [207] "RpS_d4_0_gwas.rda"
## [208] "RpS_d4_125_gwas.rda"
## [209] "RpS_d4_75_gwas.rda"
## [210] "RpS_d6_0_gwas.rda"
## [211] "RpS_d6_125_gwas.rda"
## [212] "RpS_d6_75_gwas.rda"
## [213] "RpS_d8_0_gwas.rda"
## [214] "RpS_d8_125_gwas.rda"
## [215] "RpS_d8_75_gwas.rda"
## [216] "RpS_GF_0_gwas.rda"
## [217] "RpS_GF_125_gwas.rda"
## [218] "RpS_GF_75_gwas.rda"
## [219] "S_125p0_gwas.rda"
## [220] "S_75p0_gwas.rda"
## [221] "seGF_0mM_gwas.rda"
## [222] "seGF_125mM_gwas.rda"
## [223] "seGF_75mM_gwas.rda"
## [224] "seSTART_0mM_gwas.rda"
## [225] "seSTART_125mM_gwas.rda"
## [226] "seSTART_75mM_gwas.rda"
## [227] "Shoot_d2_0mM_gwas.rda"
## [228] "Shoot_d2_125mM_gwas.rda"
## [229] "Shoot_d2_75mM_gwas.rda"
## [230] "Shoot_d4_0mM_gwas.rda"
## [231] "Shoot_d4_125mM_gwas.rda"
## [232] "Shoot_d4_75mM_gwas.rda"
## [233] "Shoot_d6_0mM_gwas.rda"
## [234] "Shoot_d6_125mM_gwas.rda"
## [235] "Shoot_d6_75mM_gwas.rda"
## [236] "Shoot_d8_0mM_gwas.rda"
## [237] "Shoot_d8_125mM_gwas.rda"
## [238] "Shoot_d8_75mM_gwas.rda"
## [239] "SpT_d2_0_gwas.rda"
## [240] "SpT_d2_125_gwas.rda"
## [241] "SpT_d2_75_gwas.rda"
## [242] "SpT_d4_0_gwas.rda"
## [243] "SpT_d4_125_gwas.rda"
## [244] "SpT_d4_75_gwas.rda"
## [245] "SpT_d6_0_gwas.rda"
## [246] "SpT_d6_125_gwas.rda"
## [247] "SpT_d6_75_gwas.rda"
## [248] "SpT_d8_0_gwas.rda"
## [249] "SpT_d8_125_gwas.rda"
## [250] "SpT_d8_75_gwas.rda"
## [251] "SpT_SHIIT_d2_125_gwas.rda"
## [252] "SpT_SHIIT_d2_75_gwas.rda"
## [253] "SpT_SHIIT_d4_125_gwas.rda"
## [254] "SpT_SHIIT_d4_75_gwas.rda"
## [255] "SpT_SHIIT_d6_125_gwas.rda"
## [256] "SpT_SHIIT_d6_75_gwas.rda"
## [257] "SpT_SHIIT_d8_125_gwas.rda"
## [258] "SpT_SHIIT_d8_75_gwas.rda"
## [259] "Tot_d2_0_gwas.rda"
## [260] "Tot_d2_125_gwas.rda"
## [261] "Tot_d2_75_gwas.rda"
## [262] "Tot_d4_0_gwas.rda"
## [263] "Tot_d4_125_gwas.rda"
## [264] "Tot_d4_75_gwas.rda"
## [265] "Tot_d6_0_gwas.rda"
## [266] "Tot_d6_125_gwas.rda"
## [267] "Tot_d6_75_gwas.rda"
## [268] "Tot_d8_0_gwas.rda"
## [269] "Tot_d8_125_gwas.rda"
## [270] "Tot_d8_75_gwas.rda"load(my_GWAS_files[1])
head(results)library(qqman)## ## For example usage please run: vignette('qqman')## ## Citation appreciated but not required:## Turner, (2018). qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. Journal of Open Source Software, 3(25), 731, https://doi.org/10.21105/joss.00731.## library(qvalue)
library(cowplot)results$Pval <-as.numeric(as.character(results$Pval))
results <- subset(results, results$MAC > 3)
a <- nrow(results)
b <- 0.05/a
head(results)-log10(b)## [1] 6.630361title <- gsub("_gwas.rda", "", my_GWAS_files[1])
manhattan(results, chr="Chr",bp="Pos", p="Pval", snp="SNP", main=title, ylim=c(0,12), col = c("aquamarine3", "bisque3"), suggestiveline=FALSE, genomewideline=-log10(b))
pdf(paste("250kManhattan.", title, ".pdf"))
manhattan(results, chr="Chr",bp="Pos", p="Pval", snp="SNP", main=title, ylim=c(0,12), col = c("aquamarine3", "bisque3"), suggestiveline=FALSE, genomewideline=-log10(b))
dev.off()## quartz_off_screen
## 2qq(results$Pval, main = title, cex = 1.5, las = 1)
pdf(paste("250kQQplot.", title, ".pdf"))
qq(results$Pval, main = title, cex = 1.5, las = 1)
dev.off()## quartz_off_screen
## 2results$Pval <-as.numeric(as.character(results$Pval))
results$LOD <- -log10(results$Pval)
results$trait <- gsub("_gwas.rda", "", my_GWAS_files[1])
LOD5 <- subset(results, results$LOD > 5)
LOD5length(my_GWAS_files)## [1] 89for(i in 2:89){
load(my_GWAS_files[i])
# Prepare table for graphing:
args <- commandArgs(trailingOnly=TRUE)
results$Pval <-as.numeric(as.character(results$Pval))
results <- subset(results, results$MAC > 3)
a <- nrow(results)
b <- 0.05/a
title <- gsub("_gwas.rda", "", my_GWAS_files[i])
pdf(paste("250kManhattan.", title, ".pdf"))
manhattan(results, chr="Chr",bp="Pos", p="Pval", snp="SNP", main=title, ylim=c(0,12), col = c("aquamarine3", "bisque3"), suggestiveline=FALSE, genomewideline=-log10(b))
dev.off()
pdf(paste("250kQQplot.", title, ".pdf"))
qq(results$Pval, main = title, cex = 1.5, las = 1)
dev.off()
results$Pval <-as.numeric(as.character(results$Pval))
results$LOD <- -log10(results$Pval)
results$trait <- gsub("_gwas.rda", "", my_GWAS_files[i])
temp <- subset(results, results$LOD > 5)
LOD5 <- rbind(LOD5, temp)
}dim(LOD5)## [1] 284 11unique(LOD5$trait)## [1] "R_125p0" "R_75p0" "reGF_0mM" "reGF_125mM"
## [5] "reGF_75mM" "reSTART_0mM" "reSTART_125mM" "reSTART_75mM"
## [9] "Root_d2_125mM" "Root_d2_75mM" "Root_d4_125mM" "Root_d4_75mM"
## [13] "Root_d6_0mM" "Root_d6_125mM" "Root_d8_0mM" "Root_d8_125mM"
## [17] "Root_d8_75mM" "RpS_125p0" "RpS_75p0" "RpS_d2_0"
## [21] "RpS_d2_75" "RpS_d4_0" "RpS_d4_125" "RpS_d4_75"
## [25] "RpS_d6_0" "RpS_d6_125" "RpS_d6_75" "RpS_d8_0"
## [29] "RpS_d8_125" "RpS_d8_75" "RpS_GF_0" "RpS_GF_125"
## [33] "RpS_GF_75" "S_125p0" "S_75p0" "seGF_125mM"
## [37] "seSTART_125mM" "seSTART_75mM" "Shoot_d2_125mM" "Shoot_d2_75mM"
## [41] "Shoot_d4_0mM" "Shoot_d4_125mM" "Shoot_d4_75mM" "Shoot_d6_0mM"
## [45] "Shoot_d6_125mM" "Shoot_d6_75mM" "Shoot_d8_0mM" "Shoot_d8_125mM"
## [49] "Shoot_d8_75mM" "SpT_d2_0" "SpT_d2_75" "SpT_d4_0"
## [53] "SpT_d4_125" "SpT_d4_75" "SpT_d6_0" "SpT_d6_125"
## [57] "SpT_d8_125" "SpT_d8_75" "SpT_SHIIT_d2_125" "SpT_SHIIT_d2_75"
## [61] "SpT_SHIIT_d4_125" "SpT_SHIIT_d4_75" "SpT_SHIIT_d6_125" "SpT_SHIIT_d6_75"
## [65] "SpT_SHIIT_d8_125" "SpT_SHIIT_d8_75" "Tot_d2_0" "Tot_d2_125"
## [69] "Tot_d2_75" "Tot_d4_0" "Tot_d4_125" "Tot_d4_75"
## [73] "Tot_d6_0" "Tot_d6_125" "Tot_d6_75" "Tot_d8_0"
## [77] "Tot_d8_125" "Tot_d8_75"LOD5write.csv(LOD5, "250kSNP.LOD5.csv", row.names=FALSE)-log10(b)## [1] 6.630361LOD_Bonf <- subset(LOD5, LOD5$LOD > -log10(b))
dim(LOD_Bonf)## [1] 6 11LOD_Bonfwrite.csv(LOD_Bonf, "250kSNP.Bonf.csv", row.names=FALSE)