How to overlay qq plots
4
0
Entering edit mode
7.0 years ago

Hi, I'd like to overlay two qq plots from GWAS (before and after adjustment for genomic inflation) in one image as in the following link.

https://openi.nlm.nih.gov/detailedresult.php?img=PMC3093379_pone.0019416.g003&req=4

I used R packages qqman and ggplot2, but both packages seem to not support this function. So, can someone help me with this?

Thanks,

qqplot • 8.0k views
1
Entering edit mode
7.0 years ago

Following on from @JM88's answer, here's the same thing in ggplot2.

# Set seed and make toy data
set.seed(666)
df <- data.frame(X = rnorm(20), Y = rnorm(20), Z = rnorm(20))

# Create plot using ggplot2
library(ggplot2)
gg <- ggplot(df) +
geom_point(data = df, aes(x = X, y = Y), colour = "red") +
geom_point(data = df, aes(x = X, y = Z), colour = "black") +
theme_bw()
print(gg)


0
Entering edit mode
7.0 years ago

You can use the qqplot function and then points

## to make this reproducible
set.seed(666)
df <- data.frame(X = rnorm(20), Y = rnorm(20), Z = rnorm(20))

## qq plot
qqplot(df$X, df$Y, xlim = range(df), ylim = range(df) ,xlab = "X", ylab = "Y and Z", pch =19)

points(sort(df$X), sort(df$Z), col = "darkred", pch =19)

## diagonal line
abline(0,1)

0
Entering edit mode
7.0 years ago
willgilks ▴ 360

I've run the code for the two answers above, and the plots do not look the same, because the R qqplot function applies a transformation to the data.

The code on this page generates a qqplot like the one in the example, although it's using lattice, not ggplot http://genome.sph.umich.edu/wiki/Code_Sample:_Generating_QQ_Plots_in_R

I was trying to work out how to calculate and plot the 95%CI on ggplot a while ago. The code is below but there's clearly something wrong with the plotted interval values.

p.values = rnorm(100, mean = 0.5, sd = 0.1)
obs1 = -log10 (sort(p.values))
exp1 = -log10 (1:length (obs1)/length (obs1))
r1 = tbl_df (data.frame(obs1, exp1))

c95 = rep(0,length(exp1))
c05 = rep(0,length(exp1))

## the jth order statistic from a uniform(0,1) sample has a beta(j,n-j+1) distribution (Casella & Berger, 2002)

for(i in 1:length(exp1)){
c95[i] <- qbeta(0.95,i,length(exp1)-i+1)
c05[i] <- qbeta(0.05,i,length(exp1)-i+1)
}

r1 = cbind(r1, c05, c95)

ggplot (r1, aes(obs1, exp1)) +
geom_point() +
geom_abline(intercept = 0, slope = 1, alpha = .5) +
geom_smooth(method = "lm", se = FALSE) +
geom_ribbon(aes (ymin = c05, ymax = c95, alpha = 0.5)) +
theme_bw()

0
Entering edit mode
5.0 years ago
Phoenix Mu ▴ 100

This is the function I used. It is modified from the original qq() function from qqman package.

qqplot <- function(pvector, col = "black", add = F, ylim = NULL){
expectedP <- -log10(ppoints(length(pvector)))
observedP <- -log10(sort(pvector, decreasing = F))
plot(x = expectedP, y = observedP, col = col,
xlab = expression(Expected ~ ~-log[10](italic(p))),
ylab = expression(Observed ~ ~-log[10](italic(p))),
pch = 16, cex = 0.7, ylim = NULL)
abline(0, 1, col = "red")
}else{
points(x = expectedP, y = observedP, col = col,
xlab = expression(Expected ~ ~-log[10](italic(p))),
ylab = expression(Observed ~ ~-log[10](italic(p))),
pch = 16, cex = 0.7, ylim = NULL)
}


Use add=T if you want to overlay another qqplot to an existing one.