--- title: "Summarise GWAS results" author: "Malachy Campbell" date: "8/7/2018" output: rmdformats::html_clean: fig_width: 6 fig_height: 6 highlight: kate thumbnails: true lightbox: true gallery: true --- ```{r setup, include=FALSE, echo = F, eval = T} knitr::opts_knit$set(root.dir = '~/Documents/Dropbox/Work/Manuscripts/2018_RandomRegressionGWAS/ThePlantGenome/Revision/New Analysis/') ``` # Background Here, we'll generate a quick summary of the GWAS results, make some manhattan plots, and compare marker effects between RR and TP approaches. # Sumarise RR and TP results Here, we determined how many SNPs showed a significant ($p < 1 \times 10^{-4}$) association at one or more time points. A significance threshold was selected based on the criteria used by [Zhao et al (2011)](https://www.nature.com/articles/ncomms1467). A very small number was added to all $p$-values to correct for those that had a $p$-value of 0. This was to ensure that $-log_{10}(p)$ would return a real number. ```{r Sumarise RR and TP results, echo=TRUE, eval=FALSE} TP.p <- as.matrix( read.table("DataAndCode/CompCluster/TP/pval.mat.txt", sep = "\t", header = T, row.names = 1) ) RR.p <- as.matrix( read.table("DataAndCode/CompCluster/RR/pval.mat.txt", sep = "\t", header = T, row.names = 1) ) MAP <- read.table("DataAndCode/CompCluster/MAP_maf0.05.txt", header = F, sep = "\t") row.names(TP.p) <- MAP$V2 row.names(RR.p) <- MAP$V2 #Take the significant SNPs from RR GWAS and TP GWAS and store the names in an array RR.array <- array() TP.array <- array() for (i in 1:20){ if(i == 1){ RR.array <- row.names(RR.p[RR.p[,i] < 1E-4 ,]) TP.array <- row.names(TP.p[TP.p[,i] < 1E-4 ,]) }else{ RR.array <- c(RR.array, row.names(RR.p[RR.p[,i] < 1E-4 ,]) ) TP.array <- c(TP.array, row.names(TP.p[TP.p[,i] < 1E-4 ,]) ) } } length(unique(RR.array)) #31 length(unique(TP.array)) #38 ``` # Compare SNP effects between RR and TP approaches The following code was used to examine the relationship between the SNP effects from the two approaches (RR and TP). ```{r compare RR and TP, echo=TRUE, eval=FALSE} library(reshape2) #load the files TP.p <- as.matrix( read.table("DataAndCode/CompCluster/TP/pval.mat.txt", sep = "\t", header = T, row.names = 1) ) RR.p <- as.matrix( read.table("DataAndCode/CompCluster/RR/pval.mat.txt", sep = "\t", header = T, row.names = 1) ) TP.beta <- as.matrix( read.csv("DataAndCode/CompCluster/TP/Eff.mat.csv") ) RR.beta <- as.matrix( read.csv("DataAndCode/CompCluster/RR/Eff.mat.csv") ) MAP <- read.table("DataAndCode/CompCluster/MAP_maf0.05.txt", header = F, sep = "\t") #add SNP names and days to p and beta matrices row.names(TP.p) <- MAP$V2 row.names(RR.p) <- MAP$V2 row.names(TP.beta) <- MAP$V2 row.names(RR.beta) <- MAP$V2 colnames(TP.p) <- paste0("D", 1:20) colnames(RR.p) <- paste0("D", 1:20) colnames(TP.beta) <- paste0("D", 1:20) colnames(RR.beta) <- paste0("D", 1:20) #Convert the p values and snp effects to long format TP.p <- melt(TP.p, id.vars = 0) TP.beta <- melt(TP.beta, id.vars = 0) RR.p <- melt(RR.p, id.vars = 0) RR.beta <- melt(RR.beta, id.vars = 0) #Make a column called Key to provide a unique identifier for the SNP and day RR.p$Key <- paste0(RR.p$Var1, "_", RR.p$Var2) RR.beta$Key <- paste0(RR.beta$Var1, "_", RR.beta$Var2) TP.p$Key <- paste0(TP.p$Var1, "_", TP.p$Var2) TP.beta$Key <- paste0(TP.beta$Var1, "_", TP.beta$Var2) #Check the correlation for all SNPs cor(abs(TP.beta$value), abs(RR.beta$value)) #0.8504224 #Select only the signifant SNPs Sig.TP <- unique(TP.p[TP.p$value < 1E-4 ,]$Var1) Sig.RR <- unique(RR.p[RR.p$value < 1E-4 ,]$Var1) #Get the min and max effects for signficant SNPs from both approaches min(RR.beta[RR.beta$Var1 %in% Sig.RR ,]$value) #-299.1024 max(RR.beta[RR.beta$Var1 %in% Sig.RR ,]$value) #246.0745 min(TP.beta[TP.beta$Var1 %in% Sig.TP ,]$value) #-99.03464 max(TP.beta[TP.beta$Var1 %in% Sig.TP ,]$value) #112.3483 min(TP.beta$value) #-104.59 max(TP.beta$value) #112.3483 min(RR.beta$value) #-299.1024 max(RR.beta$value) #295.0326 Sig.TP <- TP.p[TP.p$value < 1E-4 ,]$Key #442 Sig.RR <- RR.p[RR.p$value < 1E-4 ,]$Key #191 ``` ## Plot SNP effects for RR and TP approaches This code was used to create figure 1. ```{r code for figure 1, echo=TRUE, eval=FALSE} pdf("DataAndCode/Figures/Fig1.pdf", h=6, w=3.25, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,.75), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,1), c(2,2), c(3,3))) plot(TP.beta$value, RR.beta$value, pch = 19, cex = 0.3, xlab = expression(beta["TP"]), ylab = expression(beta["RR"]) ) points(TP.beta[TP.beta$Key %in% Sig.RR ,]$value, RR.beta[RR.beta$Key %in% Sig.RR ,]$value, pch = 19, cex = 0.3, col = "red") abline(c(0,1), col = "grey", lty = 2) mtext("A", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) ## hist(TP.beta$value, main = "", xlab = expression(beta["TP"])) mtext("B", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) ## hist(RR.beta$value, main = "", xlab = expression(beta["RR"]), xlim = c(-400, 400)) mtext("C", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) dev.off() ``` # Manhattan Plots This code was used to generate a manhattan plot of SNP effects ($|\beta|$) for the RR and TP approaches. The script uses a manhattan plot function that was modified from Stephen Turner's code he posted a few years back ([check it out here](http://www.gettinggeneticsdone.com/2011/04/annotated-manhattan-plots-and-qq-plots.html)). The modified function plots the $|\beta|$ rather than $-log_{10}(p)$. ```{r manhattan plot function, echo=TRUE, eval=FALSE} manhattan.Beta <- function(dataframe, colors = c("gray10", "gray50"), ymax = "max", xaxis.cex = 1, yaxis.cex = 1, limitchromosomes = 1:23, suggestiveline = NULL, genomewideline = NULL, annotate=NULL, Title, ...) { d=dataframe ymax=max(d$Beta) ymin=min(d$Beta) #throws error if you don't have columns named CHR, BP, and P in your data frame. if (!("CHR" %in% names(d) & "BP" %in% names(d) & "Beta" %in% names(d))) stop("Make sure your data frame contains columns CHR, BP, and Beta") # limits chromosomes to plot. (23=x, 24=y, 25=par?, 26=mito?) if (any(limitchromosomes)) d=d[d$CHR %in% limitchromosomes, ] # remove na's, sort by CHR and BP, and keep snps where 00 & P<=1)) # -log10(p-value) d$logp = -log10(d$P) # sets colors based on colors argument. colors <- rep(colors,max(d$CHR))[1:max(d$CHR)] # sets the maximum value on the y axis (on the -log10p scale). if (ymax=="max") ymax<-ceiling(max(d$logp)) if (ymax<8) ymax<-8 # creates continuous position markers for x axis for entire chromosome. also creates tick points. d$pos=NA ticks=NULL lastbase=0 numchroms=length(unique(d$CHR)) if (numchroms==1) { d$pos=d$BP ticks=floor(length(d$pos))/2+1 } else { for (i in unique(d$CHR)) { if (i==1) { d[d$CHR==i, ]$pos=d[d$CHR==i, ]$BP } else { lastbase=lastbase+tail(subset(d,CHR==i-1)$BP, 1) d[d$CHR==i, ]$pos=d[d$CHR==i, ]$BP+lastbase } ticks=c(ticks, d[d$CHR==i, ]$pos[floor(length(d[d$CHR==i, ]$pos)/2)+1]) } } # create the plot with(d, plot(pos, logp, ylim = c(0,ymax), ylab=expression(-log[10](italic(p))), xlab = "Chromosome", xaxt = "n", type = "n", cex = 0.3, yaxt = "n", main = Title, ...)) # then make an axis that has chromosome number instead of position axis(1, at = ticks, lab = unique(d$CHR), cex.axis = xaxis.cex) axis(2, cex.axis = yaxis.cex) icol=1 for (i in unique(d$CHR)) { with(d[d$CHR==i, ],points(pos, logp, col=colors[icol], cex=0.3, ...)) icol = icol+1 } # create a new data frame with rows from the original data frame where SNP is in annotate character vector. # then plot those points over the original graph, but with a larger point size and a different color. if (!is.null(annotate)) { d.annotate=d[which(d$SNP %in% annotate), ] with(d.annotate, points(pos, logp, col="red", cex=0.5, ...)) } # add threshold lines if (suggestiveline) abline(h=suggestiveline, col="blue") if (genomewideline) abline(h=genomewideline, col="red") } ``` ```{r manhattan plots, echo=TRUE, eval=FALSE} library(reshape2) TP.p <- as.matrix( read.table("DataAndCode/CompCluster/TP/pval.mat.txt", sep = "\t", header = T, row.names = 1) ) RR.p <- as.matrix( read.table("DataAndCode/CompCluster/RR/pval.mat.txt", sep = "\t", header = T, row.names = 1) ) MAP <- read.table("DataAndCode/CompCluster/MAP_maf0.05.txt", header = F, sep = "\t") TP.beta <- as.matrix( read.csv("DataAndCode/CompCluster/TP/Eff.mat.csv") ) RR.beta <- as.matrix( read.csv("DataAndCode/CompCluster/RR/Eff.mat.csv") ) row.names(TP.p) <- MAP$V2 row.names(RR.p) <- MAP$V2 row.names(TP.beta) <- MAP$V2 row.names(RR.beta) <- MAP$V2 Sig.RR.p <- melt(RR.p, id.vars = 0) Sig.RR.p$Var2 <- as.numeric(as.character(sub("V", "", Sig.RR.p$Var2))) Sig.RR.p <- Sig.RR.p[Sig.RR.p$value < 1E-04 ,] Sig.TP.p <- melt(TP.p, id.vars = 0) Sig.TP.p$Var2 <- as.numeric(as.character(sub("V", "", Sig.TP.p$Var2))) Sig.TP.p <- Sig.TP.p[Sig.TP.p$value < 1E-04 ,] for (i in 1:20){ if(i < 10){ pdf("DataAndCode/Figures/FigS1.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 1:10){ tmp.Sig <- Sig.RR.p[Sig.RR.p$Var2 == i ,] tmp.RR <- cbind(MAP, RR.beta[,i]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "Beta") manhattan.Beta(tmp.RR, Title = paste0("RR D", i), annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75) mtext(LETTERS[i], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() }else{ pdf("DataAndCode/Figures/FigS2.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 11:20){ tmp.Sig <- Sig.RR.p[Sig.RR.p$Var2 == i ,] tmp.RR <- cbind(MAP, RR.beta[,i]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "Beta") manhattan.Beta(tmp.RR, Title = paste0("RR D", i), annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75) mtext(LETTERS[i-10], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() } } for (i in 1:20){ if(i < 10){ pdf("DataAndCode/Figures/FigS3.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 1:10){ tmp.Sig <- Sig.TP.p[Sig.TP.p$Var2 == i ,] tmp.TP <- cbind(MAP, TP.beta[,i]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "Beta") manhattan.Beta(tmp.TP, Title = paste0("TP D", i), annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75) mtext(LETTERS[i], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() }else{ pdf("DataAndCode/Figures/FigS4.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 11:20){ tmp.Sig <- Sig.TP.p[Sig.TP.p$Var2 == i ,] tmp.TP <- cbind(MAP, TP.beta[,i]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "Beta") manhattan.Beta(tmp.TP, Title = paste0("TP D", i), annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75) mtext(LETTERS[i-10], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() } } pdf("DataAndCode/Figures/Fig2.pdf", h=4, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4))) tmp.Sig <- Sig.RR.p[Sig.RR.p$Var2 == 1 ,] tmp.RR <- cbind(MAP, RR.beta[,1]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "Beta") manhattan.Beta(tmp.RR, Title = "RR D1", annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75) mtext("A", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) tmp.Sig <- Sig.RR.p[Sig.RR.p$Var2 == 20 ,] tmp.RR <- cbind(MAP, RR.beta[,20]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "Beta") manhattan.Beta(tmp.RR, Title = "RR D20", annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75) mtext("B", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) tmp.Sig <- Sig.TP.p[Sig.TP.p$Var2 == 1 ,] tmp.TP <- cbind(MAP, TP.beta[,1]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "Beta") manhattan.Beta(tmp.TP, Title = "TP D1", annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75) mtext("C", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) tmp.Sig <- Sig.TP.p[Sig.TP.p$Var2 == 20 ,] tmp.TP <- cbind(MAP, TP.beta[,20]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "Beta") manhattan.Beta(tmp.TP, Title = "TP D20", annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75) mtext("D", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) dev.off() ##Manhattan plots for p-values for (i in 1:20){ if(i < 10){ pdf("DataAndCode/Figures/FigS1_alt.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 1:10){ tmp.Sig <- Sig.RR.p[Sig.RR.p$Var2 == i ,] tmp.RR <- cbind(MAP, RR.p[,i]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "P") manhattan(tmp.RR, Title = paste0("RR D", i), ymax = 8, annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75, suggestiveline = F, genomewideline = F) mtext(LETTERS[i], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() }else{ pdf("DataAndCode/Figures/FigS2_alt.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 11:20){ tmp.Sig <- Sig.RR.p[Sig.RR.p$Var2 == i ,] tmp.RR <- cbind(MAP, RR.p[,i]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "P") manhattan(tmp.RR, Title = paste0("RR D", i), ymax = 8, annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75, suggestiveline = F, genomewideline = F) mtext(LETTERS[i-10], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() } } for (i in 1:20){ if(i < 10){ pdf("DataAndCode/Figures/FigS3_alt.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 1:10){ tmp.Sig <- Sig.TP.p[Sig.TP.p$Var2 == i ,] tmp.TP <- cbind(MAP, TP.p[,i]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "P") manhattan(tmp.TP, Title = paste0("TP D", i), ymax = 8, annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75, suggestiveline = F, genomewideline = F) mtext(LETTERS[i], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() }else{ pdf("DataAndCode/Figures/FigS4_alt.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 11:20){ tmp.Sig <- Sig.TP.p[Sig.TP.p$Var2 == i ,] tmp.TP <- cbind(MAP, TP.p[,i]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "P") manhattan(tmp.TP, Title = paste0("TP D", i), ymax = 8, annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75, suggestiveline = F, genomewideline = F) mtext(LETTERS[i-10], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() } } pdf("DataAndCode/Figures/Fig2_alt.pdf", h=4, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4))) tmp.Sig <- Sig.RR.p[Sig.RR.p$Var2 == 1 ,] tmp.RR <- cbind(MAP, RR.p[,1]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "P") manhattan(tmp.RR, Title = paste0("RR D1"), ymax = 8, annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75, suggestiveline = F, genomewideline = F) mtext("A", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) tmp.Sig <- Sig.RR.p[Sig.RR.p$Var2 == 20 ,] tmp.RR <- cbind(MAP, RR.p[,20]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "P") manhattan(tmp.RR, Title = paste0("RR D20"), ymax = 8, annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75, suggestiveline = F, genomewideline = F) mtext("B", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) tmp.Sig <- Sig.TP.p[Sig.TP.p$Var2 == 1 ,] tmp.TP <- cbind(MAP, TP.p[,1]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "P") manhattan(tmp.TP, Title = paste0("TP D1"), ymax = 8, annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75, suggestiveline = F, genomewideline = F) mtext("C", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) tmp.Sig <- Sig.TP.p[Sig.TP.p$Var2 == 20 ,] tmp.TP <- cbind(MAP, TP.p[,20]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "P") manhattan(tmp.TP, Title = paste0("TP D20"), ymax = 8, annotate = tmp.Sig$Var1, xaxis.cex = 0.75, yaxis.cex = 0.75, suggestiveline = F, genomewideline = F) mtext("D", 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) dev.off() ``` ## QQ-plots ```{r, echo=TRUE, eval=FALSE} library(qqman) ##TP for (i in 1:20){ if(i < 10){ pdf("DataAndCode/Figures/FigS6.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 1:10){ tmp.TP <- cbind(MAP, TP.p[,i]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "P") qq(tmp.TP$P, main = paste0("D", i)) mtext(LETTERS[i], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() }else{ pdf("DataAndCode/Figures/FigS7.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 11:20){ tmp.TP <- cbind(MAP, TP.p[,i]) colnames(tmp.TP) <- c("CHR", "SNP", "CM", "BP", "P") qq(tmp.TP$P, main = paste0("D", i)) mtext(LETTERS[i-10], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() } } ##RR for (i in 1:20){ if(i < 10){ pdf("DataAndCode/Figures/FigS8.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 1:10){ tmp.RR <- cbind(MAP, RR.p[,i]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "P") qq(tmp.RR$P, main = paste0("D", i)) mtext(LETTERS[i], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() }else{ pdf("DataAndCode/Figures/FigS9.pdf", h=10, w=7, useDingbats = F, pointsize = 10) par(mar=c(3,3,2,1), mgp=c(1.8,0.5,0)) nf = layout(rbind(c(1,2), c(3,4), c(5,6), c(7,8), c(9,10))) for (i in 11:20){ tmp.RR <- cbind(MAP, RR.p[,i]) colnames(tmp.RR) <- c("CHR", "SNP", "CM", "BP", "P") qq(tmp.RR$P, main = paste0("D", i)) mtext(LETTERS[i-10], 2, adj=3.5, las=1, padj=-10, font=1, cex=0.8) } dev.off() } } ```