--- a
+++ b/createBulkAtlantisFlock.R
@@ -0,0 +1,491 @@
+library(gridExtra)
+library(foreach)
+library(atlantis)
+
+createBulkAtlantisFlock <- function(analysis.description, run_mode, predMat.file, targetMat.file, anno.file, featureSetToUse, kMinNumSamples=50, seed=300, batch.size=1, predFeaturesPerTarget=NULL, fitControlSettings="default", config.summary=NULL, predictOnly=NULL, permutation.count=NULL, only.permutations=F, 
+  code.dir=".",
+  test.fraction=NA,
+  predFeatureRegExps=NULL) {
+  # analysis.description: a description which will be saved in the run directory
+  # run_mode: ignored
+  # mustHave: a vector of dataset ids which must have a value present to be included in the prediction matrix
+  # niceToHave: a vector of dataset ids to include in the prediction matrix, but these may have missing values in the prediction matrix
+  # featureSetToUse: one of "all", "perTarget", "perTargetAndSelf" or "single" determining how to subset the prediction matrix when
+  #                  for a given target.
+  # predFeaturesPerTarget: if featureSetToUse = "perTarget" or "perTargetAndSelf" then this should be a dataframe with two columns: c("target", "partner") which will be used to look up which genes to use as predicitive features based on a given target
+
+  if (featureSetToUse == "perTarget" || featureSetToUse == "perTargetAndSelf") {
+	  stopifnot(!is.null(predFeaturesPerTarget) && is.data.frame(predFeaturesPerTarget) && all(colnames(predFeaturesPerTarget) == c("target", "partner")))
+  } else if(featureSetToUse == "perGroup") {
+     stopifnot(!is.null(predFeaturesPerTarget) && is.data.frame(predFeaturesPerTarget) && all(colnames(predFeaturesPerTarget) == c("group", "target", "group.name")))
+  }
+
+  # flock_run_dir is the name of the directory of where the results will be written
+  # the analysis.ID is the trailing dir name
+  analysis.ID <- basename(flock_run_dir)
+  output.dir <- paste(flock_run_dir, '/results/', sep='')
+  dir.create(output.dir)
+
+  # cache the functional related genes if we'll need them
+  predFeaturesPerTarget.file <- NULL
+  if(featureSetToUse == "perTarget" || featureSetToUse == "perTargetAndSelf" || featureSetToUse == "perGroup") {
+    stopifnot(!is.null(predFeaturesPerTarget))
+
+    predFeaturesPerTarget.file <- paste(output.dir, '/featuresPerTarget.Rdata', sep="")
+    save(predFeaturesPerTarget, file=predFeaturesPerTarget.file)
+  }
+
+  # load the target matrix to determine the targets
+  load(targetMat.file);
+  stopifnot(is.matrix(targetMat) || is.data.frame(targetMat))
+  targets <- colnames(targetMat);
+  if(is.null(permutation.count)) {
+    permutation.count <- length(targets)
+  }
+
+  if(featureSetToUse == "perGroup") {
+    # then we want to create a target per member in our predFeaturesPerTarget table
+    with.valid.targets <- predFeaturesPerTarget[paste("GS_",predFeaturesPerTarget$target,sep='') %in% targets,,drop=F]
+    targets <- lapply(seq(nrow( with.valid.targets)), function(i) { 
+	x<-with.valid.targets[i,]; 
+        list(ID=paste("GS_",as.character(x$target),sep=''), 
+             gene=as.character(x$target), 
+             group=as.character(x$group), 
+             output.prefix=paste(x$group,'-',x$target,sep=''),
+             group.name=as.character(x$group.name)) 
+        } )
+    stopifnot(! ("GS_NA" %in% sapply(targets, function(t){t$ID})))
+  } else {
+    # filter out targets which have a single output value
+    print("before filtering out singles")
+    print(length(targets))
+    is.bad.target <- sapply(seq(dim(targetMat)[[2]]), function(i) { t<-targetMat[,i]; t<-t[!is.na(t)]; length(unique(t))==1 } )
+    targets <- targets[!is.bad.target]
+    targets <- lapply(targets, function(t) { list(ID=t, output.prefix=t) })
+    print("after filtering out targets with a single target value")
+    print(length(targets))
+    stopifnot(grep("GS_NA", targets)!=0)
+  }
+
+  common <- list( 
+    fitControlSettings=fitControlSettings,
+    targetMat.file=targetMat.file,
+    predMat.file=predMat.file,
+    anno.file=anno.file,
+    output.dir=output.dir,
+    featureSetToUse=featureSetToUse,
+    kMinNumSamples=kMinNumSamples,
+    predFeaturesPerTarget.file=predFeaturesPerTarget.file,
+    seed=seed,
+    quality.dist.file=paste(output.dir, "null_quality_distribution.txt", sep=""),
+    config.summary=config.summary,
+    predictOnly=predictOnly,
+    permutation.count=permutation.count,
+    only.permutations=only.permutations,
+    batch.size=batch.size,
+    predFeatureRegExps=predFeatureRegExps
+    )
+
+  sources <- c("createBulkAtlantisFlock.R")
+
+  common$sources <- sources
+
+  # add in runs for each permuted version used for the null
+  permutations <- sample(targets, permutation.count, replace=T)
+  for(i in seq_along(permutations)) { 
+    permutations[[i]]$permute.seed <- i + seed
+  }
+
+  if(common$only.permutations) {
+    tasks <- permutations
+  } else {
+    tasks <- c(targets, permutations)
+  }
+
+  tasks <- split(tasks, as.integer(seq(length(tasks))/common$batch.size))  
+
+  if(!is.na(test.fraction)) {
+    task.count <- as.integer(ceiling(test.fraction * length(tasks)))
+    cat("Only running", task.count, "tasks out of", length(tasks), "because test.fraction was set to", test.fraction, "\n")
+    tasks <- tasks [ seq(task.count) ]
+  }
+
+  flock.run(tasks, task_function='per_target', gather_function='gather', sources=common$sources, flock_common_state=common)
+}
+
+prev.proc.time <- NULL
+
+time.checkpoint <- function(label) {
+#  now <- proc.time()
+#  if(!is.null(prev.proc.time)) {
+#    cat(sprintf("%s: %f\n", label, (now["elapsed"] - prev.proc.time["elapsed"])))
+#  } else {
+#    cat(sprintf("%s: first\n", label))
+#  }
+#  prev.proc.time <<- now
+}
+
+# given a gene symbol, returns a list of functionaly related genes based on dataframe which was persisted
+get.related.genes <- function(predFeaturesPerTarget, targetToLookup, target, include.self) {
+  matches <- as.character(predFeaturesPerTarget$target) == targetToLookup
+  partners <- as.character(predFeaturesPerTarget[matches, "partner"])
+  if(include.self) {
+    ret <- unique(c(target, partners))
+  } else {
+    ret <- setdiff(unique(partners), target)
+  }
+
+  return(ret)
+#  sapply(ret, getFeatureNameFromFeatureID)
+}
+
+get.gene.group <- function(predFeaturesPerTarget, groupToLookup, target, include.self) {
+  partners <- as.character(predFeaturesPerTarget[predFeaturesPerTarget$group == groupToLookup, "target"])
+  if(include.self) {
+    ret <- unique(c(target, partners))
+  } else {
+    ret <- setdiff(unique(partners), target)
+  }
+  return(ret)
+#  sapply(ret, getFeatureNameFromFeatureID)
+}
+
+determine.feature.genes <- function(target, featureSetToUse, predFeaturesPerTarget.file) {
+  if(!is.null(predFeaturesPerTarget.file)) {
+    load(predFeaturesPerTarget.file)
+  }
+
+  if (featureSetToUse == "all") {
+    genes <- NULL
+  } else if (featureSetToUse == "perTarget" || featureSetToUse == "perTargetAndSelf") {
+    targetName <- getFeatureNameFromFeatureID(target$ID)
+    stopifnot(!is.null(targetName))
+    genes <- get.related.genes(predFeaturesPerTarget, targetName, targetName, featureSetToUse == "perTargetAndSelf")
+    cat("found", length(genes), "gene names to be used as features for", targetName, "\n")
+  } else if (featureSetToUse == "perGroup") {
+    stopifnot(all(colnames(predFeaturesPerTarget) == c("group", "target", "group.name")))
+    group <- target$group
+    stopifnot(! is.null(group))
+    stopifnot(! is.null(target$gene))
+    genes <- get.gene.group(predFeaturesPerTarget, group, target$gene, FALSE)
+  } else {
+    stopifnot(featureSetToUse ==  "single")
+    targetName <- getFeatureNameFromFeatureID(target$ID)
+    genes <- targetName
+    stopifnot(!is.null(genes))
+  }
+
+  return(genes)
+}
+
+per_target <- function(per.task.state, common.state, output.file, run) {
+  library(foreach)
+
+  flock_run_dir <- run$dir
+
+  analysis.dir <- basename(flock_run_dir);
+  analysis.ID <- analysis.dir;
+  output.dir <- flock_common_state$output.dir
+  predMat.file <- flock_common_state$predMat.file
+  targetMat.file <- flock_common_state$targetMat.file
+  anno.file <- flock_common_state$anno
+  kMinNumSamples <- flock_common_state$kMinNumSamples
+  predFeaturesPerTarget.file <- flock_common_state$predFeaturesPerTarget.file
+  config.summary <- flock_common_state$config.summary
+  predFeatureRegExps <- flock_common_state$predFeatureRegExps
+
+  fitControlSettings = flock_common_state$fitControlSettings
+
+  # the following figures out what are the 100 most similar genes to target gene and uses only them
+  # for prediction.
+  ####
+  targetMat.file <- flock_common_state$targetMat.file
+  load(targetMat.file)
+
+  featureSetToUse <- flock_common_state$featureSetToUse
+
+  fit.single.target <- function(target) {
+    # allow target to be specified as either a string or a list of (ID, group)
+    if(is.character(target)) {
+      target <- list(ID=target, targetUsedToLookupRelated=target, permute.seed=NULL)
+    }
+
+    targetID <- target$ID
+    permute.seed <- target$permute.seed
+
+    print(targetID)
+    cat(sprintf("Generating model for %s\n", targetID))
+    set.seed(flock_common_state$seed)
+
+    if (is.numeric(targetID)) {
+      targetID <- colnames(targetMat)[targetID]
+    }
+
+    stopifnot(length(grep("NO_CURRENT", targetID)) == 0)
+
+    permuteRows <- function(seed, targetMat) {
+      saved <- .Random.seed
+      set.seed(seed)
+      permutedTargetMat <- targetMat
+      for(i in seq(ncol(targetMat))) { 
+        permutedTargetMat[,i] = sample(targetMat[,i]) 
+      }
+      .Random.seed <- saved
+      permutedTargetMat
+    }
+
+    if(featureSetToUse == "all") {
+      predFeatureRegExps <- ".*"
+      genes <- NULL
+    } else {
+      genes <- determine.feature.genes(target, featureSetToUse, predFeaturesPerTarget.file)
+    }
+
+    output.prefix <- target$output.prefix
+    group.name <- target$group.name
+    if(is.null(group.name)) {
+      group.name <- NA
+    } else {
+      config.summary <- c(config.summary, paste("group:", group.name))
+    }
+
+    # if this is for the null distribution
+    if(!is.null(permute.seed)) {
+      permutedTargetMat <- permuteRows(permute.seed, targetMat)
+
+      res <- runATLANTIS(
+        analysis.ID=analysis.ID,
+        targetID=targetID,
+        output.prefix=paste(output.prefix, "-NULL-",permute.seed, sep=''),
+        makePlot=FALSE,
+        output.dir = paste(dirname(output.dir), '/temp', sep=''),
+        fitControlSettings = fitControlSettings,
+        additionalFeatures=NULL,
+        predFeatureNamesToUse=genes,
+        predFeatureRegExps=predFeatureRegExps,
+        predMat.file = predMat.file,
+        targetMat = permutedTargetMat,
+        anno.file=anno.file,
+        kMinNumSamples=kMinNumSamples,
+        save.params=F,
+        save.model=F,
+        save.featureData=F)
+    } else {
+      res <- runATLANTIS(
+        analysis.ID=analysis.ID,
+        targetID=targetID, 
+        makePlot=TRUE,
+        output.dir = output.dir,
+        output.prefix=output.prefix,
+        fitControlSettings = fitControlSettings,
+        additionalFeatures=NULL,
+        predFeatureNamesToUse=genes,
+        predFeatureRegExps=predFeatureRegExps,
+        predMat.file = predMat.file,
+        targetMat = targetMat,
+        anno.file=anno.file,
+        kMinNumSamples=kMinNumSamples,
+        save.params=F,
+        save.model=F,
+        report.summary=config.summary,
+        predictOnly=flock_common_state$predictOnly
+        )
+    }
+
+    list(targetID=targetID,
+         result=res,
+         permute.seed=permute.seed,
+         output.prefix=output.prefix,
+         group.name=group.name)
+
+  }
+
+  task.results <- lapply(flock_per_task_state, fit.single.target)
+
+  save(task.results, file=output.file)
+}
+
+    zscore <- function(x) {
+     (x-mean(x, na.rm=T))/sd(x, na.rm=T)
+    }
+
+
+gather <- function(per.task.state, common.state, output.file, run) {
+  task.results <- do.call(c, lapply(flock_per_task_state, function(job.details) { 
+    e <- new.env()
+    load(job.details$flock_output_file, envir=e)
+    e$task.results
+  }))
+
+  is.permuted <- sapply(task.results, function(x) { !is.null(x$permute.seed) } )
+
+  null.dist.details <- foreach(task.result=task.results[is.permuted], .combine=rbind) %do% {
+    load(task.result$result$fit.file)
+    if(!is.null(fit.info$OOB$weighted.cor)) {
+      data.frame(quality=fit.info$OOB$quality, weighted.cor=fit.info$OOB$weighted.cor, weighted.cor.R2=fit.info$OOB$weighted.cor.R2, with.more.weight=sum(min(fit.info$OOB$weights) != fit.info$OOB$weights), sample.count=length(fit.info$targetVec), failure.reason=NA, targetID=task.result$targetID, nfeatures=fit.info$OOB$nfeatures)
+    } else {
+      data.frame(quality=NA, weighted.cor=NA, weighted.cor.R2=NA, with.more.weight=NA, sample.count=NA, failure.reason=fit.info$failure.reason, targetID=task.result$targetID, nfeatures=NA)
+    }
+  }
+
+  null.distribution.values <- na.omit(null.dist.details$quality[!is.na(null.dist.details$quality)])
+  pval.from.null.distribution <- function(x) { (sum(null.distribution.values >= x)+1)/(length(null.distribution.values)+1) }
+
+  write.table(null.distribution.values, file=flock_common_state$quality.dist.file, col.names=F, row.names=F)
+  write.csv(null.dist.details, file=paste(flock_common_state$quality.dist.file,"-details.csv",sep=''),row.names=F)
+
+  # summarize the runs into a table 
+
+  fits <- lapply(task.results[!is.permuted], function(job) { 
+    if(!is.null(job$result) && !is.null(job$result$fit.file)) {
+        load(job$result$fit.file)
+        # reduce memory by dropping unused fields
+        top.feature <- names(which.max(fit.info$varImp))
+        if(!is.null(top.feature)) {
+          top.feature.values <- fit.info$predData[,top.feature]
+          wt.cor.top.feature <- atlantis:::weighted.cor(top.feature.values, fit.info$targetVec, fit.info$OOB$weights)
+          target.top.pred.corr <- cor(fit.info$targetVec, top.feature.values, use='pairwise.complete')
+          lines.with.more.weight <- which(min(fit.info$OOB$weights) < fit.info$OOB$weights)
+          sensitive.lines.zmean.top.feature <- mean(zscore(top.feature.values)[lines.with.more.weight], na.rm=T)
+        } else {
+          target.top.pred.corr <- NA
+          wt.cor.top.feature <- NA
+          sensitive.lines.zmean.top.feature <- NA
+        }
+        fit.info$top.feature <- top.feature
+        fit.info$target.top.pred.corr <- target.top.pred.corr
+        fit.info$wt.cor.top.feature <- wt.cor.top.feature
+        fit.info$sensitive.lines.zmean.top.feature <- sensitive.lines.zmean.top.feature
+        fit.info$predData <- NULL
+        fit.info$predictors <- NULL
+        fit.info
+    } else {
+      NULL
+    }
+  })
+
+  taskResults <- task.results[!is.permuted]
+
+  varsPerModel <- foreach(taskResult=taskResults, fit.info=fits, .combine=rbind) %do% {
+    if(length(fit.info$varImp) == 0) {
+      NULL
+    } else {
+      d <- data.frame(variable=names(fit.info$varImp), varImportance=fit.info$varImp, targetID=taskResult$targetID, group.name=taskResult$group.name)
+      rownames(d) <- NULL
+      d
+    }
+  }
+
+  coerce.to.na <- function(x) { ifelse(is.null(x), NA, x) }
+
+  qualPerModel <- foreach(taskResult=taskResults, fit.info=fits, .combine=rbind) %do% {
+    targetID <- taskResult$targetID
+    Rsquared <- coerce.to.na(fit.info$OOB$quality)
+    if(!is.null(fit.info$failure.reason)) {
+      Rsquared <- NA
+    }
+    target.top.pred.corr <- fit.info$target.top.pred.corr
+    top.feature <- fit.info$top.feature
+    wt.Rsquared <- fit.info$OOB$weighted.cor.R2
+    Pvalue <- pval.from.null.distribution(Rsquared)
+    n <- names(fit.info$varImp[top.feature])
+    group.name <- taskResult$group.name
+    if(is.null(group.name)) {
+      group.name <- NA
+    }
+    output.prefix <- taskResult$output.prefix
+    target.min <- min(fit.info$targetVec, na.rm=T)
+    sample.count <- length(fit.info$targetVec)
+    with.more.weight <- sum(min(fit.info$OOB$weights) != fit.info$OOB$weights)
+    nfeatures=fit.info$OOB$nfeatures
+    if(is.null(nfeatures)) { 
+      nfeatures <- NA
+    }
+    print("nfeatures")
+    print(nfeatures)
+
+    #stopifnot(!is.null(Rsquared))
+    #stopifnot(!is.null(targetID))
+
+    data.frame(Rsquared=Rsquared,
+               targetID=targetID, 
+               topFeature=coerce.to.na(n), 
+               Pvalue=coerce.to.na(Pvalue), 
+               topFeatureCor=target.top.pred.corr, 
+               targetMin=target.min,
+               group.name=group.name,
+               output.prefix=output.prefix,
+               wt.Rsquared=coerce.to.na(wt.Rsquared),
+               sample.count=sample.count,
+               with.more.weight=with.more.weight,
+               failure.reason=coerce.to.na(fit.info$failure.reason),
+               nfeatures=nfeatures,
+               sensitive.lines.zmean.top.feature=fit.info$sensitive.lines.zmean.top.feature,
+               wt.cor.top.feature=fit.info$wt.cor.top.feature)
+  }
+
+  qualPerModel$percentile <- rank(qualPerModel$Rsquared,na.last=F)/nrow(qualPerModel)
+  qualPerModel$fdr <- p.adjust(qualPerModel$Pvalue, method="fdr")
+
+  cat("writing summary\n")
+  write.csv(varsPerModel, file=paste(flock_run_dir,'/results/varsPerModel.csv',sep=''), row.names=F)
+  write.csv(qualPerModel, file=paste(flock_run_dir,'/results/qualPerModel.csv',sep=''), row.names=F)
+
+  annTable.file <- flock_common_state$anno.file
+
+  # make report summary plot with:
+  # 1. BRAF, KRAS, CTNNB1, ESR1
+  # 2. random 4 out of the top 20 models of that bulk run in terms of R^2
+  # 3. random 4 models with p-value > 0.1 
+
+  if(length(qualPerModel) > 4) {
+
+    sample.at.most <- function(v, n) {
+       sample(v, min(length(v), n))
+    }
+
+    randTop4 <- sample.at.most(qualPerModel$targetID[rank(-qualPerModel$Rsquared) <= 20], 4)
+    randBad4 <- sample.at.most(qualPerModel$targetID[qualPerModel$Pvalue > 0.1], 4)
+
+    add.atlantis.plots <- function(targetIds) {
+      for(targetID in targetIds) {
+        e <- new.env()
+        fi <- sprintf("%s/results/%s_fit.info_cforest.Rdata", flock_run_dir, targetID)
+        sum.file <- sprintf("%s/results/%s_ATLANTIS_Summary.Rdata", flock_run_dir, targetID)
+        if(file.exists(fi) && file.exists(sum.file)) {
+          load(fi, envir=e)
+          PlotATLANTISresults(e$fit.info, sum.file,
+            NA, NULL, annTable.file)
+        } else {
+          # write something to file to indicate target not present
+        }
+      }
+    }
+
+    pdf(file.path(flock_run_dir, 'results/summary.pdf'), w=11*1.5, h=8.5*1.5)
+    pardefault <- par()
+    grid.table(head(qualPerModel[order(qualPerModel$Rsquared, decreasing=T),],20))
+    par(pardefault)
+    add.atlantis.plots(c("GS_BRAF", "GS_KRAS", "GS_CTNNB1", "GS_ESR1"))
+    par(pardefault)
+    add.atlantis.plots(randTop4)
+    par(pardefault)
+    add.atlantis.plots(randBad4)
+    par(pardefault)
+    if(length(null.distribution.values) > 0) {
+      hist(null.distribution.values)
+    }
+    par(pardefault)
+    dev.off()
+  }
+}
+
+# given a feature ID (a compound string that includes a feature type and a feature name),
+# returns the feature name
+getFeatureNameFromFeatureID <- function(ID) {
+  res <- strsplit(ID, "[_:]")[[1]]
+  if (length(res) < 2)
+    NULL
+  else
+    res[[2]]
+}