We need to fix the Q-value calculation. This goes awry when there are a lot - or only - very significant p-values, or some NAs. Or too little values. Currently we use the Benjamini - Hochberg FDR correction, which is fine (a bit less conservative) and quite similar to Q-value (Storey & Tibshirani, which is the least conservative).
It pertains this part of the QTL_QC.R script:
cat("\n* Least conservative correction: Storey & Tibshirani correction...\n")
### Storey & Tibshirani correction - Least conservative
### references:
### - http://en.wikipedia.org/wiki/False_discovery_rate
### - http://svitsrv25.epfl.ch/R-doc/library/qvalue/html/qvalue.html
### Requires a bioconductor package: "qvalue"
if(opt$resulttype == "NOM") {
#RESULTS$Q = qvalue(RESULTS$Nominal_P)$qvalues # original code
RESULTS$Q = "Not calculated: throws an error when p-value is infinite or NA. NEED FIXING"
} else if(opt$resulttype == "PERM") {
#RESULTS$Q = qvalue(RESULTS$Approx_Perm_P)$qvalues # original code
RESULTS$Q = ifelse(RESULTS$Approx_Perm_P > 0, qvalue(RESULTS$Approx_Perm_P)$qvalues, "NA")
} else {
cat ("\n\n*** ERROR *** Something is rotten in the City of Gotham; most likely a typo. Double back, please.\n\n",
file=stderr()) # print error messages to stder
}
# RESULTS$Q = "Currently not calculated due to an issue with the qvalue() package."
We need to fix the Q-value calculation. This goes awry when there are a lot - or only - very significant p-values, or some NAs. Or too little values. Currently we use the Benjamini - Hochberg FDR correction, which is fine (a bit less conservative) and quite similar to Q-value (Storey & Tibshirani, which is the least conservative).
It pertains this part of the
QTL_QC.Rscript: