R/amaretto_preprocess.R

In broadinstitute/ImagingAMARETTO: ImagingAMARETTO: tools for interpreting multi-omics networks for relevance to clinical outcomes and radiographic and histopathology imaging-derived biomarkers

#' AMARETTO_Preprocess
#'
#' Wrapper code that analyzes process TCGA GISTIC (CNV) and gene expression (rna-seq or microarray) data  via one call
#' @param DataSetDirectories DataSetDirectories
#' @param BatchData BatchData
#' @importFrom graphics lines par title
#' @import grDevices
#' @importFrom limma  strsplit2
#' @importFrom MultiAssayExperiment experiments
#' @importFrom stats aov  prcomp  qqline  qqnorm  qqplot  rgamma
#' @importFrom utils data read.csv  write.table install.packages
#' @return result
#' @export
#'
#' @examples
#'\dontrun{
#' TargetDirectory <-  "Downloads"  # path to data download directory
#' CancerSite <- 'CHOL'
#' DataSetDirectories <- AMARETTO_Download(CancerSite,TargetDirectory)
#' ProcessedData <- AMARETTO_Preprocess(DataSetDirectories,BatchData)
#'}
AMARETTO_Preprocess <- function(DataSetDirectories = DataSetDirectories,
                                BatchData = BatchData) {
  CancerSite <- DataSetDirectories[1]
  MinPerBatch = 5
  MAEO <- readRDS(paste0(DataSetDirectories[2], "/", 
                         CancerSite, "_RNASeq_MAEO.rds"))
  query1 <- grep("RNASeq", names(experiments(MAEO)))
  MAEO_ge <- as.matrix(assay(MAEO[[query1]]))
  MAEO_ge <- apply(MAEO_ge, 2, log2)
  MAEO_ge[is.infinite(MAEO_ge) & MAEO_ge < 0] <- NA
  cat("Loading mRNA data.\n")
  if (!is.null(MAEO)) {
    MA_TCGA = Preprocess_MAdata(CancerSite = CancerSite, MAEO_ge = MAEO_ge,BatchData = BatchData)
  } else {
    stop("No RNASeq MAEO object found.\n")
  }
  cat("Loading CNV data.\n")
  cat("\tProcessing GISTIC output.\n")
  CGH_Data = TCGA_Load_GISTICdata(DataSetDirectories$CNVdirectory)
  CNVgenes = c(CGH_Data$AMPgenes, CGH_Data$DELgenes)
  CNVgenes = gsub("\\[", "", CNVgenes)
  CNVgenes = gsub("\\]", "", CNVgenes)
  CGH_Data$CGH_Data_Segmented = CGH_Data$CGH_Data_Segmented[unique(CNVgenes), 
                                                            ]
  SampleNames = colnames(CGH_Data$CGH_Data_Segmented)
  if (CancerSite == "LAML") {
    colnames(CGH_Data$CGH_Data_Segmented) = paste(SampleNames, 
                                                  "-03", sep = "")
  } else {
    colnames(CGH_Data$CGH_Data_Segmented) = paste(SampleNames, 
                                                  "-01", sep = "")
  }
  cat("\tBatch correction.\n")
  CNV_TCGA = TCGA_BatchCorrection_MolecularData(GEN_Data = CGH_Data$CGH_Data_Segmented,
                                                BatchData =  BatchData, MinInBatch = MinPerBatch)
  CNV_TCGA = CGH_Data$CGH_Data_Segmented
  Genes = rownames(CNV_TCGA)
  SplitGenes = strsplit2(Genes, "\\|")
  rownames(CNV_TCGA) = SplitGenes[, 1]
  CNV_TCGA = TCGA_GENERIC_MergeData(unique(rownames(CNV_TCGA)), 
                                    CNV_TCGA)
  AllCancers = c("BLCA", "BRCA", "CESC", "CHOL", 
                 "COAD", "ESCA", "GBM", "HNSC", "KIRC", "KIRP", 
                 "LAML", "LGG", "LIHC", "LUAD", "LUSC", "OV", 
                 "PAAD", "PCPG", "READ", "SARC", "STAD", "THCA", 
                 "THYM", "UCEC", "COADREAD")
  if (length(intersect(CancerSite, AllCancers)) > 
      0) {
    cat("Loading MethylMix data.\n")
    cat("\tGetting MethylMix methylation states.\n")
    eval(parse(text = paste("MET_TCGA=MethylStates$", 
                            CancerSite, sep = "")), envir = environment())
    SampleNames = colnames(MET_TCGA)
    SampleNames = gsub("\\.", "-", SampleNames)
    if (CancerSite == "LAML") {
      colnames(MET_TCGA) = paste(SampleNames, 
                                 "-03", sep = "")
    } else {
      colnames(MET_TCGA) = paste(SampleNames, 
                                 "-01", sep = "")
    }
  } else {
    cat("MethylMix data for", CancerSite, "not available\n")
    MET_TCGA = matrix(0, 1, 1)
  }
  cat("Summarizing:\n")
  cat("\tFound", length(rownames(MA_TCGA)), "genes and", 
      length(colnames(MA_TCGA)), "samples for MA data.\n")
  cat("\tFound", length(rownames(CNV_TCGA)), "genes and", 
      length(colnames(CNV_TCGA)), "samples for GISTIC data before overlapping.\n")
  cat("\tFound", length(rownames(MET_TCGA)), "genes and", 
      length(colnames(MET_TCGA)), "samples for MethylMix data before overlapping.\n")
  cat("Overlapping genes and samples for GISTIC and MethylMix.\n")
  OverlapGenes = Reduce(intersect, list(rownames(MA_TCGA), 
                                        rownames(CNV_TCGA)))
  OverlapSamples = Reduce(intersect, list(colnames(MA_TCGA), 
                                          colnames(CNV_TCGA)))
  CNV_TCGA = CNV_TCGA[OverlapGenes, OverlapSamples]
  cat("\tFound", length(OverlapGenes), "overlapping genes and", 
      length(OverlapSamples), "samples for GISTIC data.\n")
  OverlapGenes = Reduce(intersect, list(rownames(MA_TCGA), 
                                        rownames(MET_TCGA)))
  OverlapSamples = Reduce(intersect, list(colnames(MA_TCGA), 
                                          colnames(MET_TCGA)))
  MET_TCGA = MET_TCGA[OverlapGenes, OverlapSamples]
  cat("\tFound", length(OverlapGenes), "overlapping genes and", 
      length(OverlapSamples), "samples for MethylMix data.\n")
  return(list(MA_TCGA = MA_TCGA, CNV_TCGA = CNV_TCGA, 
              MET_TCGA = MET_TCGA))
}

