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

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

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

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.