version 0.6.7

cran · Oct 22, 2012 · e738059 · e738059
1 parent 50123b7
commit e738059
Show file tree

Hide file tree

Showing 16 changed files with 531 additions and 45 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,17 +1,17 @@
 Package: qtlhot
-Version: 0.6.6
-Date: 2012-10-02
+Version: 0.6.7
+Date: 2012-10-22
 Author: Elias Chaibub Neto <echaibub@hotmail.com> and Brian S Yandell
         <byandell@wisc.edu>
 Title: Inference for QTL Hotspots
 Description: Functions to infer co-mapping trait hotspots and causal
         models
 Maintainer: Brian S. Yandell <yandell@stat.wisc.edu>
-Depends: R (>= 2.10), stats, qtl, lattice, corpcor, mnormt
+Depends: stats,qtl,lattice,corpcor,mnormt
 LazyLoad: yes
 LazyData: yes
 License: GPL (>= 2)
 URL: http://www.stat.wisc.edu/~yandell/statgen
-Packaged: 2012-10-02 17:33:49 UTC; yandell
+Packaged: 2012-10-22 16:25:35 UTC; yandell
 Repository: CRAN
-Date/Publication: 2012-10-02 18:22:06
+Date/Publication: 2012-10-23 06:33:22
diff --git a/MD5 b/MD5
@@ -1,12 +1,13 @@
-60e0d5a0997d40ce7e0042f0a1f6b589 *DESCRIPTION
-57e7ffef8bcbcf74d0fd2f24c96a94a2 *NAMESPACE
+f5a067410f155092f6ce9bb65547cd21 *DESCRIPTION
+812ed604c58795ede2862ce6a4e54dba *NAMESPACE
 b933207ad0c19257762b995ec2ac4730 *R/add.phenos.R
 103b12e46d558b8b37c3eb8136e76d69 *R/big.R
-bc5bd01be1400c58fbdf5bcf7306db97 *R/cmst.R
-0bcb4f8b6009685ae0d63d2139a7f70b *R/highlod.R
-74ea08f27428cbcd32935c1cd343b898 *R/hotperm.R
+382b749f1b777a0b8d3db26b8bb8dda4 *R/cmst.R
+e81baa5d1bd526eecfcf013267ca42c9 *R/highlod.R
+dec7898ba84d5b27de5d678f82e07fb7 *R/hotperm.R
 6e75c18121af3d2e3b1c133b06de311e *R/hotsize.R
 9d957d0691daafd8c8c3cae7d42c84be *R/max.R
+8afdc7d138adfa1284bd90cf63b8e605 *R/multiple_testing_sim_functions.R
 2691e0df2e1c4ad991982cee548e06db *R/neqtl.R
 5a0221e370e2ec06215b5901605af28d *R/parallel.R
 5973406c7885b09e6c49f3df11081bdf *R/plotcmst.R
@@ -22,12 +23,12 @@ dc2f482c9166cfc488dd0c0de6ee1642 *R/utilcmst.R
 39bf135b5eb185681d74d06c18549a1c *R/utilhot.R
 d41c35875c29eff957f88d263399c6ad *R/westwu.R
 d41d8cd98f00b204e9800998ecf8427e *README
-c11d97bf0278ba75468f817255c876a2 *data/CMSTCross.RData
+d9012f4bb4cc97d0cd07cf3bea96ca85 *data/CMSTCross.RData
 89bff8163dbfdb0b2e5fe28276ecd4fa *inst/deprecated.R
 d6ba35863c52c32349ce63a713b56652 *inst/doc/cmst.Rnw
-e1ff7fd37d28d308287187fd47f46d6e *inst/doc/cmst.pdf
+7eb137a2e6470471dccf4bf1cae0e85e *inst/doc/cmst.pdf
 3b4ae999433e50910389135b1d7eb850 *inst/doc/qtlhot.Rnw
-6eea94f85f5422d1d0c9275551bbab49 *inst/doc/qtlhot.pdf
+bc493486930c64899bd4508315aa3fe8 *inst/doc/qtlhot.pdf
 11704113157929b53b0e2c5ba015ad3e *inst/notes/cmst.txt
 52d8cbf0da7b5bdc67b82665b7e8c4c1 *inst/notes/code.txt
 c15b74d7121a2ddbe21306b65ccd3b64 *inst/notes/perm.txt
