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). 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.

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) )

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).

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) )

#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.

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)

##
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))

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). The modified function plots the $$|\beta|$$ rather than $$-log_{10}(p)$$.

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 0<P<=1 d = d[order(d$CHR, d$BP), ] # sets colors based on colors argument. colors <- rep(colors,max(d$CHR))[1:max(d$CHR)] # sets the maximum value on the y axis if (ymax == "max") ymax<-ceiling(max(d$Beta))

# 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))
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
# creates a blank plot
with(d, plot(pos, Beta, ylim = c(0,ymax), ylab = expression("|" ~ beta ~ "|"), 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, Beta, 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), ]
icol=1
for (i in unique(d.annotate$CHR)) { with(d.annotate[d.annotate$CHR==i, ], points(pos, Beta, col = "red", cex=0.5, pch = 20, ...))
icol = icol+1
}
}

if (!is.null(suggestiveline)) abline(h=suggestiveline, col="blue")
if (!is.null(genomewideline)) abline(h=genomewideline, col="red")
}

manhattan <- function(dataframe, colors=c("gray10", "gray50"), ymax="max", xaxis.cex=1, yaxis.cex = 1, limitchromosomes=1:23, suggestiveline=-log10(1e-5), genomewideline=-log10(5e-8), annotate=NULL, Title, ...) {

d=dataframe

#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) & "P" %in% names(d))) stop("Make sure your data frame contains columns CHR, BP, and P")

# 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 0<P<=1 d=subset(na.omit(d[order(d$CHR, d$BP), ]), (P>0 & 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") } 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 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))
}
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))
}
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))
}