# An R script for creating a QQ-plot of the test statistics in the PRESTO p-value file.
#
# Example: If you want to create a QQ-plot using the 3rd column of a PRESTO output p-value file
# called "presto.pval", perform the following 3 steps:
#
# 1) Set the working directory of R to a directory containing this file and the PRESTO p-value file
#    (see the R "setwd" and "getwd" commands)
# 2) At the R prompt enter: source("qq.plot.r")
# 3) At the R prompt enter: presto.qq.plot("presto.pval", 3)
#

#---------------------------------------------------------------------------------------------------

# presto.qq.plot: function for creating QQ-plots from PRESTO output p-value file
#
# Parameters:
#
# filename:   filename of PRESTO output pvalue (.pval) file.
#
# columns:    column number (counting from 1) of the p-value file to use in the QQ-plot.  The column number 
#             should correspond to one of the following five columns: "allelic_X2", "rec_X2", "dom_X2", 
#             "overdom_X2", or "trend_X2".  At least one of these five columns will be present in 
#             PRESTO's output p-value (.pval) file.
#
# max.chisq:  upper bound for observed chi-square test statistics (default = 50.0).  All chi-square test
#             statistics greater than max.chisq are represented by triangles with y-coordinate equal to max-chisq.
#
# Reference: Pawitan, Yudi (2001) "In all Likelihood: Statistical modelling and inference using likelihood."
#            Oxford: Clarendon Press.  pp 91-92.
#
presto.qq.plot <- function(filename, column, max.chisq=50.0) {
    d <- read.table(filename, header=T)
    statistic <- names(d)[column]
    observed <- sort(d[[column]])
        
    n <- length(observed)
    u <- ((1:n) - 1/3) / (n + 1/3)
    expected <- qchisq(u, df=1)
    
    indices <- observed<max.chisq
    sum.indices <- sum(indices)
    if (sum.indices==0) {
        stop("all test statistics > max.chisq")
    }
    plot(expected[indices], observed[indices], pch=20, xlab="Expected chi-square value", ylab="Observed test statistic",
        xlim=c(0.0, max(expected)), ylim=c(0.0, max.chisq))
    if (sum.indices < length(indices)) {
        points(expected[!indices], rep(max.chisq, length(expected[!indices])), pch=24)
    }
    title(paste("QQ-plot for", statistic, "statistic"))
    lines(expected, expected)
    grid()  
}