#' TCGA_Load_GISTICdata
#' 
#' @return result
#' @keywords internal
TCGA_Load_GISTICdata <- function(GisticDirectory) {
  GenesFile = paste(GisticDirectory, "all_data_by_genes.txt", 
                    sep = "")
  CGH_Data = read.csv(GenesFile, sep = "\t", row.names = 1, 
                      header = TRUE, na.strings = "NA")
  SampleNames = colnames(CGH_Data)
  SampleNames = gsub("\\.", "-", SampleNames)
  colnames(CGH_Data) = SampleNames
  Samplegroups = TCGA_GENERIC_GetSampleGroups(colnames(CGH_Data))
  CGH_Data = TCGA_GENERIC_CleanUpSampleNames(CGH_Data, 
                                             12)
  CGH_Data = CGH_Data[, -1]
  CGH_Data = CGH_Data[, -1]
  CGH_Data = as.matrix(CGH_Data)
  GenesFile = paste(GisticDirectory, "all_thresholded.by_genes.txt", 
                    sep = "")
  CGH_Data_Thresholded = read.csv(GenesFile, sep = "\t", 
                                  row.names = 1, header = TRUE, na.strings = "NA")
  SampleNames = colnames(CGH_Data_Thresholded)
  SampleNames = gsub("\\.", "-", SampleNames)
  colnames(CGH_Data_Thresholded) = SampleNames
  CGH_Data_Thresholded = TCGA_GENERIC_CleanUpSampleNames(CGH_Data_Thresholded, 
                                                         12)
  CGH_Data_Thresholded = CGH_Data_Thresholded[, -1]
  CGH_Data_Thresholded = CGH_Data_Thresholded[, -1]
  CGH_Data_Thresholded = as.matrix(CGH_Data_Thresholded)
  AMPfile = paste(GisticDirectory, "amp_genes.conf_99.txt", 
                  sep = "")
  AMPtable = read.csv(AMPfile, sep = "\t", header = TRUE, 
                      na.strings = "NA")
  AMPgenes <- as.vector(sapply(AMPtable[4:nrow(AMPtable), 
                                        2:ncol(AMPtable)], as.character))
  AMPgenes <- unique(AMPgenes[!AMPgenes %in% c("", 
                                               NA)])
  AMPgenes <- AMPgenes[order(AMPgenes)]
  DELfile = paste(GisticDirectory, "del_genes.conf_99.txt", 
                  sep = "")
  DELtable = read.csv(DELfile, sep = "\t", header = TRUE, 
                      na.strings = "NA")
  DELgenes <- as.vector(sapply(DELtable[4:nrow(DELtable), 
                                        2:ncol(DELtable)], as.character))
  DELgenes <- unique(DELgenes[!DELgenes %in% c("", 
                                               NA)])
  DELgenes <- DELgenes[order(DELgenes)]
  cat("There are", length(AMPgenes), "AMP genes and", 
      length(DELgenes), "DEL genes.\n")
  return(list(CGH_Data_Segmented = CGH_Data, CGH_Data_Thresholded = CGH_Data_Thresholded, 
              AMPgenes = AMPgenes, DELgenes = DELgenes))
}

