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) )
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)) #38Compare 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)
#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 #191Plot 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)
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). 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
}
}
# add threshold lines
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))
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()
}
}