How to overlay qq plots
4
0
Entering edit mode
7.8 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.6k views
ADD COMMENT
1
Entering edit mode
7.8 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)

https://github.com/AndrewSkelton/Biostars-Answers/blob/master/p237051.R

ADD COMMENT
0
Entering edit mode
7.8 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)

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

## diagonal line
abline(0,1)

Credits to: http://stackoverflow.com/questions/20299374/qq-plot-more-than-two-data

ADD COMMENT
0
Entering edit mode
7.8 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()
ADD COMMENT
0
Entering edit mode
5.8 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))
  if (add == 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.

ADD COMMENT

Login before adding your answer.

Traffic: 2413 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6