#' Preprocess_MAdata
#'
#' @return result
#' @keywords internal
Preprocess_MAdata <- function(CancerSite=CancerSite, MAEO_ge=MAEO_ge,BatchData=BatchData) {
  MinPerBatch = 5
  cat("\tMissing value estimation.\n")
  MA_TCGA = TCGA_Load_MolecularData(MAEO_ge)
  Samplegroups = TCGA_GENERIC_GetSampleGroups(colnames(MA_TCGA))
  if (CancerSite == "LAML") {
    MA_TCGA = MA_TCGA[, Samplegroups$PeripheralBloodCancer]
  } else {
    MA_TCGA = MA_TCGA[, Samplegroups$Primary]
  }
  cat("\tBatch correction.\n")
  MA_TCGA = TCGA_BatchCorrection_MolecularData(MA_TCGA, 
                                               BatchData, MinPerBatch)
  cat("\tProcessing gene ids and merging.\n")
  Genes = rownames(MA_TCGA)
  SplitGenes = strsplit2(Genes, "\\|")
  rownames(MA_TCGA) = SplitGenes[, 1]
  MA_TCGA = MA_TCGA[!rownames(MA_TCGA) %in% "?", 
                    ]
  MA_TCGA = TCGA_GENERIC_MergeData(unique(rownames(MA_TCGA)), 
                                   MA_TCGA)
  return(MA_TCGA = MA_TCGA)
}


#' TCGA_Load_MolecularData
#'
#' @return result
#' @importFrom impute impute.knn
#' @keywords internal
TCGA_Load_MolecularData <- function(MAEO_ge) {
  MET_Data <- MAEO_ge
  if (rownames(MET_Data)[1] == "Composite Element REF") {
    cat("Removing first row with text stuff.\n")
    MET_Data = MET_Data[-1, ]
    Genes = rownames(MET_Data)
    MET_Data = apply(MET_Data, 2, as.numeric)
    rownames(MET_Data) = Genes
  }
  SampleNames = colnames(MET_Data)
  SampleNames = gsub("\\.", "-", SampleNames)
  colnames(MET_Data) = SampleNames
  MissingValueThreshold = 0.1
  NrMissingsPerGene = apply(MET_Data, 1, function(x) sum(is.na(x)))/ncol(MET_Data)
  cat("Removing", sum(NrMissingsPerGene > MissingValueThreshold), 
      "genes with more than 10% missing values.\n")
  if (sum(NrMissingsPerGene > MissingValueThreshold) > 
      0) 
    MET_Data = MET_Data[NrMissingsPerGene < MissingValueThreshold, 
                        ]
  NrMissingsPerSample = apply(MET_Data, 2, function(x) sum(is.na(x)))/nrow(MET_Data)
  cat("Removing", sum(NrMissingsPerSample > MissingValueThreshold), 
      "patients with more than 10% missing values.\n")
  if (sum(NrMissingsPerSample > MissingValueThreshold) > 
      0) 
    MET_Data = MET_Data[, NrMissingsPerSample < 
                          MissingValueThreshold]
  if (length(colnames(MET_Data)) > 1) {
    k = 15
    KNNresults = impute.knn(as.matrix(MET_Data), 
                            k)
    MET_Data_KNN = KNNresults$data
    MET_Data_KNN_Clean = TCGA_GENERIC_CleanUpSampleNames(MET_Data_KNN, 
                                                         15)
    return(MET_Data_KNN_Clean)
  } else {
    MET_Data_Clean = TCGA_GENERIC_CleanUpSampleNames(MET_Data, 
                                                     15)
    return(MET_Data_Clean)
  }
}