@@ -43,16 +44,16 @@ a5c3774fa0c5b298161fd419ec5bcb50 *inst/parallel/SOAR.post.R
 0eac83f946b60f56fe187ddcd4a7760f *inst/parallel/cross.RData
 8b9712f8ce05f9005356e0b06798357c *inst/parallel/params.txt
 be0c94b6582c2f5f503730bb94833b21 *inst/parallel/runparallel.sh
-0d8ae29e6c25592f1fb7937975a0be9b *man/CMSTtests.Rd
-19953211b4de88ac4c743b769fa4207a *man/GetCandReg.Rd
+3ef4b3690d7b3436b5895a5a6fbf396f *man/CMSTtests.Rd
+b26b709dc559c8d77ba4f51ba19d3156 *man/GetCandReg.Rd
 a0d9688690f56d672a90cec7117b0a65 *man/GetCommonQtls.Rd
-e16a5f2844aa0b97f10e69e6543d9bd5 *man/PrecTpFpMatrix.Rd
-cdb8e3173159732116f20ff08da18515 *man/SimCrossCausal.Rd
+5b700880f6ffe2b9d90946a4731a22b3 *man/PrecTpFpMatrix.Rd
+efb799895e471cdd70bc2d2f30545c7a *man/SimCrossCausal.Rd
 692111bb7328f17b403106ab8dde4907 *man/add.phenos.Rd
 e8cdb877d98ca4a5f27ad7fd42cd79e1 *man/filter.threshold.Rd
 fb6907ded1db3d56c14d94389360759b *man/highlod.Rd
-7ebe43fede89e9a04ca6a81e0e17b8d9 *man/hotperm.Rd
-a7faaf43638df9f53c350a4aedfb5d84 *man/hotsize.Rd
+2f25cffd80841479c6b010ec9eeb2145 *man/hotperm.Rd
+7780b3bf6da7bac3831b7dbdb600f72b *man/hotsize.Rd
 c7ca914d7fc9cf321a7e07750fcf1f32 *man/parallel.qtlhot.Rd
 3da742be31fe06f70772508c1a51a2d6 *man/sim.hotspot.Rd
 42d951c0229a8e83cc216a2e834e204f *man/westwu.Rd

diff --git a/NAMESPACE b/NAMESPACE
@@ -60,5 +60,6 @@ export(GetCoMappingTraits)
 export(FitAllTests)
 export(JoinTestOutputs)
 export(PrecTpFpMatrix)
+export(p.adjust.np)
 
 ## export(mySimulations)
diff --git a/R/cmst.R b/R/cmst.R
@@ -1,4 +1,5 @@
 ## Modify to use names element of highlod object?
+## Fix so call with length(pheno2) == 1 and > 1 give same results; see JoinTestOutputs.
 CMSTtests <- function(cross, 
                       pheno1, 
                       pheno2,
@@ -917,15 +918,15 @@ FitAllTests <- function(cross, pheno1, pheno2, Q.chr, Q.pos, verbose = TRUE)
                      NULL, NULL, NULL, NULL, "all", "both", FALSE)
 
   nms <- pheno2
+  ntests <- length(pheno2)
   out$pvals.cit <- matrix(NA, ntests, 2, dimnames = list(nms, c("pval.1", "pval.2")))
 
