R/IDWAS.R

In oncoPredict: Drug Response Modeling and Biomarker Discovery

#'This function maps cnv data to genes.
#'@param Cnvs The cnv data. A table with the following colnames: Sample (named using the TCGA patient barcode), Chromosome, Start, End, Num_Probes, and Segment_Mean.
#'@param folder If TRUE, write a map.RData file containing the mapped CNV matrix and tumor sample indices. The default is FALSE.
#'@keywords Map CNV data to genes
#'@return A list containing theCnvQuantVecList_mat (rows are genes, columns are samples) and tumorSamps (primary tumor/01A sample column indices).
#'@import utils
#'@import stats
#'@importFrom BiocGenerics toTable
#'@importFrom GenomicFeatures genes
#'@importFrom GenomicRanges GRanges
#'@importFrom IRanges IRanges subsetByOverlaps findOverlaps countOverlaps
#'@importFrom S4Vectors Rle
#'@examples
#' if (requireNamespace("org.Hs.eg.db", quietly=TRUE) &&
#'     requireNamespace("TxDb.Hsapiens.UCSC.hg19.knownGene", quietly=TRUE)) {
#'   cnvData <- data.frame(Sample="TCGA-01-0001-01A-01D-0000-01",
#'                         Chromosome=7,
#'                         Start=55000000,
#'                         End=55200000,
#'                         Num_Probes=10,
#'                         Segment_Mean=1.5)
#'   mapped <- map_cnv(cnvData)
#'   names(mapped)
#' }
#'@export
map_cnv<-function(Cnvs, folder=FALSE)
{
  #Check colnames() of the cnv data.
  #This is important because depending on what method you use to obtain the cnv data, colnames() might slightly differ (aka Segment.Mean and not Segment_Mean).
  if (('Segment_Mean' %in% colnames(Cnvs)) == FALSE)
    stop("\nERROR: Check colnames() of cnv data. colnames() must include Sample, Chromosome, Start, End, and Segment_Mean")

  if (!requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
    stop("Package 'org.Hs.eg.db' is required for map_cnv(). Install it with BiocManager::install('org.Hs.eg.db').",
         call. = FALSE)
  }
  if (!requireNamespace("TxDb.Hsapiens.UCSC.hg19.knownGene", quietly = TRUE)) {
    stop("Package 'TxDb.Hsapiens.UCSC.hg19.knownGene' is required for map_cnv(). Install it with BiocManager::install('TxDb.Hsapiens.UCSC.hg19.knownGene').",
         call. = FALSE)
  }

  #Load the gene ranges for HG19 using.
  txdb <- getExportedValue("TxDb.Hsapiens.UCSC.hg19.knownGene", "TxDb.Hsapiens.UCSC.hg19.knownGene")
  geneRanges <- genes(txdb) #seqnames (chr), ranges, gene id.
  e2s = toTable(getExportedValue("org.Hs.eg.db", "org.Hs.egSYMBOL")) #gene id, gene symbol.
  syms <- e2s[, "symbol"] #Gene symbols for HG19.
  names(syms) <- e2s[, "gene_id"] #Gene symbols with gene id.
  theGeneSymsOrd <- syms[as.character(geneRanges$gene_id)]

  #path<-paste(theRootDir, '/mut.txt', sep="")
  #Cnvs<-read.table(path, header=TRUE, row.names=1) #The cnv data downloaded with columns sample, chromosome, start, end, num probes, segment mean.
  CnvsList <- split(Cnvs, Cnvs[, "Sample"]) #Put each patient data into list such that each element in the list is for a patient.

  cnvEdGenesList <- list()
  ampGenesList <- list()
  delGenesList <- list()
  numGenesQuantifid <- numeric()
  theCnvQuantVecList <- list()

  for(i in 1:length(CnvsList)){ #Iterate for each patient in the list CnvsList

    chrs <- paste("chr", CnvsList[[i]]$Chromosome, sep="") #chrs contains all the chromosome numbers for these samples in patient i.
    starts <- CnvsList[[i]]$Start #starts contains all the starts.
    ends <- CnvsList[[i]]$End #ends contains all the ends.
    grCnvs <- GRanges(seqnames=Rle(chrs),ranges=IRanges(starts, ends), segMeans=CnvsList[[i]]$Segment_Mean) #chr start-end

    #Amp or del
    grCnvs_ampDel <- grCnvs[grCnvs$segMeans > 1 | grCnvs$segMeans < -1]
    cnvedGenes <- subsetByOverlaps(geneRanges, grCnvs_ampDel, type="within")
    cnvEdGenesList[[i]] <- cnvedGenes$gene_sym

    #Amps
    grCnvs_amp <- grCnvs[grCnvs$segMeans > 1]
    ampedGenes <- subsetByOverlaps(geneRanges, grCnvs_amp, type="within")
    ampGenesList[[i]] <- ampedGenes$gene_sym

    #Dels
    grCnvs_Del <- grCnvs[grCnvs$segMeans < -1]
    deledGenes <- subsetByOverlaps(geneRanges, grCnvs_Del, type="within")
    delGenesList[[i]] <- deledGenes$gene_sym

    #Continuous gene level
    #Use count overlaps to find genes that unambiguously overlap a single peak. Give it an NA it it doesn't overlap a single peak. Assign it the value of the peak if it unambiguously overlaps a peak. PC.
    numOverlaps <- countOverlaps(geneRanges, grCnvs)
    numGenesQuantifid[i] <- sum(numOverlaps == 1)
    inCnv <- which(numOverlaps == 1) # take only gene unambiguously overlaping a peak, this is usually most genes.

    theCnvQuantVec <- rep(NA, length(geneRanges))
    olaps <- findOverlaps(geneRanges, grCnvs, type="within")

    #theCnvQuantVec[olaps@queryHits] <- grCnvs$segMeans[olaps@subjectHits]

    #No slot of name 'subjectHits' for this object of class 'SortedByQueryHits'
    #This problem is caused by a small change in the GenomicAlignments package.
    theCnvQuantVec[olaps@from] <- grCnvs$segMeans[olaps@to]

    theCnvQuantVecList[[i]] <- theCnvQuantVec
    #names(theCnvQuantVecList[[i]]) <- geneRanges$gene_sym
    names(theCnvQuantVecList[[i]]) <- theGeneSymsOrd
  }

  names(theCnvQuantVecList) <- names(CnvsList)
  theCnvQuantVecList_mat <- do.call(cbind, theCnvQuantVecList)
  siteVec <- sapply(strsplit(names(CnvsList), "-"), function(l)return(l[4]))
  tumorSamps <- which(siteVec == "01A")
  output <- list(theCnvQuantVecList_mat=theCnvQuantVecList_mat, tumorSamps=tumorSamps)
  if (folder) {
    save(theCnvQuantVecList_mat, tumorSamps, file="./map.RData") # Save these RData files for use by other scripts.
    return(invisible(output))
  }

  return(output)
}
#'This function will test every drug against CNA amplifications or somatic mutations for your cancer type.
#'@param drug_prediction The drug prediction data. Must be a data frame. rownames are samples, colnames are drugs. Make sure sample names are of the same form as the sample names in your cnv or mutation data. e.g. if the rownames() are TCGA barcodes of the form TCGA-##-####-###, make sure your cnv/mutation data also uses samples in the form TCGA-##-####-###
#'@param data The CNA amplification or mutation data. Must be a data frame. If you wish to use CNA data, use the output from map_cnv(), where rows are genes and columns are samples, or use data of similar form. If you wish to use mutation data, use the method for downloading mutation data outlined in the vignette, and make sure the TCGA barcodes use '-' instead of '.'; if you use another dataset (and don't download data from TCGA), make sure your data file includes the following columns: 'Variant_Classification', 'Hugo_Symbol', 'Tumor_Sample_Barcode'.
#'@param n The minimum number of samples you want CNA amplifications or mutations to occur in. The default is 10 (arbitrarily chosen).
#'@param cnv TRUE or FALSE. Indicate whether or not you would like to test CNA amplification data. If TRUE, you will test CNA amplifications. If FALSE, you will test mutation data.
#'@param folder If TRUE, write IDWAS results to CSV files in the current working directory. The default is FALSE.
#'@keywords Test CNA amplification or mutation data to genes.
#'@import utils
#'@import stats
#'@import parallel
#'@return Raw p-value and beta-values for CNA amplifications or somatic mutations.
#'@examples
#' drugPrediction <- data.frame(drug1=c(1, 2, 3),
#'                              drug2=c(2, 1, 3),
#'                              row.names=c("S1", "S2", "S3"))
#' mutationData <- data.frame(Variant_Classification=c("Missense_Mutation",
#'                                                     "Missense_Mutation",
#'                                                     "Missense_Mutation"),
#'                            Tumor_Sample_Barcode=c("S1", "S2", "S3"),
#'                            Hugo_Symbol=c("GENE1", "GENE1", "GENE2"))
#' idwas(drugPrediction, mutationData, n=1, cnv=FALSE)
#'@export
idwas<-function(drug_prediction, data, n=10, cnv, folder=FALSE){
  #Check parameters.
  #_____________________________________________________________________________
  if (!is.data.frame(drug_prediction))
    stop("\nERROR: \"drug_prediction\" must be a data frame")
  if (!is.data.frame(data))
    stop("\nERROR: \"data\" must be a data frame")

  if (cnv == FALSE){
    drug_prediction<-t(drug_prediction)
    #This will make it such that rows are drugs and columns are samples.
  }

  #If TCGA is in my colnames() (as it would if you got cnv data from map_cnv() OR
  #if TCGA is in the column of your mutation data mutation$Tumor_Sample_Barcode, then you have TCGA samples
  #and you want to make sure you only use 01A samples).
  cols<-colnames(data)

  #_____________________________________________________________________________

  if( (sum(grepl('TCGA', cols) == TRUE) > 0) | (TRUE %in% (grepl('TCGA', data$Tumor_Sample_Barcode))) ){ #If TCGA...
    if(cnv){ #If TCGA CNV...
      #If cnv=TRUE, then you have cnv data, so proceed with this testing...
      #Find the indices of the 01A/primary tumor samples and the non-duplicates.
      indices <- which(sapply(strsplit(cols, "-"), function(a)a[4]) == "01A")
      #Create a matrix using only the primary tumor patients as columns.
      matrix <- data[, indices] #Make a matrix of rows/genes and columns/TCGA O1A samples.

      #Collect all the patient names/rows of the drug prediction matrix that are 01A primary tumors.
      drugs_01A <- rownames(drug_prediction)[which(sapply(strsplit(rownames(drug_prediction), "-"), function(a)a[4]) == "01A")]
      #Collect the patient IDs of all the patients in drugs_01A (the patients IDS for the 01A samples you have drug prediction data for). These are the 4 digits after the bacode TCGA-## so TCGA-##-????
      drugs_01Aids <- sapply(strsplit(drugs_01A, "-"), function(a)a[3])
      #Index the drug prediction matrix so that its rows/patients are the primary tumor patients.
      matrix2 <- drug_prediction[drugs_01A,, drop=FALSE]
      #Rename the rows of the drug prediction matrix to the patient id (the 4 digits).

      #Make sure you only use unique ids (sometimes there are duplicates). If there are duplicates, remove.
      indices<-match(unique(drugs_01Aids), drugs_01Aids)
      matrix2<-matrix2[indices,, drop=FALSE] #Only keep the rows that have unique ids.
      rownames(matrix2) <- drugs_01Aids[indices]

      #Collect all the patient IDs of the patients in columns of matrix
      ids <- sapply(strsplit(colnames(matrix), "-"), function(a)a[3])

      #Collect the patients that are common to both prediction and cnv/mut data (using their 4 digit ids).
      overlapping_ids<-intersect(rownames(matrix2), ids)

      #Index the drug prediction matrix so that its rows only contains patients it has in common with the cnv/mut data.
      drug_prediction2<-matrix2[overlapping_ids,, drop=FALSE]

      indices<-match(overlapping_ids, ids)
      matrix3<-as.matrix(matrix[,indices, drop=FALSE])

      MatCommonPats_amps <- apply(matrix3, 2, function(theCol)return(as.numeric(theCol > 1))) #Apply this function to each column.
      rownames(MatCommonPats_amps) <- rownames(matrix3)

      MatCommonPats_dels <- apply(matrix3, 2, function(theCol)return(as.numeric(theCol < -1)))

      rownames(MatCommonPats_dels) <- rownames(matrix3)

      theFun <- function(j){
        pVals <- numeric()
        betaVal <- numeric()
        if(sum(na.omit(MatCommonPats_amps[j, ])) > n){ #Make sure the gene is amplifed at least 50 times
          for(i in 1:ncol(drug_prediction2)){
            theMod <- coef(summary(lm(drug_prediction2[,i]~MatCommonPats_amps[j, ]))) #Linear regression betwween each amplified CNV (j/row) and each drug (i/column).
            pVals[i] <- theMod[2,4] #P-values for each model.
            betaVal[i] <- theMod[2,1] #Effect size for each model.
          }
          names(pVals) <- colnames(drug_prediction2)
          names(betaVal) <- colnames(drug_prediction2)
        }
        return(list(pVals, betaVal))
      }

      #mclapply is a parallelized vector of lapply.
      #It's structure is mclapply(X, FUN...) where X is a vector and FUN is the function applied to each element of X.
      #Here, the function is applied to each value of X which is j (each CNV).
      allCors <- mclapply(1:nrow(matrix3), theFun) #You may want to change this number based on the number of available cores.
      #allCors'is a list where each element/drug consists of 2 elements (the p value for that drug's association with the CNV and the beta value).

      #Name each element after each CNV.
      names(allCors) <- rownames(matrix3)

      hasAmps <- apply(MatCommonPats_amps, 1, function(theRow)return(sum(na.omit(theRow)) > n)) #Restrict analysis to CNAs that occur in 50 or more samples.

      allCors_hasAmps <- allCors[hasAmps]
      if(length(allCors_hasAmps) == 0){
        stop((paste("\nERROR: A gene was not amplified in at least", n, "patients. Recommend decreasing the n parameter.", sep=" ")))
      }else{
        pVals <- sapply(allCors_hasAmps, function(item)return(item[[1]]))
        betas <- sapply(allCors_hasAmps, function(item)return(item[[2]]))
        output <- list(pVals=pVals, betas=betas)
        if (folder) {
          write.csv(pVals, file='./CnvTestOutput_pVals.csv')
          write.csv(betas, file='./CnvTestOutput_betas.csv')
          return(invisible(output))
        }
        return(output)
      }
    }else{ #If TCGA MUT...
      #Obtain a list for the patients you have data for and initiate empty lists to fill.
      #_______________________________________
      tcgaIds<-data$Tumor_Sample_Barcode
      unique<-unique(tcgaIds)
      mutsListAll <- list() #A list of all the mutation occurring in each sample.
      proteinChangingMutations <- c("Missense_Mutation", "Nonsense_Mutation", "Frame_Shift_Del", "In_Frame_Del",
                                    "Frame_Shift_Ins", "In_Frame_Ins", "Nonstop_Mutation", "De_novo_Start_OutOfFrame",
                                    "De_novo_Start_InFrame", "Missense", "Read-through", "Indel")
      genesWithProteinChangeList <- list()

      #Fill those lists.
      #_______________________________________
      for(i in 1:length(unique)){
        #print(paste(i, "of", length(unique)), sep='')
        indices<-unique[i] == tcgaIds

        #Now for each sample, pull out a list of genes with somatic mutations
        variantType <- data[indices,"Variant_Classification"]
        theGenes <- data[indices,"Hugo_Symbol"]
        names(variantType) <- theGenes
        mutsListAll[[i]] <- variantType #A list of all the mutation occurring in each sample.
        genesWithProteinChangeList[[i]] <- unique(names(variantType[variantType %in% proteinChangingMutations]))
      }

      allMutatedGenes <- unique(names(unlist(mutsListAll))) #All the unique genes mutated across all patients.
      mutationTypes <- table(unlist(mutsListAll)) #Frequency of mutation type across all patients.
      mutsListAll_unlist <- unlist(mutsListAll, recursive=F) #Mutation type and gene association.
      genesWithProteinChangeList_unlist <- unlist(genesWithProteinChangeList, recursive=F) #Genes with protein changing mutations across all patients.

      #From mutsListAll we can then create a matrix indicating if the gene has a coding mutation.
      #_______________________________________
      mutMat <- numeric((length(unique)*length(allMutatedGenes)))
      dim(mutMat) <- c(length(allMutatedGenes), length(unique))
      rownames(mutMat) <- allMutatedGenes
      colnames(mutMat) <- unique

      #Now populate this matrix with the relevant information about what kind of mutation each gene has in each sample.
      #_______________________________________
      for(i in 1:length(unique)){
        #print(paste(i, "of", length(unique)), sep='')
        mutMat[genesWithProteinChangeList[[i]], i] <- rep(1, length(genesWithProteinChangeList[[i]]))
      }

      tumorTypeId <- sapply(strsplit(colnames(mutMat), "-", fixed=TRUE), function(l)return(l[4]))

      #Lets remove everything but the "Primary Solid Tumors (i.e. "01")".
      #_______________________________________
      mutMat_only01 <- mutMat[, tumorTypeId == "01A", drop=FALSE] #Genes that were mutated in 01A samples.
      theIds <- colnames(mutMat_only01) #Patient ids that were 01A.
      mutIds <- sapply(strsplit(theIds, "-", fixed=T), function(l)return(l[3])) #The TCGA #### part of the barcode.
      colnames(mutMat_only01) <- mutIds

      #Extract the 01a samples from the drug prediction data, i.e. tumor samples.
      #_______________________________________
      all01ASamples <- colnames(drug_prediction)[which(sapply(strsplit(colnames(drug_prediction), "-", fixed=T), function(a)a[4]) == "01A")]
      preds01a <- drug_prediction[, all01ASamples, drop=FALSE] #Predictions for 01A samples.
      sampIds01a <- sapply(strsplit(all01ASamples, "-", fixed=T), function(l)return(l[3])) #The TCGA #### digit number.
      colnames(preds01a) <- sampIds01a
      inPredAndMutData <- sampIds01a[sampIds01a %in% mutIds] #Samples for which we have both predicted drug response and mutation calls

      if (length(inPredAndMutData) == 0){
        stop((paste("\nERROR: Samples in drug_prediction and data don't overlap. Make sure their sample identifiers are of similar form")))
      }

      #Run the associations between all genes and drugs, for drugs with at least 50 mutations.
      #_______________________________________
      preds01a_filt_ord <- as.matrix(preds01a[, inPredAndMutData, drop=FALSE]) #The preds for the 01A samples we have both prediction and mutation data for.
      mutMat_nodups_ordFilt <- mutMat_only01[, inPredAndMutData, drop=FALSE]
      commonMuts <- apply(mutMat_nodups_ordFilt, 1, sum)
      if (length(which(commonMuts >= n)) == 0){
        stop((paste("\nERROR: A gene was not mutated in at least", n, "patients. Recommend decreasing the n parameter.", sep=" ")))
      }
      commonlyMutated <- mutMat_nodups_ordFilt[which(commonMuts >= n), , drop=FALSE] #Rows are genes that are mutated in n+ patients. Cols are patients.

      #If there are gene entries with an unknown HUGO ID, remove it.
      #_______________________________________
      if("Unknown" %in% rownames(commonlyMutated)){
        commonlyMutated<-commonlyMutated[rownames(commonlyMutated) != "Unknown",, drop=FALSE]
        if (nrow(commonlyMutated) == 0) {
          stop("\nERROR: No known genes remain after removing Unknown HUGO IDs.", call. = FALSE)
        }
      }

      #Get p values and beta values. Keeping commonlyMutated as a matrix lets the same
      #code handle one passing gene and many passing genes.
      #_______________________________________
      pValList <- list()
      betaValList <- list()

      suppressWarnings(for(i in 1:nrow(preds01a_filt_ord)){ #For each drug...
        pValList[[i]] <- numeric()
        betaValList[[i]] <- numeric()
        for(j in 1:nrow(commonlyMutated)) #For each row/gene that is mutated in n+ patients...
        {
          thecoefs <- coef(summary(lm(preds01a_filt_ord[i,]~commonlyMutated[j,])))
          pValList[[i]][[j]] <- thecoefs[2,4]
          betaValList[[i]][[j]] <- thecoefs[2,1]
        }
      })
      #Get the adjusted p-value for each gene-drug combination, pull out the significant associations
      #and create a supplementary table that lists these for "predictable" drugs?.
      #_______________________________________
      sigPs <- list()
      pAdjListCantype <- list()
      for(i in 1:length(pValList)){
        names(pValList[[i]]) <- rownames(commonlyMutated)
        names(betaValList[[i]]) <- rownames(commonlyMutated)
        padj <- p.adjust(pValList[[i]], method="BH")
        sigPs[[i]] <- padj[padj < 0.05]
        pAdjListCantype[[i]] <- padj
      }

      names(sigPs) <- rownames(preds01a_filt_ord)
      names(pValList) <- rownames(preds01a_filt_ord)
      names(betaValList) <- rownames(preds01a_filt_ord)
      names(pAdjListCantype) <- rownames(preds01a_filt_ord)

      pVal<-unlist(pValList)
      betaVal<-unlist(betaValList)

      output <- cbind(pVal, betaVal)
      if (folder) {
        write.csv(output, file='./MutationTestOutput_pVal_and_betaVal.csv')
        return(invisible(output))
      }
      return(output)
    }#The end of the first else statement.

    #This code is for when you don't have TCGA barcoded samples (it's similar to above)
  }else{ #If non TCGA...
    if(cnv){ #If non TCGA CNV...
      overlapping_samples<-intersect(rownames(drug_prediction), colnames(data))

      #Index the drug prediction matrix so that it contains the overlapping samples.
      indices<-match(overlapping_samples, rownames(drug_prediction))
      drug_prediction2<-drug_prediction[indices,]
      #Index the cnv.mutation matrix so that it contains the overlapping samples.
      indices<-match(overlapping_samples, colnames(data))
      matrix3<-as.matrix(data[,indices])

      MatCommonPats_amps <- apply(matrix3, 2, function(theCol)return(as.numeric(theCol > 1))) #Apply this function to each column.
      rownames(MatCommonPats_amps) <- rownames(matrix3)

      MatCommonPats_dels <- apply(matrix3, 2, function(theCol)return(as.numeric(theCol < -1)))

      rownames(MatCommonPats_dels) <- rownames(matrix3)

      theFun <- function(j){
        pVals <- numeric()
        betaVal <- numeric()

        if(sum(na.omit(MatCommonPats_amps[j, ])) > n){ # Make sure the gene is amplifed at least 50 times
          for(i in 1:ncol(drug_prediction2)){
            theMod <- coef(summary(lm(drug_prediction2[,i]~MatCommonPats_amps[j, ]))) #Linear regression betwween each amplified CNV (j/row) and each drug (i/column).
            pVals[i] <- theMod[2,4] #P-values for each model.
            betaVal[i] <- theMod[2,1] #Effect size for each model.
          }
          names(pVals) <- colnames(drug_prediction2)
          names(betaVal) <- colnames(drug_prediction2)
        }
        return(list(pVals, betaVal))
      }

      #mclapply is a parallelized vector of lapply.
      #It's structure is mclapply(X, FUN...) where X is a vector and FUN is the function applied to each element of X.
      #Here, the function is applied to each value of X which is j (each CNV).
      allCors <- mclapply(1:nrow(matrix3), theFun) #You may want to change this number based on the number of available cores.
      #allCors'is a list where each element/drug consists of 2 elements (the p value for that drug's association with the CNV and the beta value).

      #Name each element after each CNV.
      names(allCors) <- rownames(matrix3)

      hasAmps <- apply(MatCommonPats_amps, 1, function(theRow)return(sum(na.omit(theRow)) > n)) #Restrict analysis to CNAs that occur in 50 or more samples.

      allCors_hasAmps <- allCors[hasAmps]
      if(length(allCors_hasAmps) == 0){
        stop((paste("\nERROR: A gene was not amplified in at least", n, "patients. Recommend decreasing the n parameter.", sep=" ")))
      }else{
        pVals <- sapply(allCors_hasAmps, function(item)return(item[[1]]))
        betas <- sapply(allCors_hasAmps, function(item)return(item[[2]]))
        output <- list(pVals=pVals, betas=betas)
        if (folder) {
          write.csv(pVals, file='./CnvTestOutput_pVals.csv')
          write.csv(betas, file='./CnvTestOutput_betas.csv')
          return(invisible(output))
        }
        return(output)
      }

    }else{ #If non TCGA Mut...

      drug_prediction<-t(drug_prediction)

      #Obtain a list for the patients you have data for and initiate empty lists to fill.
      #_______________________________________
      sampIds<-data$Tumor_Sample_Barcode
      unique<-unique(sampIds)
      mutsListAll <- list() #A list of all the mutation occurring in each sample.
      proteinChangingMutations <- c("Missense_Mutation", "Nonsense_Mutation", "Frame_Shift_Del", "In_Frame_Del",
                                    "Frame_Shift_Ins", "In_Frame_Ins", "Nonstop_Mutation", "De_novo_Start_OutOfFrame",
                                    "De_novo_Start_InFrame", "Missense", "Read-through", "Indel")
      genesWithProteinChangeList <- list()

      #Fill those lists.
      #_______________________________________
      for(i in 1:length(unique)){ #For each unique sample in your mut dataset...
        #print(paste(i, "of", length(unique)), sep='')
        indices<-unique[i] == sampIds #Find the index/row of that unique sample in all the samples (in the mut data)

        #Now for each sample, pull out a list of genes with somatic mutations
        variantType <- data[indices,"Variant_Classification"] #Get its variant classification...
        theGenes <- data[indices,"Hugo_Symbol"] #Get its gene
        names(variantType) <- theGenes
        mutsListAll[[i]] <- variantType #A list of all the mutation occurring in each sample.
        genesWithProteinChangeList[[i]] <- unique(names(variantType[variantType %in% proteinChangingMutations]))
      }

      allMutatedGenes <- unique(names(unlist(mutsListAll))) #All the unique genes mutated across all patients.
      mutationTypes <- table(unlist(mutsListAll)) #Frequency of mutation type across all patients.
      mutsListAll_unlist <- unlist(mutsListAll, recursive=F) #Mutation type and gene association.
      genesWithProteinChangeList_unlist <- unlist(genesWithProteinChangeList, recursive=F) #Genes with protein changing mutations across all patients.

      #From mutsListAll we can then create a matrix indicating if the gene has a coding mutation.
      #_______________________________________
      mutMat <- numeric((length(unique)*length(allMutatedGenes)))
      dim(mutMat) <- c(length(allMutatedGenes), length(unique)) #genes/muts/rows by samples.
      rownames(mutMat) <- allMutatedGenes
      colnames(mutMat) <- unique

      #Now populate this matrix with the relevant information about what kind of mutation each gene has in each sample.
      #_______________________________________
      for(i in 1:length(unique)){
        #print(paste(i, "of", length(unique)), sep='')
        mutMat[genesWithProteinChangeList[[i]], i] <- rep(1, length(genesWithProteinChangeList[[i]]))
      }

      #Identify samples that occur in drug and mutation data.
      #_______________________________________
      drug_samps<-rownames(drug_prediction)
      inPredAndMutData <- drug_samps[drug_samps %in% unique] #Samples for which we have both predicted drug response and mutation calls

      if (length(inPredAndMutData) == 0){
        stop((paste("\nERROR: Samples in drug_prediction and data don't overlap. Make sure their sample identifiers are of similar form")))
      }

      #Run the associations between all genes and drugs, for drugs with at least 50 mutations.
      #_______________________________________
      drug_prediction_filt_ord <- as.matrix(drug_prediction[inPredAndMutData,, drop=FALSE]) #The preds for the samples we have both prediction and mutation data for.
      mutMat_nodups_ordFilt <- mutMat[, inPredAndMutData, drop=FALSE]
      commonMuts <- apply(mutMat_nodups_ordFilt, 1, sum)
      if (length(which(commonMuts >= n)) == 0){
        stop((paste("\nERROR: A gene was not mutated in at least", n, "patients. Recommend decreasing the n parameter.", sep=" ")))
      }
      #dim(mutMat_nodups_ordFilt) #genes by samples.
      commonlyMutated <- mutMat_nodups_ordFilt[which(commonMuts >= n), , drop=FALSE]

      #If there are gene entries with an unknown HUGO ID, remove it.
      #_______________________________________
      if("Unknown" %in% rownames(commonlyMutated)){
        commonlyMutated<-commonlyMutated[rownames(commonlyMutated) != "Unknown",, drop=FALSE]
        if (nrow(commonlyMutated) == 0) {
          stop("\nERROR: No known genes remain after removing Unknown HUGO IDs.", call. = FALSE)
        }
      }

      #Get p values and beta values. Keeping commonlyMutated as a matrix lets the same
      #code handle one passing gene and many passing genes.
      #_______________________________________
      pValList <- list()
      betaValList <- list()

      suppressWarnings(for(i in 1:ncol(drug_prediction_filt_ord)){ #For each drug...drug_prediction_filt_ord is samples/rows and drugs/columns. 208 11
        pValList[[i]] <- numeric()
        betaValList[[i]] <- numeric()
        for(j in 1:nrow(commonlyMutated)) #For each mut gene...commonlyMutated is mut genes/rows and samples/columns.3 208
        {
          thecoefs <- coef(summary(lm(drug_prediction_filt_ord[,i]~commonlyMutated[j,])))
          pValList[[i]][[j]] <- thecoefs[2,4]
          betaValList[[i]][[j]] <- thecoefs[2,1]
        }
      })

      #Get the adjusted p-value for each gene-drug combination, pull out the significant associations
      #and create a supplementary table that lists these for "predictable" drugs?.
      #_______________________________________
      sigPs <- list()
      pAdjListCantype <- list()
      for(i in 1:length(pValList))
      {
        names(pValList[[i]]) <- rownames(commonlyMutated)
        names(betaValList[[i]]) <- rownames(commonlyMutated)
        padj <- p.adjust(pValList[[i]], method="BH")
        sigPs[[i]] <- padj[padj < 0.05]
        pAdjListCantype[[i]] <- padj
      }

      names(sigPs) <- colnames(drug_prediction_filt_ord)
      names(pValList) <- colnames(drug_prediction_filt_ord)
      names(betaValList) <- colnames(drug_prediction_filt_ord)
      names(pAdjListCantype) <- colnames(drug_prediction_filt_ord)

      drugs=names(pValList)
      genes=paste(rownames(commonlyMutated), ":", sep="")
      mut_drug<-do.call(paste0,expand.grid(genes,drugs))
      pvalues=as.vector(unlist(pValList)) #All p values as a vector.
      betavalues=as.vector(unlist(betaValList)) #All beta values as a vector.

      df=data.frame(pvalues, betavalues)
      rownames(df)<-mut_drug

      if (folder) {
        write.csv(df, file='./MutationTestOutput_pVal_and_betaVal.csv')
        return(invisible(df))
      }
      return(df)
    }
  }
}