#' TCGA_GENERIC_CleanUpSampleNames
#' 
#' @return result
#' @keywords internal
TCGA_GENERIC_CleanUpSampleNames <- function(GEN_Data=GEN_Data, 
                                            IDlength = 12) {
  SampleNames = colnames(GEN_Data)
  SampleNamesShort = as.character(apply(as.matrix(SampleNames), 
                                        2, substr, 1, IDlength))
  if (length(SampleNamesShort) != length(unique(SampleNamesShort))) {
    Counts = table(SampleNamesShort)
    Doubles = rownames(Counts)[which(Counts > 1)]
    cat("Removing doubles for", length(Doubles), 
        "samples.\n")
    for (i in 1:length(Doubles)) {
      pos = grep(Doubles[i], SampleNames)
      GEN_Data = GEN_Data[, -pos[2:length(pos)]]
      SampleNames = colnames(GEN_Data)
    }
    SampleNames = colnames(GEN_Data)
    SampleNamesShort = as.character(apply(as.matrix(SampleNames), 
                                          2, substr, 1, IDlength))
    colnames(GEN_Data) = SampleNamesShort
  } else {
    colnames(GEN_Data) = SampleNamesShort
  }
  return(GEN_Data)
}

#' TCGA_GENERIC_GetSampleGroups
#'
#' @return result
#' @keywords internal
TCGA_GENERIC_GetSampleGroups <- function(SampleNames) {
  SampleGroups = list()
  Matches = regexpr("TCGA[.|-]\\w\\w[.|-]\\w\\w\\w\\w[.|-]01[.|-]*", 
                    SampleNames, perl = FALSE, useBytes = FALSE)
  SampleGroups$Primary = SampleNames[Matches == 1]
  Matches = regexpr("TCGA[.|-]\\w\\w[.|-]\\w\\w\\w\\w[.|-]02[.|-]*", 
                    SampleNames, perl = FALSE, useBytes = FALSE)
  SampleGroups$Recurrent = SampleNames[Matches == 
                                         1]
  Matches = regexpr("TCGA[.|-]\\w\\w[.|-]\\w\\w\\w\\w[.|-]03[.|-]*", 
                    SampleNames, perl = FALSE, useBytes = FALSE)
  SampleGroups$PeripheralBloodCancer = SampleNames[Matches == 
                                                     1]
  Matches = regexpr("TCGA[.|-]\\w\\w[.|-]\\w\\w\\w\\w[.|-]10[.|-]*", 
                    SampleNames, perl = FALSE, useBytes = FALSE)
  SampleGroups$BloodNormal = SampleNames[Matches == 
                                           1]
  Matches = regexpr("TCGA[.|-]\\w\\w[.|-]\\w\\w\\w\\w[.|-]11[.|-]*", 
                    SampleNames, perl = FALSE, useBytes = FALSE)
  SampleGroups$SolidNormal = SampleNames[Matches == 
                                           1]
  Matches = regexpr("TCGA[.|-]\\w\\w[.|-]\\w\\w\\w\\w[.|-]20[.|-]*", 
                    SampleNames, perl = FALSE, useBytes = FALSE)
  SampleGroups$CellLines = SampleNames[Matches == 
                                         1]
  Matches = regexpr("TCGA[.|-]\\w\\w[.|-]\\w\\w\\w\\w[.|-]06[.|-]*", 
                    SampleNames, perl = FALSE, useBytes = FALSE)
  SampleGroups$Metastatic = SampleNames[Matches == 
                                          1]
  return(SampleGroups)
}

#' TCGA_BatchCorrection_MolecularData
#' 
#' @return result
#' @keywords internal
TCGA_BatchCorrection_MolecularData <- function(GEN_Data =GEN_Data,
                                               BatchData = BatchData,
                                               MinInBatch = MinInBatch) {
  if (length(-which(BatchData[, 3] == 0)) > 0) {
    BatchData = BatchData[-which(BatchData[, 3] == 
                                   0), ]
  }
  PresentSamples = is.element(BatchData[, 1], colnames(GEN_Data))
  BatchDataSelected = BatchData[PresentSamples, ]
  if (sum(PresentSamples) != length(colnames(GEN_Data))) 
    BatchDataSelected = BatchData[-which(PresentSamples == 
                                           FALSE), ]
  BatchDataSelected$Batch <- factor(BatchDataSelected$Batch)
  NrPerBatch = table(BatchDataSelected$Batch)
  SmallBatches = NrPerBatch < MinInBatch
  BatchesToBeRemoved = names(SmallBatches)[which(SmallBatches == 
                                                   TRUE)]
  SamplesToBeRemoved = as.character(BatchDataSelected[which(BatchDataSelected$Batch %in% 
                                                              BatchesToBeRemoved), 1])
  if (length(colnames(GEN_Data)) - length(which(colnames(GEN_Data) %in% 
                                                SamplesToBeRemoved)) > 5) {
    # just checking if we have enough samples after
    # removing the too small batches
    if (length(which(colnames(GEN_Data) %in% SamplesToBeRemoved)) > 
        0) {
      cat("Removing", length(which(colnames(GEN_Data) %in% 
                                     SamplesToBeRemoved)), "samples because their batches are too small.\n")
      GEN_Data = GEN_Data[, -which(colnames(GEN_Data) %in% 
                                     SamplesToBeRemoved)]
    }
    BatchCheck = TCGA_GENERIC_CheckBatchEffect(GEN_Data, 
                                               BatchData)
    if (is.list(BatchCheck)) {
      GEN_Data_Corrected = TCGA_GENERIC_BatchCorrection(GEN_Data, 
                                                        BatchData)
      BatchCheck = TCGA_GENERIC_CheckBatchEffect(GEN_Data_Corrected, 
                                                 BatchData)
      return(GEN_Data_Corrected)
    } else {
      cat("Only one batch, no batch correction possible.\n")
      return(GEN_Data)
    }
  } else {
    cat("The nr of samples becomes to small, no batch correction possible.\n")
    return(GEN_Data)
  }
}

