Mercurial > repos > morinlab > maftools
view oncodriveclust_merge.R @ 2:c3ed203d814d draft
Uploaded
| author | morinlab |
|---|---|
| date | Sun, 04 Dec 2016 20:47:23 -0500 |
| parents | |
| children |
line wrap: on
line source
require(maftools); library(argparse); require(data.table); ### parser <- ArgumentParser(description="Create a Gene Lollipop using Maftools"); parser$add_argument( "--input_maf", "-maf", required="True", help="Input Variants in MAF format" ); parser$add_argument( "--gene_blacklist", "-gl", help="Input gene list with separated by newline" ); parser$add_argument( "--min_mut", "-mm", default=5, help="Minimum number of mutations seen in the gene for it to be included in the calculation"); parser$add_argument( "--fdr", "-f", default=0.1, help="FDR threshold to use in plots and returned gene list"); parser$add_argument( "--aacol", "-ac", help="Optionally provide the name of the column that contains the amino acid annotation in your MAF file"); parser$add_argument( "--output_detail", "-o", required="True", help="Output text file for oncodriveclust detail" ) parser$add_argument( "--output_plot", "-p", required="True", help="Output pdf file for oncodriveclust detail" ) args <- parser$parse_args(); ### aacol = 'HGVSp_Short' if(!is.null(args$aacol)){ aacol = args$aacol } min_mut = as.integer(args$min_mut) #--------------------- based on binaomial distribution, estimate threshhold. get_threshold = function(gene_muts, gene_length){ th = which(unlist(lapply(X = 2:gene_muts, FUN = function(x) dbinom(x = x, size = gene_muts, prob = 1/gene_length) )) < 0.01)[1] return(th+1) } #-------------------- end of function. parse_prot_fix = function(dat, AACol, gl, m, calBg = FALSE, nBg){ if(is.null(AACol)){ if(! 'AAChange' %in% colnames(dat)){ message('Available fields:') print(colnames(dat)) stop('AAChange field not found in MAF. Use argument AACol to manually specifiy field name containing protein changes.') } }else{ colnames(dat)[which(colnames(dat) == AACol)] = 'AAChange' } all.prot.dat = dat[,.(Hugo_Symbol, Variant_Classification, AAChange)] all.prot.dat = all.prot.dat[Variant_Classification != 'Splice_Site'] #parse AAchanges to get postion prot.spl = strsplit(x = as.character(all.prot.dat$AAChange), split = '.', fixed = TRUE) prot.conv = sapply(prot.spl, function(x) x[length(x)]) all.prot.dat[,conv := prot.conv] all.prot.dat = all.prot.dat[!conv == 'NULL'] #If conversions are in HGVSp_long (default HGVSp) format, we will remove strings Ter followed by anything (e.g; p.Asn1986GlnfsTer13) pos = gsub(pattern = 'Ter.*', replacement = '',x = all.prot.dat$conv) #Following parsing takes care of most of HGVSp_short and HGVSp_long format pos = gsub(pattern = '[[:alpha:]]', replacement = '', x = pos) pos = gsub(pattern = '\\*$', replacement = '', x = pos) #Remove * if nonsense mutation ends with * pos = gsub(pattern = '^\\*', replacement = '', x = pos) #Remove * if nonsense mutation starts with * pos = gsub(pattern = '\\*.*', replacement = '', x = pos) #Remove * followed by position e.g, p.C229Lfs*18 pos = suppressWarnings( as.numeric(sapply(strsplit(x = pos, split = '_', fixed = TRUE), '[', 1)) ) all.prot.dat[,pos := pos] if(nrow( all.prot.dat[is.na(all.prot.dat$pos),]) > 0){ #message(paste('Removed', nrow( all.prot.dat[is.na(all.prot.dat$pos),]), 'mutations for which AA position was not available', sep = ' ')) #print(prot.dat[is.na(prot.dat$pos),]) all.prot.dat = all.prot.dat[!is.na(all.prot.dat$pos),] } gene.sum = summarizeMaf_fix(maf = dat)$gene.summary #gene.sum = merge.data.frame(x = gene.sum, y = gl, by = 'Hugo_Symbol', all.x = TRUE) gene.sum = merge(x = gene.sum, y = gl, by = 'Hugo_Symbol', all.x = TRUE) #gene.sum = gene.sum[!is.na(gene.sum$aa.length),] gene.sum = gene.sum[!is.na(gene.sum$aa.length)] num_mut_colIndex = which(colnames(gene.sum) == 'total') aalen_colIndex = which(colnames(gene.sum) == 'aa.length') #Get background threshold gene.sum$th = apply(gene.sum, 1, function(x) get_threshold(gene_muts = as.numeric(x[num_mut_colIndex]), gene_length = as.numeric(x[aalen_colIndex]))) #use only genes with atleast 2 (or m ) mutations. gene.sum = gene.sum[total >= m] if(calBg){ if(nrow(gene.sum) < nBg){ #message("Not enough genes to build background. Using predefined values. (Mean = 0.279; SD = 0.13)") return(NULL) } else{ syn.res = c() pb <- txtProgressBar(min = 0, max = nrow(gene.sum), style = 3) #progress bar for(i in 1:nrow(gene.sum)){ prot.dat = all.prot.dat[Hugo_Symbol == gene.sum[i, "Hugo_Symbol"]] syn.res = rbind(syn.res, cluster_prot_fix(prot.dat = prot.dat, gene = gene.sum[i, "Hugo_Symbol"], th = gene.sum[i,"th"], protLen = gene.sum[i,"aa.length"])) setTxtProgressBar(pb, i) } return(syn.res) } } else{ nonsyn.res = c() pb <- txtProgressBar(min = 0, max = nrow(gene.sum), style = 3) #progress bar for(i in 1:nrow(gene.sum)){ hs = gene.sum[i, Hugo_Symbol] #print(hs) prot.dat = all.prot.dat[Hugo_Symbol %in% hs] nonsyn.res = rbind(nonsyn.res, cluster_prot_fix(prot.dat = prot.dat, gene = hs, th = gene.sum[Hugo_Symbol %in% hs, th], protLen = gene.sum[Hugo_Symbol %in% hs, aa.length])) setTxtProgressBar(pb, i) } return(nonsyn.res) } } cluster_prot_fix = function(prot.dat, gene, th, protLen){ mergeDist = 5 #hard coded inter event distance. #prot.dat = all.prot.dat[Hugo_Symbol == gene] #Summarise counts per position pos.counts = prot.dat[,.N,pos] pos.counts = pos.counts[order(pos)] #classify position as meaningful if its greater than background threshhold. pos.counts$cluster = ifelse(test = pos.counts$N >= th, yes = 'meaningful', no = 'nonMeaningful') #Just choose meaningful positions clust.tbl = pos.counts[cluster %in% 'meaningful'] nonclust.tbl = pos.counts[cluster %in% 'nonMeaningful'] if(nrow(clust.tbl) == 0){ #message(paste('No meaningful positions found for', gene, sep=' ')) return(NULL) } clust.tbl$distance = c(0,diff(clust.tbl$pos)) #calculate inter event distance. #If more than one meaningful positions are found within a 5 aa distance, join them to form a cluster. if(nrow(clust.tbl) > 1){ #initialize variables. cstart = end = clust.tbl[1,pos] n = clust.tbl[1,N] cdf = c() cluster = 1 #Go through entire table and update variables. for(i in 2:nrow(clust.tbl)){ pos = clust.tbl[i,pos] d = clust.tbl[i,distance] if(d < mergeDist){ end = pos n = n + clust.tbl[i,N] }else{ tempdf = data.frame(cluster = paste('cluster', cluster, sep='_'), start = cstart, end = end ,N = n) cdf = rbind(cdf, tempdf) cstart = end = pos n = clust.tbl[i,N] cluster = cluster + 1 } } cdf = rbind(cdf, data.frame(cluster = paste('cluster', cluster, sep='_'), start = cstart, end = end ,N = n)) } else { cdf = data.frame(cluster = 'cluster_1', start = clust.tbl$pos, end = clust.tbl$pos ,N = clust.tbl$N) } #merge adjacent variants to clusters. for(i in 1:nrow(cdf)){ tempcdf = cdf[i,] nonclust.tbl$startDist = nonclust.tbl$pos - tempcdf$start nonclust.tbl$endDist = nonclust.tbl$pos - tempcdf$end merge.adj.to.start = nonclust.tbl[startDist >= -5 & startDist <= 0] if(nrow(merge.adj.to.start) > 0){ tempcdf$start = merge.adj.to.start[which(merge.adj.to.start$startDist == min(merge.adj.to.start$startDist)),pos] tempcdf$N = tempcdf$N + sum(merge.adj.to.start$N) } merge.adj.to.end = nonclust.tbl[endDist <= 5 & endDist >= 0] if(nrow(merge.adj.to.end) > 0){ tempcdf$end = merge.adj.to.end[which(merge.adj.to.end$endDist == max(merge.adj.to.end$endDist)),pos] tempcdf$N = tempcdf$N + sum(merge.adj.to.end$N) } cdf[i,] = tempcdf } cdf$Hugo_Symbol = gene #Calcluate cluster score. total.muts = nrow(prot.dat) #total variants for this gene. clusterScores = c() for(i in 1:nrow(cdf)){ temp.prot.dat = prot.dat[pos >= as.numeric(cdf$start[i]) & pos <= as.numeric(cdf$end[i])] temp.prot.dat.summary = temp.prot.dat[,.N, pos] temp.prot.dat.summary[,fraction:= N/total.muts] peak = temp.prot.dat.summary[N == max(N), pos] posVector = as.numeric(temp.prot.dat.summary[,pos]) fractionMutVector = unlist(lapply(posVector, FUN = function(x) temp.prot.dat.summary[pos == x, fraction])) distanceVector = suppressWarnings(abs(posVector - peak)) clusterScores = c(clusterScores, sum( fractionMutVector / (sqrt(2)^ distanceVector))) } cdf$clusterScore = clusterScores gene.clust.res = data.frame(Hugo_Symbol = gene, clusters = nrow(cdf), muts_in_clusters = sum(cdf$N), clusterScores = sum(cdf$clusterScore), protLen = protLen) return(gene.clust.res) } createOncoMatrix<- function(maf){ message('Creating oncomatrix (this might take a while)..') oncomat = data.table::dcast(data = maf[,.(Hugo_Symbol, Variant_Classification, Tumor_Sample_Barcode)], formula = Hugo_Symbol ~ Tumor_Sample_Barcode, fun.aggregate = function(x) {ifelse(test = length(as.character(x))>1 , no = as.character(x), yes = vcr(x, gis = FALSE)) }, value.var = 'Variant_Classification', fill = '') #If maf contains only one sample converting to matrix is not trivial. if(ncol(oncomat) == 2){ genes = oncomat[,Hugo_Symbol] sampleId = colnames(oncomat)[2] oncomat = as.matrix(data.frame(row.names = genes, sample = oncomat[,2, with =FALSE])) }else if(nrow(oncomat) == 1){ #If MAF has only one gene gene = oncomat[,Hugo_Symbol] oncomat[,Hugo_Symbol:= NULL] oncomat = as.matrix(oncomat) rownames(oncomat) = gene sampleID = colnames(oncomat) }else{ oncomat = as.matrix(oncomat) rownames(oncomat) = oncomat[,1] oncomat = oncomat[,-1] } variant.classes = as.character(unique(maf[,Variant_Classification])) variant.classes = c('',variant.classes, 'Multi_Hit') names(variant.classes) = 0:(length(variant.classes)-1) #Complex variant classes will be assigned a single integer. vc.onc = unique(unlist(apply(oncomat, 2, unique))) vc.onc = vc.onc[!vc.onc %in% names(variant.classes)] names(vc.onc) = rep(as.character(as.numeric(names(variant.classes)[length(variant.classes)])+1), length(vc.onc)) variant.classes2 = c(variant.classes, vc.onc) oncomat.copy <- oncomat #Make a numeric coded matrix for(i in 1:length(variant.classes2)){ oncomat[oncomat == variant.classes2[i]] = names(variant.classes2)[i] } #If maf has only one gene if(nrow(oncomat) == 1){ mdf = t(matrix(as.numeric(oncomat))) rownames(mdf) = gene colnames(mdf) = sampleID return(list(oncomat = oncomat.copy, nummat = mdf, vc = variant.classes)) } #convert from character to numeric mdf = as.matrix(apply(oncomat, 2, function(x) as.numeric(as.character(x)))) rownames(mdf) = rownames(oncomat.copy) message('Sorting..') #If MAF file contains a single sample, simple sorting is enuf. if(ncol(mdf) == 1){ mdf = as.matrix(mdf[order(mdf, decreasing = TRUE),]) colnames(mdf) = sampleId oncomat.copy = as.matrix(oncomat.copy[rownames(mdf),]) colnames(oncomat.copy) = sampleId return(list(oncomat = oncomat.copy, nummat = mdf, vc = variant.classes)) } else{ #Sort by rows as well columns if >1 samples present in MAF #Add total variants per gene mdf = cbind(mdf, variants = apply(mdf, 1, function(x) { length(x[x != "0"]) })) #Sort by total variants mdf = mdf[order(mdf[, ncol(mdf)], decreasing = TRUE), ] colnames(mdf) = gsub(pattern = "^X", replacement = "", colnames(mdf)) nMut = mdf[, ncol(mdf)] mdf = mdf[, -ncol(mdf)] mdf.temp.copy = mdf #temp copy of original unsorted numeric coded matrix mdf[mdf != 0] = 1 #replacing all non-zero integers with 1 improves sorting (& grouping) tmdf = t(mdf) #transposematrix mdf = t(tmdf[do.call(order, c(as.list(as.data.frame(tmdf)), decreasing = TRUE)), ]) #sort mdf.temp.copy = mdf.temp.copy[rownames(mdf),] #organise original matrix into sorted matrix mdf.temp.copy = mdf.temp.copy[,colnames(mdf)] mdf = mdf.temp.copy #organise original character matrix into sorted matrix oncomat.copy <- oncomat.copy[,colnames(mdf)] oncomat.copy <- oncomat.copy[rownames(mdf),] return(list(oncomat = oncomat.copy, nummat = mdf, vc = variant.classes)) } } validateMaf<-function(maf, rdup = TRUE, isTCGA = isTCGA){ #necessary fields. required.fields = c('Hugo_Symbol', 'Chromosome', 'Start_Position', 'End_Position', 'Reference_Allele', 'Tumor_Seq_Allele2', 'Variant_Classification', 'Variant_Type', 'Tumor_Sample_Barcode') #Change column names to standard names; i.e, camel case for(i in 1:length(required.fields)){ colId = suppressWarnings(grep(pattern = required.fields[i], x = colnames(maf), ignore.case = TRUE)) if(length(colId) > 0){ colnames(maf)[colId] = required.fields[i] } } missing.fileds = required.fields[!required.fields %in% colnames(maf)] #check if any of them are missing if(length(missing.fileds) > 0){ missing.fileds = paste(missing.fileds[1], sep = ',', collapse = ', ') stop(paste('missing required fields from MAF:', missing.fileds)) #stop if any of required.fields are missing } #convert "-" to "." in "Tumor_Sample_Barcode" to avoid complexity in naming maf$Tumor_Sample_Barcode = gsub(pattern = '-', replacement = '.', x = as.character(maf$Tumor_Sample_Barcode)) if(rdup){ maf = maf[, variantId := paste(Chromosome, Start_Position, Tumor_Sample_Barcode, sep = ':')] if(nrow(maf[duplicated(variantId)]) > 0){ message("NOTE: Removed ", nrow(maf[duplicated(variantId)]) ," duplicated variants") maf = maf[!duplicated(variantId)] } maf[,variantId := NULL] } if(nrow(maf[Hugo_Symbol %in% ""]) > 0){ message('NOTE: Found ', nrow(maf[Hugo_Symbol %in% ""]), ' variants with no Gene Symbols.') print(maf[Hugo_Symbol %in% "", required.fields, with = FALSE]) message("Annotating them as 'UnknownGene' for convenience") maf$Hugo_Symbol = ifelse(test = maf$Hugo_Symbol == "", yes = 'UnknownGene', no = maf$Hugo_Symbol) } if(nrow(maf[is.na(Hugo_Symbol)]) > 0){ message('NOTE: Found ', nrow(maf[is.na(Hugo_Symbol) > 0]), ' variants with no Gene Symbols.') print(maf[is.na(Hugo_Symbol), required.fields, with =FALSE]) message("Annotating them as 'UnknownGene' for convenience") maf$Hugo_Symbol = ifelse(test = is.na(maf$Hugo_Symbol), yes = 'UnknownGene', no = maf$Hugo_Symbol) } if(isTCGA){ maf$Tumor_Sample_Barcode = substr(x = maf$Tumor_Sample_Barcode, start = 1, stop = 12) } return(maf) } read.maf_fix = function(maf, removeSilent = TRUE, useAll = TRUE, gisticAllLesionsFile = NULL, gisticAmpGenesFile = NULL, gisticDelGenesFile = NULL, cnTable = NULL, removeDuplicatedVariants = TRUE, isTCGA = FALSE){ message('reading maf..') if(as.logical(length(grep(pattern = 'gz$', x = maf, fixed = FALSE)))){ #If system is Linux use fread, else use gz connection to read gz file. if(Sys.info()[['sysname']] == 'Windows'){ maf.gz = gzfile(description = maf, open = 'r') suppressWarnings(maf <- data.table(read.csv(file = maf.gz, header = TRUE, sep = '\t', stringsAsFactors = FALSE))) close(maf.gz) } else{ maf = suppressWarnings(data.table::fread(input = paste('zcat <', maf), sep = '\t', stringsAsFactors = FALSE, verbose = FALSE, data.table = TRUE, showProgress = TRUE, header = TRUE)) } } else{ suppressWarnings(maf <- data.table::fread(input = maf, sep = "\t", stringsAsFactors = FALSE, verbose = FALSE, data.table = TRUE, showProgress = TRUE, header = TRUE)) } #validate MAF file maf = validateMaf(maf = maf, isTCGA = isTCGA, rdup = removeDuplicatedVariants) #validation check for variants classified as Somatic in Mutation_Status field. if(length(colnames(maf)[colnames(x = maf) %in% 'Mutation_Status']) > 0){ if(!useAll){ message('Using only Somatic variants from Mutation_Status. Switch on useAll to include everything.') maf = maf[Mutation_Status %in% "Somatic"] if(nrow(maf) == 0){ stop('No more Somatic mutations left after filtering for Mutation_Status! Maybe set useAll to TRUE ?') } #maf = subset(maf, Mutation_Status == 'Somatic') }else { message('Using all variants.') } }else{ message('Mutation_Status not found. Assuming all variants are Somatic and validated.') } #Variant Classification with Low/Modifier variant consequences. http://asia.ensembl.org/Help/Glossary?id=535 silent = c("3'UTR", "5'UTR", "3'Flank", "Targeted_Region", "Silent", "Intron", "RNA", "IGR", "Splice_Region", "5'Flank", "lincRNA") #Variant Classification with High/Moderate variant consequences. http://asia.ensembl.org/Help/Glossary?id=535 vc.nonSilent = c("Frame_Shift_Del", "Frame_Shift_Ins", "Splice_Site", "Translation_Start_Site", "Nonsense_Mutation", "Nonstop_Mutation", "In_Frame_Del", "In_Frame_Ins", "Missense_Mutation") maf.silent = maf[Variant_Classification %in% silent] if(removeSilent){ if(nrow(maf.silent) > 0){ maf.silent.vc = maf.silent[,.N, .(Tumor_Sample_Barcode, Variant_Classification)] maf.silent.vc.cast = data.table::dcast(data = maf.silent.vc, formula = Tumor_Sample_Barcode ~ Variant_Classification, fill = 0, value.var = 'N') #why dcast is not returning it as data.table ? summary.silent = data.table(ID = c('Samples',colnames(maf.silent.vc.cast)[2:ncol(maf.silent.vc.cast)]), N = c(nrow(maf.silent.vc.cast), colSums(maf.silent.vc.cast[,2:ncol(maf.silent.vc.cast), with = FALSE]))) maf = maf[!Variant_Classification %in% silent] #Remove silent variants from main table message(paste('Excluding',nrow(maf.silent), 'silent variants.')) print(summary.silent) } else{ message(message(paste('Excluding',nrow(maf.silent), 'silent variants.'))) } }else{ message('Silent variants are being kept!') } if(!is.null(gisticAllLesionsFile)){ gisticIp = readGistic(gisticAllLesionsFile = gisticAllLesionsFile, gisticAmpGenesFile = gisticAmpGenesFile, gisticDelGenesFile = gisticDelGenesFile, isTCGA = isTCGA) gisticIp = gisticIp@data gisticIp[, id := paste(Hugo_Symbol, Tumor_Sample_Barcode, sep=':')] gisticIp = gisticIp[!duplicated(id)] gisticIp[,id := NULL] maf = rbind(maf, gisticIp, fill =TRUE) oncomat = createOncoMatrix(maf) }else if(!is.null(cnTable)){ message('Processing copy number data..') cnDat = data.table::fread(input = cnTable, sep = '\t', stringsAsFactors = FALSE, header = TRUE, colClasses = 'character') colnames(cnDat) = c('Hugo_Symbol', 'Tumor_Sample_Barcode', 'Variant_Classification') cnDat$Variant_Type = 'CNV' suppressWarnings(cnDat[, id := paste(Hugo_Symbol, Tumor_Sample_Barcode, sep=':')]) cnDat = cnDat[!duplicated(id)] cnDat[,id := NULL] maf = rbind(maf, cnDat, fill =TRUE) oncomat = createOncoMatrix(maf) }else{ oncomat = createOncoMatrix(maf) } #convert to factors maf$Variant_Type = as.factor(as.character(maf$Variant_Type)) maf$Variant_Classification = as.factor(as.character(maf$Variant_Classification)) maf$Tumor_Sample_Barcode = as.factor(as.character(maf$Tumor_Sample_Barcode)) message('Summarizing..') mafSummary = summarizeMaf_fix(maf = maf) #Create MAF object m = MAF(data = maf, variants.per.sample = mafSummary$variants.per.sample, variant.type.summary = mafSummary$variant.type.summary, variant.classification.summary = mafSummary$variant.classification.summary,gene.summary = mafSummary$gene.summary, oncoMatrix = oncomat$oncomat, numericMatrix = oncomat$nummat, summary = mafSummary$summary, classCode = oncomat$vc, maf.silent = maf.silent) message('Done !') return(m) } #' Class MAF #' @description S4 class for storing summarized MAF. #' @slot data data.table of original MAF file. #' @slot variants.per.sample table containing variants per sample #' @slot variant.type.summary table containing variant types per sample #' @slot variant.classification.summary table containing variant classification per sample #' @slot gene.summary table containing variant classification per gene #' @slot oncoMatrix character matrix of dimension n*m where n is number of genes and m is number of variants #' @slot numericMatrix numeric matrix of dimension n*m where n is number of genes and m is number of variants #' @slot summary table with basic MAF summary stats #' @slot classCode mapping between numeric values in numericMatrix and Variant Classification #' @slot maf.silent subset of main MAF containing only silent variants #' @exportClass MAF #' @import methods #' @seealso \code{\link{getGeneSummary}} \code{\link{getSampleSummary}} \code{\link{getFields}} ## MAF object MAF <- setClass(Class = 'MAF', slots = c(data = 'data.table', variants.per.sample = 'data.table', variant.type.summary = 'data.table', variant.classification.summary = 'data.table', gene.summary = 'data.table', oncoMatrix = 'matrix', numericMatrix = 'matrix', summary = 'data.table', classCode = 'character', maf.silent = 'data.table')) setMethod(f = 'show', signature = 'MAF', definition = function(object){ cat(paste('An object of class ', class(object), "\n")) print(object@summary) }) summarizeMaf_fix = function(maf){ if('NCBI_Build' %in% colnames(maf)){ NCBI_Build = unique(maf[!Variant_Type %in% 'CNV', NCBI_Build]) NCBI_Build = NCBI_Build[!is.na(NCBI_Build)] if(length(NCBI_Build) > 1){ message('NOTE: Mutiple reference builds found!') NCBI_Build = do.call(paste, c(as.list(NCBI_Build), sep=";")) message(NCBI_Build) } }else{ NCBI_Build = NA } if('Center' %in% colnames(maf)){ Center = unique(maf[!Variant_Type %in% 'CNV', Center]) #Center = Center[is.na(Center)] if(length(Center) > 1){ message('Mutiple centers found.') Center = do.call(paste, c(as.list(Center), sep=";")) print(Center) } }else{ Center = NA } #nGenes nGenes = length(unique(maf[,Hugo_Symbol])) #Top 20 FLAGS - https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4267152/ flags = c("TTN", "MUC16", "OBSCN", "AHNAK2", "SYNE1", "FLG", "MUC5B", "DNAH17", "PLEC", "DST", "SYNE2", "NEB", "HSPG2", "LAMA5", "AHNAK", "HMCN1", "USH2A", "DNAH11", "MACF1", "MUC17") #Variants per TSB tsb = maf[,.N, Tumor_Sample_Barcode] colnames(tsb)[2] = 'Variants' tsb = tsb[order(tsb$Variants, decreasing = TRUE),] #summarise and casting by 'Variant_Classification' vc = maf[,.N, .(Tumor_Sample_Barcode, Variant_Classification )] vc.cast = data.table::dcast(data = vc, formula = Tumor_Sample_Barcode ~ Variant_Classification, fill = 0, value.var = 'N') if(any(colnames(vc.cast) %in% c('Amp', 'Del'))){ vc.cast.cnv = vc.cast[,colnames(vc.cast)[colnames(vc.cast) %in% c('Amp', 'Del')], with =FALSE] vc.cast.cnv$CNV_total = rowSums(x = vc.cast.cnv) vc.cast = vc.cast[,!colnames(vc.cast)[colnames(vc.cast) %in% c('Amp', 'Del')], with =FALSE] vc.cast[,total:=rowSums(vc.cast[,2:ncol(vc.cast), with = FALSE])] vc.cast = cbind(vc.cast, vc.cast.cnv) vc.cast = vc.cast[order(total, CNV_total, decreasing = TRUE)] vc.mean = as.numeric(as.character(c(NA, NA, NA, NA, apply(vc.cast[,2:ncol(vc.cast), with = FALSE], 2, mean)))) vc.median = as.numeric(as.character(c(NA, NA, NA, NA, apply(vc.cast[,2:ncol(vc.cast), with = FALSE], 2, median)))) }else{ vc.cast[,total:=rowSums(vc.cast[,2:ncol(vc.cast), with = FALSE])] vc.cast = vc.cast[order(total, decreasing = TRUE)] vc.mean = as.numeric(as.character(c(NA, NA, NA, NA, apply(vc.cast[,2:ncol(vc.cast), with = FALSE], 2, mean)))) vc.median = as.numeric(as.character(c(NA, NA, NA, NA, apply(vc.cast[,2:ncol(vc.cast), with = FALSE], 2, median)))) } #summarise and casting by 'Variant_Type' vt = maf[,.N, .(Tumor_Sample_Barcode, Variant_Type )] vt.cast = data.table::dcast(data = vt, formula = Tumor_Sample_Barcode ~ Variant_Type, value.var = 'N', fill = 0) if(any(colnames(vt.cast) %in% c('CNV'))){ vt.cast.cnv = vt.cast[,colnames(vt.cast)[colnames(vt.cast) %in% c('CNV')], with =FALSE] vt.cast = vt.cast[,!colnames(vt.cast)[colnames(vt.cast) %in% c('CNV')], with =FALSE] vt.cast[,total:=rowSums(vt.cast[,2:ncol(vt.cast), with = FALSE])] vt.cast = vt.cast[order(total, decreasing = TRUE)] vt.cast = cbind(vt.cast, vt.cast.cnv) vt.cast[order(total, CNV, decreasing = TRUE)] }else{ vt.cast[,total:=rowSums(vt.cast[,2:ncol(vt.cast), with = FALSE])] vt.cast = vt.cast[order(total, decreasing = TRUE)] } #summarise and casting by 'Hugo_Symbol' hs = maf[,.N, .(Hugo_Symbol, Variant_Classification)] hs.cast = data.table::dcast(data = hs, formula = Hugo_Symbol ~Variant_Classification, fill = 0, value.var = 'N') #---- if(any(colnames(hs.cast) %in% c('Amp', 'Del'))){ hs.cast.cnv = hs.cast[,colnames(hs.cast)[colnames(hs.cast) %in% c('Amp', 'Del')], with =FALSE] hs.cast.cnv$CNV_total = rowSums(x = hs.cast.cnv) hs.cast = hs.cast[,!colnames(hs.cast)[colnames(hs.cast) %in% c('Amp', 'Del')], with =FALSE] hs.cast[,total:=rowSums(hs.cast[,2:ncol(hs.cast), with = FALSE])] hs.cast = cbind(hs.cast, hs.cast.cnv) hs.cast = hs.cast[order(total, CNV_total, decreasing = TRUE)] }else{ hs.cast[,total:=rowSums(hs.cast[,2:ncol(hs.cast), with = FALSE])] hs.cast = hs.cast[order(total, decreasing = TRUE)] } #Get in how many samples a gene ismutated numMutatedSamples = maf[!Variant_Type %in% 'CNV', .(MutatedSamples = length(unique(Tumor_Sample_Barcode))), by = Hugo_Symbol] #Merge and sort hs.cast = merge(hs.cast, numMutatedSamples, by = 'Hugo_Symbol', all = TRUE) hs.cast = hs.cast[order(MutatedSamples, total, decreasing = TRUE)] #Make a summarized table summary = data.table::data.table(ID = c('NCBI_Build', 'Center','Samples', 'nGenes',colnames(vc.cast)[2:ncol(vc.cast)]), summary = c(NCBI_Build, Center, nrow(vc.cast), nGenes, colSums(vc.cast[,2:ncol(vc.cast), with =FALSE]))) summary[,Mean := vc.mean] summary[,Median := vc.median] print(summary) message("Frequently mutated genes..") print(hs.cast) #Check for flags. if(nrow(hs.cast) > 10){ topten = hs.cast[1:10, Hugo_Symbol] topten = topten[topten %in% flags] if(length(topten) > 0){ message('NOTE: Possible FLAGS among top ten genes:') print(topten) } } return(list(variants.per.sample = tsb, variant.type.summary = vt.cast, variant.classification.summary = vc.cast, gene.summary = hs.cast, summary = summary)) } oncodrive_fix = function(maf, AACol = NULL, minMut = 5, pvalMethod = 'zscore', nBgGenes = 100, bgEstimate = TRUE, ignoreGenes = NULL){ #Proetin Length source gl = system.file('extdata', 'prot_len.txt.gz', package = 'maftools') if(Sys.info()[['sysname']] == 'Windows'){ gl.gz = gzfile(description = gl, open = 'r') gl <- suppressWarnings( data.table(read.csv( file = gl.gz, header = TRUE, sep = '\t', stringsAsFactors = FALSE)) ) close(gl.gz) } else{ gl = data.table::fread(input = paste('zcat <', gl), sep = '\t', stringsAsFactors = FALSE) } pval.options = c('zscore', 'poisson', 'combined') if(!pvalMethod %in% pval.options){ stop('pvalMethod can only be either zscore, poisson or combined') } if(length(pvalMethod) > 1){ stop('pvalMethod can only be either zscore, poisson or combined') } #syn variants for background syn.maf = maf@maf.silent #number of samples in maf numSamples = as.numeric(maf@summary[3,summary]) #Perform clustering and calculate background scores. if(bgEstimate){ if(nrow(syn.maf) == 0){ message('No syn mutations found! Skipping background estimation. Using predefined values. (Mean = 0.279; SD = 0.13)') bg.mean = 0.279 bg.sd = 0.13 }else{ message('Estimating background scores from synonymous variants..') syn.bg.scores = parse_prot_fix(dat = syn.maf, AACol = AACol, gl, m = minMut, calBg = TRUE, nBg = nBgGenes) #If number of genes to calculate background scores is not enough, use predefined scores. if(is.null(syn.bg.scores)){ message("Not enough genes to build background. Using predefined values. (Mean = 0.279; SD = 0.13)") bg.mean = 0.279 bg.sd = 0.13 }else { if(nrow(syn.bg.scores) < nBgGenes){ message("Not enough genes to build background. Using predefined values. (Mean = 0.279; SD = 0.13)") bg.mean = 0.279 bg.sd = 0.13 }else{ bg.mean = mean(syn.bg.scores$clusterScores) bg.sd = sd(syn.bg.scores$clusterScores) message(paste('Estimated background mean: ', bg.mean)) message(paste('Estimated background SD: ', bg.sd)) } } } }else{ message("Using predefined values for background. (Mean = 0.279; SD = 0.13)") bg.mean = 0.279 bg.sd = 0.13 } #non-syn variants non.syn.maf = maf@data #Variant Classification with Low/Modifier variant consequences. http://asia.ensembl.org/Help/Glossary?id=535 silent = c("3'UTR", "5'UTR", "3'Flank", "Targeted_Region", "Silent", "Intron", "RNA", "IGR", "Splice_Region", "5'Flank", "lincRNA", "Amp", "Del") non.syn.maf = non.syn.maf[!Variant_Classification %in% silent] #Remove silent variants from main table #Remove genes to ignore if(!is.null(ignoreGenes)){ ignoreGenes.count = nrow(non.syn.maf[Hugo_Symbol %in% ignoreGenes]) message(paste('Removed', ignoreGenes.count, 'variants belonging to', paste(ignoreGenes, collapse = ', ', sep=','))) non.syn.maf = non.syn.maf[!Hugo_Symbol %in% ignoreGenes] } #Perform clustering and calculate cluster scores for nonsyn variants. message('Estimating cluster scores from non-syn variants..') nonsyn.scores = parse_prot_fix(dat = non.syn.maf, AACol = AACol, gl = gl, m = minMut, calBg = FALSE, nBg = nBgGenes) if(pvalMethod == 'combined'){ message('Comapring with background model and estimating p-values..') nonsyn.scores$zscore = (nonsyn.scores$clusterScores - bg.mean) / bg.sd nonsyn.scores$tPval = 1- pnorm(nonsyn.scores$zscore) nonsyn.scores$tFdr = p.adjust(nonsyn.scores$tPval, method = 'fdr') nonsyn.scores = merge(getGeneSummary(maf), nonsyn.scores, by = 'Hugo_Symbol') nonsyn.scores[,fract_muts_in_clusters := muts_in_clusters/total] counts.glm = glm(formula = total ~ protLen+clusters, family = poisson(link = identity), data = nonsyn.scores) #Poisson model nonsyn.scores$Expected = counts.glm$fitted.values #Get expected number of events (mutations) from the model observed_mut_colIndex = which(colnames(nonsyn.scores) == 'total') expected_mut_colIndex = which(colnames(nonsyn.scores) == 'Expected') #Poisson test to caluclate difference (p-value) nonsyn.scores$poissonPval = apply(nonsyn.scores, 1, function(x) { poisson.test(as.numeric(x[observed_mut_colIndex]), as.numeric(x[expected_mut_colIndex]))$p.value }) nonsyn.scores$poissonFdr = p.adjust(nonsyn.scores$poissonPval) nonsyn.scores = nonsyn.scores[order(poissonFdr)] nonsyn.scores$fdr = apply(nonsyn.scores[,.(tFdr, poissonFdr)], MARGIN = 1, FUN = min) } else if(pvalMethod == 'zscore'){ #Oncodrive clust way of caluclating pvalues #Calculate z scores; compare it to bg scores and estimate z-score, pvalues, corrected pvalues (fdr) (assumes normal distribution) message('Comapring with background model and estimating p-values..') nonsyn.scores$zscore = (nonsyn.scores$clusterScores - bg.mean) / bg.sd nonsyn.scores$pval = 1- pnorm(nonsyn.scores$zscore) nonsyn.scores$fdr = p.adjust(nonsyn.scores$pval, method = 'fdr') nonsyn.scores = merge(getGeneSummary(maf), nonsyn.scores, by = 'Hugo_Symbol') nonsyn.scores[,fract_muts_in_clusters := muts_in_clusters/total] #nonsyn.scores[,fract_MutatedSamples := MutatedSamples/numSamples] nonsyn.scores = nonsyn.scores[order(fdr)] }else{ #Assuming poisson distribution of mutation counts #Now model observed number of mutations as a function of number of clusters and protein length. Calculate expected number of events based on poisson distribution. nonsyn.scores = merge(getGeneSummary(maf), nonsyn.scores, by = 'Hugo_Symbol') nonsyn.scores[,fract_muts_in_clusters := muts_in_clusters/total] counts.glm = glm(formula = total ~ protLen+clusters, family = poisson(link = identity), data = nonsyn.scores) #Poisson model nonsyn.scores$Expected = counts.glm$fitted.values #Get expected number of events (mutations) from the model observed_mut_colIndex = which(colnames(nonsyn.scores) == 'total') expected_mut_colIndex = which(colnames(nonsyn.scores) == 'Expected') #Poisson test to caluclate difference (p-value) nonsyn.scores$pval = apply(nonsyn.scores, 1, function(x) { poisson.test(as.numeric(x[observed_mut_colIndex]), as.numeric(x[expected_mut_colIndex]))$p.value }) nonsyn.scores$fdr = p.adjust(nonsyn.scores$pval) nonsyn.scores = nonsyn.scores[order(fdr)] } message('Done !') return(nonsyn.scores) } laml = read.maf(maf = args$input_maf, removeSilent = F, useAll = T) if(is.null(args$gene_blacklist)){ laml.sig = oncodrive(maf =laml, AACol = aacol, pvalMethod = 'zscore',minMut = min_mut) write.table(laml.sig,file=args$output_detail, quote=FALSE,row.names=FALSE,sep="\t") pdf(args$output_plot) plotOncodrive(res=laml.sig,fdrCutOff=as.numeric(args$fdr),useFraction=TRUE) dev.off() }else{ all_genes <- read.table(args$gene_blacklist, stringsAsFactors=FALSE)[,1] laml.sig = oncodrive(maf =laml, AACol = aacol, pvalMethod = 'zscore',minMut = min_mut,ignoreGenes=all_genes) write.table(laml.sig,file=args$output_detail, quote=FALSE,row.names=FALSE,sep="\t") pdf(args$output_plot) plotOncodrive(res=laml.sig,fdrCutOff=as.numeric(args$fdr),useFraction=TRUE) dev.off() }