-  ntests <- length(pheno2)
   for(k in 1 : ntests) {
     cit.mar <- find.marker(cross, Q.chr, Q.pos)
     LL <- pull.geno(cross)[, cit.mar]
     GG <- cross$pheno[, pheno1]
     TT <- cross$pheno[, pheno2[k]]
-    aux2 <- try(CitTests(LL, GG, TT), silent = TRUE)
+    aux2 <- try(qtlhot:::CitTests(LL, GG, TT), silent = TRUE)
     if(class(aux2) != "try-error") {
       out$pvals.cit[k,] <- aux2
     }
@@ -1197,6 +1198,78 @@ counts <- function(out, alpha, method=c("aic","bic","cit",
   ###
   output
 }
+##############################################################################
+PerformanceSummariesKo <- function(alpha, nms, val.targets, all.orfs,
+                                   tests, cis.index)
+{
+  rnms <- c("aic", "bic", "j.bic", "p.bic", "np.bic", "j.aic", "p.aic", "np.aic", "cit")
+  nt <- length(nms)
+  TP <- FP <- TN <- FN <- NC <- data.frame(matrix(0, 9, nt))
+  tar <- rep(NA, nt)
+  names(TP) <- names(FP) <- names(TN) <- names(FN) <- names(NC) <- nms
+  row.names(TP) <- row.names(FP) <- row.names(TN) <- row.names(FN) <- row.names(NC) <- rnms
+  Causal <- NotCausal <- vector(mode = "list", length = 9)
+  for (k in 1 : nt) {
+      aux <- which(tests[[11]][,1] == nms[k])
+      aux.nms <- tests[[11]][aux, 2]
+      tar[k] <- length(aux)
+      for (i in 2 : 3) {
+        aux.rank <- apply(tests[[i]][aux, 1:4, drop = F], 1, rank)
+        aux.best <- apply(aux.rank, 2, function(x) which.min(x))
+        aux.index <- as.numeric(which(aux.best == 1))
+        Causal[[i - 1]] <- aux.nms[aux.index]
+        NotCausal[[i - 1]] <- aux.nms[-aux.index]
+      }
+      for (i in 4 : 9) {
+        Causal[[i - 1]] <- aux.nms[which(tests[[i]][aux, 1] <= alpha)]
+        NotCausal[[i - 1]] <-
+          aux.nms[c(which(tests[[i]][aux, 2] <= alpha),
+                    which(tests[[i]][aux, 3] <= alpha),
+                    which(tests[[i]][aux, 4] <= alpha))]
+      }
+      Causal[[9]] <-
+        aux.nms[which(tests[[10]][aux, 1] <= alpha & tests[[10]][aux, 2] > alpha)]
+      NotCausal[[9]] <-
+        aux.nms[c(which(tests[[10]][aux, 1] > alpha & tests[[10]][aux, 2] <= alpha),
+                  which(tests[[10]][aux, 1] >= alpha & tests[[10]][aux, 2] >= alpha))]
+      val <- val.targets[[match(nms[k], names(val.targets))]]
+      not.val <- all.orfs[-match(unique(c(nms[k], val)), all.orfs)]
+      for (i in 1 : 9) {
+        TP[i, k] <- length(!is.na(intersect(Causal[[i]], val)))
+        FP[i, k] <- length(!is.na(intersect(Causal[[i]], not.val)))
+        TN[i, k] <- length(!is.na(intersect(NotCausal[[i]], not.val)))
+        FN[i, k] <- length(!is.na(intersect(NotCausal[[i]], val)))
+      }
+      for (i in 4 : 9) {
+        NC[i - 1, k] <- length(c(which(tests[[i]][aux, 1] > alpha),
+                             which(tests[[i]][aux, 2] > alpha),
+                             which(tests[[i]][aux, 3] > alpha),
+                             which(tests[[i]][aux, 4] > alpha)))
+      }
+      NC[9, k] <- length(c(which(tests[[10]][aux, 1] < alpha),
+                           which(tests[[10]][aux, 2] < alpha)))
+  }
+  tp <- apply(TP, 1, sum)
+  fp <- apply(FP, 1, sum)
+  tn <- apply(TN, 1, sum)
+  fn <- apply(FN, 1, sum)
+  nc <- apply(NC, 1, sum)
+  prec <- tp/(tp + fp)
+  tpr <- tp/(tp + fn)
+  fpr <- fp/(fp + tn)
+  overall.1 <- data.frame(prec, tp, fp, tpr, fpr, tn, fn, nc)
+  tp <- apply(TP[, cis.index], 1, sum)
+  fp <- apply(FP[, cis.index], 1, sum)
+  tn <- apply(TN[, cis.index], 1, sum)
+  fn <- apply(FN[, cis.index], 1, sum)
+  nc <- apply(NC[, cis.index], 1, sum)
+  prec <- tp/(tp + fp)
+  tpr <- tp/(tp + fn)
+  fpr <- fp/(fp + tn)
+  overall.2 <- data.frame(prec, tp, fp, tpr, fpr, tn, fn, nc)
+  list(overall.1, overall.2, tar)
+}
+
 ##############################################################################
 performance.summaries.cmst <- function(out, model, alpha=0.05, method)
 {
@@ -1353,6 +1426,19 @@ JoinTestOutputs <- function(comap)
 
   for (k in 2 : length(comap)) {
     load(paste("output_ko_validation", reg.nms[k], "Rdata", sep="."))
+
+    if (length(out) != 10) { ## CMSTtests if length(pheno2) == 1
+      tmp <- t(out$Z.aic)
+      AIC.stats <- c(out$AICs, tmp[!is.na(tmp)])
+      tmp <- t(out$Z.bic)
+      BIC.stats <- c(out$BICs, tmp[!is.na(tmp)])      
+      out <- list(R2s = out$R2, AIC.stats = AIC.stats, BIC.stats = BIC.stats, 
+                  pvals.j.BIC = out$pvals.j.BIC, pvals.p.BIC = out$pvals.p.BIC, 
+                  pvals.np.BIC = out$pvals.np.BIC,  pvals.j.AIC = out$pvals.j.AIC, 
+                  pvals.p.AIC = out$pvals.p.AIC,  pvals.np.AIC = out$pvals.np.AIC,
+                  pvals.cit = out$pvals.cit)
+    }
+
     for (i in 1 : 10) {
       join.out[[i]] <- rbind(join.out[[i]], out[[i]])
     }
@@ -1362,6 +1448,16 @@ JoinTestOutputs <- function(comap)
  join.out
 }
 ##############################################################################
+p.adjust.np <- function(tests, method = "BH")
+{
+  for(test.name in paste("pvals.np", c("AIC","BIC"), sep = "."))
+    for (i in seq(ncol(tests[[test.name]])))
+      tests[[test.name]][, i] <- p.adjust(tests[[test.name]][, i], method = method)
+
+  tests
+
+}
+##############################################################################
 GetCis <- function(x, window = 10) {
   xx <- x[x[, 2] == x[, 4],]
   xx <- xx[abs(xx[, 3] - xx[, 5]) <= window, ]
@@ -1409,7 +1505,8 @@ GetCisCandReg <- function(highobj, cand.reg, lod.thr = NULL)
   out <- data.frame(cand.reg, peak.pos)
   is.cis <- (out$phys.pos >= out$peak.pos.lower &
              out$phys.pos <= out$peak.pos.upper)
-  out <- out[is.cis,, drop = FALSE]
+  ## Keep cis traits, but leave off peak.chr (since it == phys.chr).
+  out <- out[is.cis, -4, drop = FALSE]
   if(nrow(out))
     attr(out, "cis.index") <- match(out[, 1], highobj$names)
   out

diff --git a/R/highlod.R b/R/highlod.R
@@ -99,14 +99,7 @@ highlod <- function(scans, lod.thr = 0, drop.lod = 1.5,
 print.highlod <- function(x, ...) print(summary(x, ...))
 summary.highlod <- function(object, ...)
 {
-  cat("frequency of high LOD entries by chromosome:\n")
-  print(addmargins(table(object$chr.pos$chr[object$highlod$row])))
-  cat("\nfrequency of high LOD entries per phenotype:\n")
-  tbl <- table(table(object$names[object$highlod$pheno]))
-  if(!is.null(object$names))
-    names(tbl) <- object$names[as.numeric(names(tbl))]
-  print(tbl)
-  invisible()
+  summary(hotsize(object, ...))
 }
 plot.highlod <- function(x, ..., quant.level = NULL, sliding = FALSE)
 {

diff --git a/R/hotperm.R b/R/hotperm.R
@@ -37,6 +37,10 @@ hotperm <- function(cross, n.quant, n.perm, lod.thrs, alpha.levels, drop.lod = 1
   n.phe <- nphe(cross)
   n.ind <- nind(cross)
   n.quant <- min(n.quant, n.phe)
+
+  tmp <- table(sapply(cross$pheno, class))
+  if(length(tmp) > 1 | names(tmp)[1] != "numeric")
+    stop("all phenotypes in cross object must be numeric")
 
   s.quant <- seq(n.quant)
   quants <- 1 - (s.quant - 1) / n.phe