#' TCGA_GENERIC_BatchCorrection
#'
#' @return result
#' @keywords internal
TCGA_GENERIC_BatchCorrection <- function(GEN_Data = GEN_Data, 
                                         BatchData = BatchData) {
  WithBatchSamples = is.element(colnames(GEN_Data), 
                                BatchData[, 1])
  if (length(which(WithBatchSamples == FALSE)) > 
      0) 
    GEN_Data = GEN_Data[, -which(WithBatchSamples == 
                                   FALSE)]
  PresentSamples = is.element(BatchData[, 1], colnames(GEN_Data))
  BatchDataSelected = BatchData
  if (sum(PresentSamples) != length(colnames(GEN_Data))) 
    BatchDataSelected = BatchData[-which(PresentSamples == 
                                           FALSE), ]
  BatchDataSelected$Batch <- factor(BatchDataSelected$Batch)
  BatchDataSelected$ArrayName <- factor(BatchDataSelected$ArrayName)
  order <- match(colnames(GEN_Data), BatchDataSelected[, 
                                                       1])
  BatchDataSelected = BatchDataSelected[order, ]
  BatchDataSelected$Batch <- factor(BatchDataSelected$Batch)
  CombatResults = ComBat_NoFiles(GEN_Data, BatchDataSelected)
  GEN_Data_Corrected = CombatResults[, -1]
  class(GEN_Data_Corrected) <- "numeric"
  return(GEN_Data_Corrected)
}

#' ComBat_NoFiles
#' 
#' @return result
#' @keywords internal
ComBat_NoFiles <- function(dat, saminfo, type = "txt", 
                           write = FALSE, covariates = "all", par.prior = TRUE, 
                           filter = FALSE, skip = 0, prior.plots = FALSE) {
  cat("Reading Sample Information File\n")
  if (sum(colnames(saminfo) == "Batch") != 1) {
    return("ERROR: Sample Information File does not have a Batch column!")
  }
  expression_xls = "exp.txt"
  cat("Reading Expression Data File\n")
  dat <- dat[, trim.dat(dat)]
  if (skip > 0) {
    geneinfo <- as.matrix(dat[, 1:skip])
    dat <- dat[, -c(1:skip)]
  } else {
    geneinfo = NULL
  }
  if (filter) {
    ngenes <- nrow(dat)
    col <- ncol(dat)/2
    present <- apply(dat, 1, filter.absent, filter)
    dat <- dat[present, -(2 * (1:col))]
    if (skip > 0) {
      geneinfo <- geneinfo[present, ]
    }
    cat("Filtered genes absent in more than", filter, 
        "of samples. Genes remaining:", nrow(dat), 
        "; Genes filtered:", ngenes - nrow(dat), 
        "\n")
  }
  if (any(apply(dat, 2, mode) != "numeric")) {
    return("ERROR: Array expression columns contain non-numeric values! (Check your .xls file for non-numeric values and if this is not the problem, make a .csv file and use the type=csv option)")
  }
  tmp <- match(colnames(dat), saminfo[, 1])
  if (any(is.na(tmp))) {
    return("ERROR: Sample Information File and Data Array Names are not the same!")
  }
  tmp1 <- match(saminfo[, 1], colnames(dat))
  saminfo <- saminfo[tmp1[!is.na(tmp1)], ]
  if (any(covariates != "all")) {
    saminfo <- saminfo[, c(1:2, covariates)]
  }
  design <- design.mat(saminfo)
  batches <- list.batch(saminfo)
  n.batch <- length(batches)
  n.batches <- sapply(batches, length)
  n.array <- sum(n.batches)
  NAs = any(is.na(dat))
  if (NAs) {
    cat(c("Found", sum(is.na(dat)), "Missing Data Values\n"), 
        sep = " ")
  }
  cat("Standardizing Data across genes\n")
  if (!NAs) {
    B.hat <- solve(t(design) %*% design) %*% t(design) %*% 
      t(as.matrix(dat))
  } else {
    B.hat = apply(dat, 1, Beta.NA, design)
  }
  grand.mean <- t(n.batches/n.array) %*% B.hat[1:n.batch, 
                                               ]
  if (!NAs) {
    var.pooled <- ((dat - t(design %*% B.hat))^2) %*% 
      rep(1/n.array, n.array)
  } else {
    var.pooled <- apply(dat - t(design %*% B.hat), 
                        1, var, na.rm = TRUE)
  }
  stand.mean <- t(grand.mean) %*% t(rep(1, n.array))
  if (!is.null(design)) {
    tmp <- design
    tmp[, c(1:n.batch)] <- 0
    stand.mean <- stand.mean + t(tmp %*% B.hat)
  }
  s.data <- (dat - stand.mean)/(sqrt(var.pooled) %*% 
                                  t(rep(1, n.array)))
  cat("Fitting L/S model and finding priors\n")
  batch.design <- design[, 1:n.batch]
  if (!NAs) {
    gamma.hat <- solve(t(batch.design) %*% batch.design) %*% 
      t(batch.design) %*% t(as.matrix(s.data))
  } else {
    gamma.hat = apply(s.data, 1, Beta.NA, batch.design)
  }
  delta.hat <- NULL
  for (i in batches) {
    delta.hat <- rbind(delta.hat, apply(s.data[, 
                                               i], 1, var, na.rm = TRUE))
  }
  gamma.bar <- apply(gamma.hat, 1, mean)
  t2 <- apply(gamma.hat, 1, var)
  a.prior <- apply(delta.hat, 1, aprior)
  b.prior <- apply(delta.hat, 1, bprior)
  if (prior.plots & par.prior) {
    pdf(file = "prior_plots.pdf")
    par(mfrow = c(2, 2))
    tmp <- density(gamma.hat[1, ])
    plot(tmp, type = "l", main = "Density Plot")
    xx <- seq(min(tmp$x), max(tmp$x), length = 100)
    lines(xx, dnorm(xx, gamma.bar[1], sqrt(t2[1])), 
          col = 2)
    qqnorm(gamma.hat[1, ])
    qqline(gamma.hat[1, ], col = 2)
    tmp <- density(delta.hat[1, ])
    invgam <- 1/rgamma(ncol(delta.hat), a.prior[1], 
                       b.prior[1])
    tmp1 <- density(invgam)
    plot(tmp, typ = "l", main = "Density Plot", 
         ylim = c(0, max(tmp$y, tmp1$y)))
    lines(tmp1, col = 2)
    qqplot(delta.hat[1, ], invgam, xlab = "Sample Quantiles", 
           ylab = "Theoretical Quantiles")
    lines(c(0, max(invgam)), c(0, max(invgam)), 
          col = 2)
    title("Q-Q Plot")
    dev.off()
  }
  gamma.star <- delta.star <- NULL
  if (par.prior) {
    cat("Finding parametric adjustments\n")
    for (i in 1:n.batch) {
      temp <- it.sol(s.data[, batches[[i]]], 
                     gamma.hat[i, ], delta.hat[i, ], gamma.bar[i], 
                     t2[i], a.prior[i], b.prior[i])
      gamma.star <- rbind(gamma.star, temp[1, 
                                           ])
      delta.star <- rbind(delta.star, temp[2, 
                                           ])
    }
  } else {
    cat("Finding nonparametric adjustments\n")
    for (i in 1:n.batch) {
      temp <- int.eprior(as.matrix(s.data[, batches[[i]]]), 
                         gamma.hat[i, ], delta.hat[i, ])
      gamma.star <- rbind(gamma.star, temp[1, 
                                           ])
      delta.star <- rbind(delta.star, temp[2, 
                                           ])
    }
  }
  cat("Adjusting the Data\n")
  bayesdata <- s.data
  j <- 1
  for (i in batches) {
    bayesdata[, i] <- (bayesdata[, i] - t(batch.design[i, 
                                                       ] %*% gamma.star))/(sqrt(delta.star[j, 
                                                                                           ]) %*% t(rep(1, n.batches[j])))
    j <- j + 1
  }
  bayesdata <- (bayesdata * (sqrt(var.pooled) %*% 
                               t(rep(1, n.array)))) + stand.mean
  if (write) {
    output_file <- paste(expression_xls, "Adjusted", 
                         ".txt", sep = "_")
    outdata <- cbind(ProbeID = rownames(dat), bayesdata)
    write.table(outdata, file = output_file, sep = "\t")
    cat("Adjusted data saved in file:", output_file, 
        "\n")
  } else {
    return(cbind(rownames(dat), bayesdata))
  }
}

#' filter.absent
#'
#' filters data based on presence/absence call
#' @return result
#' @keywords internal
filter.absent <- function(x, pct) {
  present <- TRUE
  col <- length(x)/2
  pct.absent <- (sum(x[2 * (1:col)] == "A") + sum(x[2 * 
                                                      (1:col)] == "M"))/col
  if (pct.absent > pct) {
    present <- FALSE
  }
  return(present)
}

#' build.design
#'
#' Next two functions make the design matrix (X) from the sample info file
#' @return result
#' @keywords internal
build.design <- function(vec, des = NULL, start = 2) {
  tmp <- matrix(0, length(vec), nlevels(vec) - start + 
                  1)
  for (i in 1:ncol(tmp)) {
    tmp[, i] <- vec == levels(vec)[i + start - 
                                     1]
  }
  return(cbind(des, tmp))
}

#' design.mat
#'
#' @return result
#' @keywords internal
design.mat <- function(saminfo) {
  tmp <- which(colnames(saminfo) == "Batch")
  tmp1 <- as.factor(saminfo[, tmp])
  cat("Found", nlevels(tmp1), "batches\n")
  design <- build.design(tmp1, start = 1)
  ncov <- ncol(as.matrix(saminfo[, -c(1:2, tmp)]))
  cat("Found", ncov, "covariate(s)\n")
  if (ncov > 0) {
    for (j in 1:ncov) {
      tmp1 <- as.factor(as.matrix(saminfo[, -c(1:2, 
                                               tmp)])[, j])
      design <- build.design(tmp1, des = design)
    }
  }
  return(design)
}

#' list.batch
#'
#' Makes a list with elements pointing to which array belongs to which batch
#' @return result
#' @keywords internal
list.batch <- function(saminfo) {
  tmp1 <- as.factor(saminfo[, which(colnames(saminfo) == 
                                      "Batch")])
  batches <- lapply(1:nlevels(tmp1), function(x) (1:length(tmp1))[tmp1 == 
                                                                    levels(tmp1)[x]])
  return(batches)
}

#' trim.dat
#'
#' Trims the data of extra columns, note your array names cannot be named 'X' or start with 'X.'
#' @return result
#' @keywords internal
trim.dat <- function(dat) {
  tmp <- strsplit(colnames(dat), "\\.")
  tr <- sapply(1:length(tmp), function(x) tmp[[x]][1] != 
                 "X")
  return(tr)
}

#' aprior
#'
#' Following four find empirical hyper-prior values
#' @return result
#' @keywords internal
aprior <- function(gamma.hat) {
  m = mean(gamma.hat)
  s2 = var(gamma.hat)
  return((2 * s2 + m^2)/s2)
}

#' bprior
#'
#' @return result
#' @keywords internal
bprior <- function(gamma.hat) {
  m = mean(gamma.hat)
  s2 = var(gamma.hat)
  return((m * s2 + m^3)/s2)
}

#' postmean
#'
#' @return result
#' @keywords internal
postmean <- function(g.hat, g.bar, n, d.star, t2) {
  return((t2 * n * g.hat + d.star * g.bar)/(t2 * 
                                              n + d.star))
}

#' postvar
#'
#' @return result
#' @keywords internal
postvar <- function(sum2, n, a, b) {
  return((0.5 * sum2 + b)/(n/2 + a - 1))
}

#' it.sol
#'
#' Pass in entire data set, the design matrix for the entire data, the batch means, the batch variances, priors (m, t2, a, b), columns of the data  matrix for the batch. Uses the EM to find the parametric batch adjustments
#' 
#' @return result
#' @keywords internal
it.sol <- function(sdat, g.hat, d.hat, g.bar, t2, a, 
                   b, conv = 1e-04) {
  n <- apply(!is.na(sdat), 1, sum)
  g.old <- g.hat
  d.old <- d.hat
  change <- 1
  count <- 0
  while (change > conv) {
    g.new <- postmean(g.hat, g.bar, n, d.old, t2)
    sum2 <- apply((sdat - g.new %*% t(rep(1, ncol(sdat))))^2, 
                  1, sum, na.rm = TRUE)
    d.new <- postvar(sum2, n, a, b)
    change <- max(abs(g.new - g.old)/g.old, abs(d.new - 
                                                  d.old)/d.old)
    g.old <- g.new
    d.old <- d.new
    count <- count + 1
  }
  adjust <- rbind(g.new, d.new)
  rownames(adjust) <- c("g.star", "d.star")
  return(adjust)
}

#' L
#'
#' likelihood function
#' 
#' @return result
#' @keywords internal
L <- function(x, g.hat, d.hat) {
  return(prod(dnorm(x, g.hat, sqrt(d.hat))))
}

#' int.eprior
#'
#'Monte Carlo integration function to find the nonparametric adjustments
#' @return result
#' @keywords internal
int.eprior <- function(sdat, g.hat, d.hat) {
  g.star <- d.star <- NULL
  r <- nrow(sdat)
  for (i in 1:r) {
    g <- g.hat[-i]
    d <- d.hat[-i]
    x <- sdat[i, !is.na(sdat[i, ])]
    n <- length(x)
    j <- numeric(n) + 1
    dat <- matrix(as.numeric(x), length(g), n, 
                  byrow = TRUE)
    resid2 <- (dat - g)^2
    sum2 <- resid2 %*% j
    LH <- 1/(2 * pi * d)^(n/2) * exp(-sum2/(2 * 
                                              d))
    LH[LH == "NaN"] = 0
    g.star <- c(g.star, sum(g * LH)/sum(LH))
    d.star <- c(d.star, sum(d * LH)/sum(LH))
  }
  adjust <- rbind(g.star, d.star)
  rownames(adjust) <- c("g.star", "d.star")
  return(adjust)
}


#' Beta.NA
#'
#' @return result
#' @keywords internal
Beta.NA = function(y, X) {
  des = X[!is.na(y), ]
  y1 = y[!is.na(y)]
  B <- solve(t(des) %*% des) %*% t(des) %*% y1
  return(B)
}

#' TCGA_GENERIC_MergeData
#'
#' @return result
#' @keywords internal
TCGA_GENERIC_MergeData <- function(NewIDListUnique, 
                                   DataMatrix, MergeMethod) {
  NrUniqueGenes = length(NewIDListUnique)
  MergedData = matrix(0, NrUniqueGenes, length(colnames(DataMatrix)))
  for (i in 1:NrUniqueGenes) {
    tmpData = DataMatrix[which(rownames(DataMatrix) %in% 
                                 NewIDListUnique[i]), ]
    ifelse(length(rownames(tmpData)) > 1, MergedData[i, 
                                                     ] <- colMeans(tmpData), MergedData[i, ] <- tmpData)
  }
  rownames(MergedData) = NewIDListUnique
  colnames(MergedData) = colnames(DataMatrix)
  return(MergedData)
}

#' TCGA_GENERIC_CheckBatchEffect
#'
#' @return result
#' @keywords internal
TCGA_GENERIC_CheckBatchEffect <- function(GEN_Data, 
                                          BatchData) {
  Order = match(colnames(GEN_Data), BatchData[, 1])
  BatchDataSelected = BatchData[Order, ]
  BatchDataSelected$Batch <- factor(BatchDataSelected$Batch)
  PCAanalysis = prcomp(t(GEN_Data))
  PCdata = PCAanalysis$x
  #plot(PCdata[, 1] ~ BatchDataSelected[, 3])
  if (length(unique(BatchDataSelected$Batch)) > 1) {
    tmp = aov(PCdata[, 1] ~ BatchDataSelected[, 
                                              3])
    return(list(Pvalues = summary(tmp), PCA = PCdata, 
                BatchDataSelected = BatchDataSelected))
  } else {
    return(-1)
  }
}

#' Save_CancerSite
#'
#' @return result
#' @keywords internal
Save_CancerSite <- function(CancerSite, TargetDirectory, 
                            DataSetDirectories, ProcessedData) {
  if (length(DataSetDirectories$MAdirectory) > 1) {
    tmpMAdir = DataSetDirectories$MAdirectory[1]
  } else {
    tmpMAdir = DataSetDirectories$MAdirectory
  }
  MAdir = strsplit2(tmpMAdir, "/")
  tmpPos = grep("gdac", MAdir)
  MAversion = gsub("gdac_", "", MAdir[tmpPos[1]])
  CNVdir = strsplit2(DataSetDirectories$CNVdirectory, 
                     "/")
  CNVversion = gsub("gdac_", "", CNVdir[tmpPos[1]])
  METversion = "MethylMix2015"
  SaveFile = paste(TargetDirectory, "TCGA_", CancerSite, 
                   "_ProcessedData_MA", MAversion, "_CNV", CNVversion, 
                   "_MET", METversion, ".RData", sep = "")
  save(file = SaveFile, ProcessedData)
  return(SaveFile)
}