R/oncoplot.R
In PoisonAlien/maftools: Summarize, Analyze and Visualize MAF Files
Defines functions
oncoplot
Documented in
oncoplot
#' draw an oncoplot
#' @description takes output generated by read.maf and draws an oncoplot
#'
#' @details
#' Takes an \code{\link{MAF}} object as an input and plots it as a matrix. Any desired clincal features can be added at the bottom of the oncoplot by providing \code{clinicalFeatures}.
#' Oncoplot can be sorted either by mutations or by clinicalFeatures using arguments \code{sortByMutation} and \code{sortByAnnotation} respectively.
#'
#'By setting `pathways` argument either `sigpw` or `smgbp` - cohort can be summarized by altered pathways. pathways argument also accepts a custom pathway list in the form of a two column tsv file or a data.frame containing gene names and their corresponding pathway.
#'
#' @references Bailey, Matthew H et al. “Comprehensive Characterization of Cancer Driver Genes and Mutations.” Cell vol. 173,2 (2018): 371-385.e18. doi:10.1016/j.cell.2018.02.060
#' Sanchez-Vega, Francisco et al. “Oncogenic Signaling Pathways in The Cancer Genome Atlas.” Cell vol. 173,2 (2018): 321-337.e10. doi:10.1016/j.cell.2018.03.035
#'
#' @param maf an \code{\link{MAF}} object generated by \code{\link{read.maf}}
#' @param top how many top genes to be drawn. defaults to 20.
#' @param minMut draw all genes with `min` number of mutations. Can be an integer or fraction (of samples mutated), Default NULL
#' @param genes Just draw oncoplot for these genes. Default NULL.
#' @param altered Default FALSE. Chooses top genes based on muatation status. If \code{TRUE} chooses top genes based alterations (CNV or mutation).
#' @param drawColBar logical plots top barplot for each sample. Default \code{TRUE}.
#' @param drawRowBar logical. Plots righ barplot for each gene. Default \code{TRUE}.
#' @param leftBarData Data for leftside barplot. Must be a data.frame with two columns containing gene names and values. Default `NULL`
#' @param leftBarLims limits for `leftBarData`. Default `NULL`.
#' @param leftBarVline Draw vertical lines at these values. Default `NULL`.
#' @param leftBarVlineCol Line color for `leftBarVline` Default gray70
#' @param topBarData Default `NULL` which draws absolute number of mutation load for each sample. Can be overridden by choosing one clinical indicator(Numeric) or by providing a two column data.frame containing sample names and values for each sample. This option is applicable when only `drawColBar` is TRUE.
#' @param topBarLims limits for `topBarData`. Default `NULL`.
#' @param topBarHline Draw horizontal lines at these values. Default `NULL`.
#' @param topBarHlineCol Line color for `topBarHline.` Default gray70
#' @param rightBarData Data for rightside barplot. Must be a data.frame with two columns containing to gene names and values. Default `NULL` which draws distribution by variant classification. This option is applicable when only `drawRowBar` is TRUE.
#' @param rightBarLims limits for `rightBarData`. Default `NULL`.
#' @param rightBarVline Draw vertical lines at these values. Default `NULL`.
#' @param rightBarVlineCol Line color for `rightBarVline` Default gray70
#' @param logColBar Plot top bar plot on log10 scale. Default \code{FALSE}.
#' @param includeColBarCN Whether to include CN in column bar plot. Default TRUE
#' @param clinicalFeatures columns names from `clinical.data` slot of \code{MAF} to be drawn in the plot. Default NULL.
#' @param annotationColor Custom colors to use for `clinicalFeatures`. Must be a named list containing a named vector of colors. Default NULL. See example for more info.
#' @param annotationDat If MAF file was read without clinical data, provide a custom \code{data.frame} with a column \code{Tumor_Sample_Barcode} containing sample names along with rest of columns with annotations.
#' You can specify which columns to be drawn using `clinicalFeatures` argument.
#' @param pathways Default `NULL`. Can be `sigpw`, `smgbp`, or a two column data.frame/tsv-file with genes and corresponding pathway mappings.`
#' @param topPathways Top most altered pathways to draw. Default 3. Mutually exclusive with `selectedPathways`
#' @param path_order Default `NULL` Manually specify the order of pathways
#' @param selectedPathways Manually provide the subset of pathway names to be selected from `pathways`. Default NULL. In case `pathways` is `auto` draws top 3 altered pathways.
#' @param collapsePathway Shows only rows corresponding to the pathways. Default FALSE.
#' @param pwLineCol Color for the box around the pathways Default #535c68
#' @param pwLineWd Line width for the box around the pathways Default Default 1
#' @param draw_titv logical Includes TiTv plot. \code{FALSE}
#' @param titv_col named vector of colors for each transition and transversion classes. Should be of length six with the names "C>T" "C>G" "C>A" "T>A" "T>C" "T>G". Default NULL.
#' @param showTumorSampleBarcodes logical to include sample names.
#' @param tsbToPIDs Custom names for Tumor_Sample_Barcodes. Can be a column name in clinicaldata or a 2 column data.frame of Tumor_Sample_Barcodes to patient ID mappings. Applicable only when `showTumorSampleBarcodes = TRUE`. Default NULL.
#' @param barcode_mar Margin width for sample names. Default 4
#' @param barcodeSrt Rotate sample labels. Default 90.
#' @param gene_mar Margin width for gene names. Default 5
#' @param anno_height Height of plotting area for sample annotations. Default 1
#' @param legend_height Height of plotting area for legend. Default 4
#' @param sortByAnnotation logical sort oncomatrix (samples) by provided `clinicalFeatures`. Sorts based on first `clinicalFeatures`. Defaults to FALSE. column-sort
#' @param groupAnnotationBySize Further group `sortByAnnotation` orders by their size. Defaults to TRUE. Largest groups comes first.
#' @param annotationOrder Manually specify order for annotations. Works only for first `clinicalFeatures`. Default NULL.
#' @param sortByMutation Force sort matrix according mutations. Helpful in case of MAF was read along with copy number data. Default FALSE.
#' @param keepGeneOrder logical whether to keep order of given genes. Default FALSE, order according to mutation frequency
#' @param GeneOrderSort logical this is applicable when `keepGeneOrder` is TRUE. Default TRUE
#' @param sampleOrder Manually speify sample names for oncolplot ordering. Default NULL.
#' @param additionalFeature a vector of length two indicating column name in the MAF and the factor level to be highlighted. Provide a list of values for highlighting more than one features
#' @param additionalFeaturePch Default 20
#' @param additionalFeatureCol Default "gray70"
#' @param additionalFeatureCex Default 0.9
#' @param genesToIgnore do not show these genes in Oncoplot. Default NULL.
#' @param removeNonMutated Logical. If \code{TRUE} removes samples with no mutations in the oncoplot for better visualization. Default \code{FALSE}.
#' @param fill Logical. If \code{TRUE} draws genes and samples as blank grids even when they are not altered.
#' @param cohortSize Number of sequenced samples in the cohort. Default all samples from Cohort. You can manually specify the cohort size. Default \code{NULL}
#' @param colors named vector of colors for each Variant_Classification.
#' @param cBioPortal Adds annotations similar to cBioPortals MutationMapper and collapse Variants into Truncating and rest.
#' @param bgCol Background grid color for wild-type (not-mutated) samples. Default "#ecf0f1"
#' @param borderCol border grid color (not-mutated) samples. Default 'white'.
#' @param annoBorderCol border grid color for annotations. Default NA.
#' @param numericAnnoCol color palette used for numeric annotations. Default 'YlOrBr' from RColorBrewer
#' @param drawBox logical whether to draw a box around main matrix. Default FALSE
#' @param fontSize font size for gene names. Default 0.8.
#' @param SampleNamefontSize font size for sample names. Default 1
#' @param titleFontSize font size for title. Default 1.5
#' @param legendFontSize font size for legend. Default 1.2
#' @param annotationFontSize font size for annotations. Default 1.2
#' @param sepwd_genes size of lines seperating genes. Default 0.5
#' @param sepwd_samples size of lines seperating samples. Default 0.25
#' @param writeMatrix writes character coded matrix used to generate the plot to an output file.
#' @param colbar_pathway Draw top column bar with respect to diplayed pathway. Default FALSE.
#' @param showTitle Default TRUE
#' @param titleText Custom title. Default `NULL`
#' @param showPct Default TRUE. Shows percent altered to the right side of the plot.
#' @returns Invisibly returns a list with components
#' 1. `oncomatrix` A matrix used for drawing the oncoplot. Values are numeric coded for each variant classification
#' 2. `vc_legend` A mapping of variant classification to numeric values in the oncomatrix
#' 3. `vc_color` Color coding used for each variant classification
#' @examples
#' laml.maf <- system.file("extdata", "tcga_laml.maf.gz", package = "maftools")
#' laml.clin = system.file('extdata', 'tcga_laml_annot.tsv', package = 'maftools')
#' laml <- read.maf(maf = laml.maf, clinicalData = laml.clin)
#' #Basic onocplot
#' oncoplot(maf = laml, top = 3)
#' #Changing colors for variant classifications (You can use any colors, here in this example we will use a color palette from RColorBrewer)
#' col = RColorBrewer::brewer.pal(n = 8, name = 'Paired')
#' names(col) = c('Frame_Shift_Del','Missense_Mutation', 'Nonsense_Mutation', 'Multi_Hit', 'Frame_Shift_Ins',
#' 'In_Frame_Ins', 'Splice_Site', 'In_Frame_Del')
#' #Color coding for FAB classification; try getAnnotations(x = laml) to see available annotations.
#' fabcolors = RColorBrewer::brewer.pal(n = 8,name = 'Spectral')
#' names(fabcolors) = c("M0", "M1", "M2", "M3", "M4", "M5", "M6", "M7")
#' fabcolors = list(FAB_classification = fabcolors)
#' oncoplot(maf = laml, colors = col, clinicalFeatures = 'FAB_classification', sortByAnnotation = TRUE, annotationColor = fabcolors)
#' @seealso \code{\link{pathways}}
#' @export
oncoplot = oncoplot = function(maf, top = 20, minMut = NULL, genes = NULL, altered = FALSE,
drawRowBar = TRUE, drawColBar = TRUE,
leftBarData = NULL, leftBarLims = NULL, leftBarVline = NULL, leftBarVlineCol = 'gray70',
rightBarData = NULL, rightBarLims = NULL,rightBarVline = NULL, rightBarVlineCol = 'gray70',
topBarData = NULL, topBarLims = NULL, topBarHline = NULL, topBarHlineCol = 'gray70', logColBar = FALSE, includeColBarCN = TRUE,
clinicalFeatures = NULL, annotationColor = NULL, annotationDat = NULL,
pathways = NULL, topPathways = 3,path_order = NULL, selectedPathways = NULL, collapsePathway = FALSE, pwLineCol = "#535c68", pwLineWd = 1, draw_titv = FALSE, titv_col = NULL,
showTumorSampleBarcodes = FALSE, tsbToPIDs = NULL, barcode_mar = 4, barcodeSrt = 90, gene_mar = 5,
anno_height = 1, legend_height = 4,
sortByAnnotation = FALSE, groupAnnotationBySize = TRUE, annotationOrder = NULL,
sortByMutation = FALSE, keepGeneOrder = FALSE, GeneOrderSort = TRUE, sampleOrder = NULL,
additionalFeature = NULL, additionalFeaturePch = 20, additionalFeatureCol = "gray70", additionalFeatureCex = 0.9,
genesToIgnore = NULL, removeNonMutated = FALSE, fill = TRUE, cohortSize = NULL,
colors = NULL, cBioPortal = FALSE, bgCol = "#ecf0f1", borderCol = 'white', annoBorderCol = NA, numericAnnoCol = NULL,
drawBox = FALSE, fontSize = 0.8, SampleNamefontSize = 1, titleFontSize = 1.5, legendFontSize = 1.2, annotationFontSize = 1.2,
sepwd_genes = 0.5, sepwd_samples = 0.25, writeMatrix = FALSE, colbar_pathway = FALSE, showTitle = TRUE, titleText = NULL, showPct = TRUE){
#Total samples
if(is.null(cohortSize)){
totSamps = as.numeric(maf@summary[3,summary])
}else{
totSamps = cohortSize
}
if(!is.null(genes)){ #If user provides a gene list
genes = unique(as.character(genes))
om = createOncoMatrix(m = maf, g = genes, add_missing = fill, cbio = cBioPortal)
numMat = om$numericMatrix
mat_origin = om$oncoMatrix
if(!is.null(pathways)){
stop("Please use genes or pathways.")
}
} else if(!is.null(pathways)){ #If user provides pathway list
if(is(pathways, 'data.frame')){
pathways = data.table::copy(pathways)
colnames(pathways)[1:2] = c('Gene', 'Pathway')
data.table::setDT(x = pathways)
}else if(file.exists(pathways)){
pathways = data.table::fread(file = pathways)
colnames(pathways)[1:2] = c('Gene', 'Pathway')
data.table::setDT(x = pathways)
}else{
pathways = match.arg(arg = pathways, choices = c("sigpw", "smgbp"))
if(pathways == "sigpw"){
pathdb <- system.file("extdata", "oncogenic_sig_patwhays.tsv", package = "maftools")
#message("Summarizing signalling pathways [Sanchez-Vega et al., https://doi.org/10.1016/j.cell.2018.03.035]")
}else{
#message("Summarizing known driver genes [Bailey et al., https://doi.org/10.1016/j.cell.2018.02.060]")
pathdb <- system.file("extdata", "BP_SMGs.txt.gz", package = "maftools")
}
pathways = data.table::fread(file = pathdb, skip = "Gene")
}
#Convert to characters in case of factors
pathways$Gene = as.character(pathways$Gene)
pathways$Pathway = as.character(pathways$Pathway)
if(is.null(selectedPathways)){
pathwayLoad = get_pw_summary(maf = maf, pathways = pathways[,.(Pathway, Gene)])
message("Drawing upto top ", topPathways, " mutated pathways")
if(ncol(pathwayLoad) >= topPathways){
pathways = pathways[Pathway %in% pathwayLoad[1:topPathways, Pathway]]
}else{
pathways = pathways[Pathway %in% pathwayLoad[, Pathway]]
}
}else{
pws_missing = setdiff(selectedPathways, pathways[,.N,Pathway][,Pathway])
if(length(pws_missing) > 0){
warning("Could not find the following pathways: ", paste(pws_missing, collapse = ", "), immediate. = TRUE)
}
pathways = pathways[Pathway %in% selectedPathways]
if(nrow(pathways) == 0){
stop("Zero pathways to plot!")
}
}
#Make sure pathway names are distinct from the gene names
dup_pw_names = intersect(pathways$Pathway, pathways$Gene)
if(length(dup_pw_names) > 0){
warning("Found following pathways with the same name as gene symbols. Renamed them with the suffix _pw\n", paste(dup_pw_names, collapse = ", "), immediate. = TRUE)
unique_pw_names = setdiff(pathways$Pathway, dup_pw_names)
pw_named_list = setNames(nm = c(unique_pw_names, dup_pw_names), object = c(unique_pw_names, paste0(dup_pw_names, "_pw")))
pathways$Pathway = unname(obj = pw_named_list[pathways$Pathway])
}
if(nrow(pathways[duplicated(Gene)]) > 0){
warning("Duplicated genes found across multiple pathways! This might cause issue while plotting.", immediate. = TRUE)
# pathways = pathways[!duplicated(Gene)]
}
genes = as.character(pathways$Gene)
om = createOncoMatrix(m = maf, g = genes, cbio = cBioPortal)
numMat = om$numericMatrix
mat_origin = om$oncoMatrix
#Check for any missing genes and ignore them if necessary
if(length(genes[!genes %in% rownames(numMat)]) > 0){
#warning('Missing following genes from provided pathways ', paste(genes[!genes %in% rownames(numMat)], collapse = ", "))
genes = genes[genes %in% rownames(numMat)]
}
}else if(!is.null(minMut)){
if(altered){
genes = getGeneSummary(x = maf)[order(AlteredSamples, decreasing = TRUE)][,.(Hugo_Symbol, AlteredSamples)]
}else{
genes = getGeneSummary(x = maf)[order(MutatedSamples, decreasing = TRUE)][,.(Hugo_Symbol, MutatedSamples)]
}
colnames(genes)[2] = "mutload"
genes[,fractMutated := mutload/totSamps]
if(minMut > 1){
genes = genes[mutload >= minMut, Hugo_Symbol]
}else{
genes = genes[fractMutated >= minMut, Hugo_Symbol]
}
om = createOncoMatrix(m = maf, g = genes, cbio = cBioPortal)
numMat = om$numericMatrix
mat_origin = om$oncoMatrix
}else { #If user does not provide gene list or MutSig results, draw TOP (default 20) genes
if(altered){
genes = getGeneSummary(x = maf)[order(AlteredSamples, decreasing = TRUE)][1:top, Hugo_Symbol]
}else{
genes = getGeneSummary(x = maf)[1:top, Hugo_Symbol]
}
om = createOncoMatrix(m = maf, g = genes, cbio = cBioPortal)
numMat = om$numericMatrix
mat_origin = om$oncoMatrix
}
#---remove genes from genesToIgnore if any
if(!is.null(genesToIgnore)){
numMat = numMat[!rownames(numMat) %in% genesToIgnore,]
mat_origin = mat_origin[!rownames(mat_origin) %in% genesToIgnore,]
}
tsbs = levels(getSampleSummary(x = maf)[,Tumor_Sample_Barcode])
if(!removeNonMutated){
tsb.include = matrix(data = 0, nrow = nrow(numMat),
ncol = length(tsbs[!tsbs %in% colnames(numMat)]))
colnames(tsb.include) = tsbs[!tsbs %in% colnames(numMat)]
rownames(tsb.include) = rownames(numMat)
numMat = cbind(numMat, tsb.include)
mat_origin = cbind(mat_origin, tsb.include)
}
#If user wannts to keep given gene order
if(keepGeneOrder){
if(fill){
numMat = numMat[genes, , drop = FALSE]
}else{
if(GeneOrderSort){
numMat = sortByGeneOrder(m = numMat, g = genes)
}else{
numMat = numMat[genes, , drop = FALSE]
}
}
mat = mat_origin[rownames(numMat), , drop = FALSE]
mat = mat[,colnames(numMat), drop = FALSE]
}
if(sortByMutation){
numMat_temp = sortByMutation(numMat = numMat, maf = maf)
numMat = numMat[rownames(numMat_temp), colnames(numMat_temp), drop = FALSE]
}
samp_sum = data.table::copy(getSampleSummary(x = maf))
samp_sum[,total := NULL]
if(cBioPortal){
samp_sum = data.table::melt(data = samp_sum, id.vars = "Tumor_Sample_Barcode")
colnames(samp_sum)[2] = "Variant_Classification"
vc = c("Nonstop_Mutation", "Frame_Shift_Del", "Missense_Mutation",
"Nonsense_Mutation", "Splice_Site", "Frame_Shift_Ins", "In_Frame_Del", "In_Frame_Ins")
vc.cbio = c("Truncating", "Truncating", "Missense", "Truncating", "Truncating", "Truncating",
"In-frame", "In-frame")
names(vc.cbio) = vc
samp_sum[,Variant_Classification_temp := vc.cbio[as.character(samp_sum$Variant_Classification)]]
samp_sum$Variant_Classification_temp = ifelse(test = is.na(samp_sum$Variant_Classification_temp), yes = as.character(samp_sum$Variant_Classification), no = samp_sum$Variant_Classification_temp)
samp_sum[,Variant_Classification := as.factor(as.character(Variant_Classification_temp))]
samp_sum[,Variant_Classification_temp := NULL]
samp_sum = samp_sum[,sum(value), .(Tumor_Sample_Barcode, Variant_Classification)]
samp_sum = data.table::dcast(data = samp_sum, formula = Tumor_Sample_Barcode~Variant_Classification, value.var = "V1")
}
if("CNV_total" %in% colnames(samp_sum)){
samp_sum[,CNV_total := NULL]
}
if(!includeColBarCN){
suppressWarnings(samp_sum[,Amp := NULL])
suppressWarnings(samp_sum[,Del := NULL])
}
data.table::setDF(x = samp_sum, rownames = as.character(samp_sum$Tumor_Sample_Barcode))
samp_sum = samp_sum[,-1, drop = FALSE]
#Parse annotations
if(!is.null(clinicalFeatures)){
if(is.null(annotationDat)){
annotation = parse_annotation_dat(annotationDat = maf, clinicalFeatures = clinicalFeatures)
}else{
annotation = parse_annotation_dat(annotationDat = annotationDat, clinicalFeatures = clinicalFeatures)
}
annotation = annotation[colnames(numMat),, drop = FALSE]
if(sortByAnnotation){
numMat = sortByAnnotation(numMat = numMat, maf = maf, anno = annotation, annoOrder = annotationOrder, group = groupAnnotationBySize, isNumeric = FALSE)
}
}
if(!is.null(sampleOrder)){
sampleOrder = as.character(sampleOrder)
sampleOrder = sampleOrder[sampleOrder %in% colnames(numMat)]
if(length(sampleOrder) == 0){
stop("None of the provided samples are present in the input MAF")
}
numMat = numMat[,sampleOrder, drop = FALSE]
}
gene_sum = apply(numMat, 1, function(x) length(x[x!=0]))
percent_alt = paste0(round(100*(apply(numMat, 1, function(x) length(x[x!=0]))/totSamps)), "%")
nonmut_samps = colnames(numMat)[!colnames(numMat) %in% rownames(samp_sum)]
if(length(nonmut_samps) > 0){
samp_sum_missing_samps = matrix(data = 0, nrow = length(nonmut_samps), ncol = ncol(samp_sum))
colnames(samp_sum_missing_samps) = colnames(samp_sum)
rownames(samp_sum_missing_samps) = nonmut_samps
samp_sum = rbind(samp_sum, samp_sum_missing_samps)
}
if(colbar_pathway){
samp_sum = getSampleSummary(subsetMaf(maf = maf, genes = genes))
samp_sum[,total := NULL]
if("CNV_total" %in% colnames(samp_sum)){
samp_sum[,CNV_total := NULL]
}
if(!includeColBarCN){
suppressWarnings(samp_sum[,Amp := NULL])
suppressWarnings(samp_sum[,Del := NULL])
}
data.table::setDF(x = samp_sum, rownames = as.character(samp_sum$Tumor_Sample_Barcode))
samp_sum = samp_sum[,-1]
top_bar_data = t(samp_sum[colnames(numMat),, drop = FALSE])
}else{
top_bar_data = t(samp_sum[colnames(numMat),, drop = FALSE])
}
#VC codes
if(cBioPortal){
if(is.null(colors)){
vc_col = grDevices::adjustcolor(col = c("black", "#33A02C", "brown"), alpha.f = 0.7)
vc_col = c('Truncating' = vc_col[1], 'Missense' = vc_col[2], 'In-frame' = vc_col[3], "Multi_Hit" = "#9b59b6", "pathway" = "#d35400")
}else{
vc_col = colors
if(!"pathway" %in% names(vc_col)){
vc_col = c(vc_col, "pathway" = "#535C68FF")
}
}
}else{
if(is.null(colors)){
vc_col = get_vcColors(websafe = FALSE)
}else{
vc_col = colors
if(!"pathway" %in% names(vc_col)){
vc_col = c(vc_col, "pathway" = "#535C68FF")
}
}
}
vc_codes = update_vc_codes(om_op = om)
vc_col = update_colors(x = vc_codes, y = vc_col)
if(nrow(numMat) == 1){
stop("Oncoplot requires at-least two genes for plottng.")
}
if(!is.null(pathways)){
vc_codes = c(vc_codes, '99' = 'pathway')
#Order genes in matrix by most frequently altered gene, followed by pathway size
temp_dat = data.table::data.table(Gene = rownames(numMat), percent_alt)
temp_dat = merge(temp_dat, pathways, all.x = TRUE)
temp_dat[is.na(temp_dat)] = "Unknown"
temp_dat$pct_alt = as.numeric(gsub(pattern = "%", replacement = "", x = temp_dat$percent_alt))
temp_dat = split(temp_dat, as.factor(as.character(temp_dat$Pathway)))
temp_dat_freq = lapply(temp_dat, function(x){
data.table::data.table(n = nrow(x), max_alt = max(x$pct_alt))
})
temp_dat_freq = data.table::rbindlist(l = temp_dat_freq, use.names = TRUE, fill = TRUE, idcol = "Pathway")
temp_dat_ord = temp_dat_freq[order(max_alt, n, decreasing = TRUE)][, Pathway]
temp_dat = lapply(temp_dat, function(x){x[order(pct_alt, decreasing = TRUE)]})[temp_dat_ord]
if("Unknown" %in% names(temp_dat)){
temp_dat_ord = c(grep(pattern = "Unknown", x = names(temp_dat), invert = TRUE, value = TRUE), "Unknown")
temp_dat = temp_dat[temp_dat_ord]
}
if(!is.null(path_order)){
path_order = intersect(x = path_order, names(temp_dat) )
temp_dat = temp_dat[path_order]
if(length(temp_dat) == 0){
stop("None of the pathways match the provided order!")
}
}
nm = lapply(seq_along(temp_dat), function(i){
x = numMat[temp_dat[[i]]$Gene,, drop = FALSE]
x = rbind(x, apply(x, 2, function(y) ifelse(test = sum(y) == 0, yes = 0, no = 99)))
if(names(temp_dat)[i] %in% rownames(x)){
stop("pathway name ", names(temp_dat)[i], " is same as one of its member genes! Please rename it.\ne.g, ", names(temp_dat)[i], " to ", names(temp_dat)[i], "_pw")
}
rownames(x)[nrow(x)] = names(temp_dat)[i]
x
})
numMat = do.call(rbind, nm)
#Add pathway information to the character matrix
mat_origin_path = rownames(numMat)[!rownames(numMat) %in% rownames(mat_origin)]
mat_origin_path = numMat[mat_origin_path,, drop = FALSE]
mat_origin_path[mat_origin_path == 0] = ""
mat_origin_path[mat_origin_path == "99"] = "pathway"
mat_origin_path = mat_origin_path[,colnames(mat_origin), drop = FALSE]
if(collapsePathway){
mat_origin = mat_origin_path
numMat = numMat[rownames(mat_origin),, drop = FALSE]
row_ord = names(sort(apply(numMat, 1, function(x) length(x[x == 0])), decreasing = FALSE))
numMat = numMat[row_ord,, drop = FALSE]
tnumMat = t(numMat) #transposematrix
numMat = t(tnumMat[do.call(order, c(as.list(as.data.frame(tnumMat)), decreasing = TRUE)), ]) #sort
if(sortByAnnotation){
numMat = sortByAnnotation(numMat = numMat, maf = maf, anno = annotation, annoOrder = annotationOrder, group = groupAnnotationBySize, isNumeric = FALSE)
}
}else{
mat_origin = rbind(mat_origin, mat_origin_path)
}
mat_origin = mat_origin[rownames(numMat), colnames(numMat), drop = FALSE]
}
#Plot layout
plot_layout(clinicalFeatures = clinicalFeatures, drawRowBar = drawRowBar, anno_height = anno_height,
drawColBar = drawColBar, draw_titv = draw_titv, exprsTbl = leftBarData, legend_height = legend_height)
#01: Draw scale axis for left barplot table
leftBarTitle = NULL
if(!is.null(leftBarData)){
leftBarTitle = colnames(leftBarData)[2]
colnames(leftBarData) = c('genes', 'exprn')
data.table::setDF(x = leftBarData, rownames = as.character(leftBarData$genes))
missing_samps = rownames(numMat)[!rownames(numMat) %in% rownames(leftBarData)]
if(length(missing_samps) > 0){
temp_data = data.frame(
row.names = missing_samps,
genes = missing_samps,
exprn = rep(0, length(missing_samps)),
stringsAsFactors = FALSE
)
leftBarData = rbind(leftBarData, temp_data)
}
leftBarData = leftBarData[rownames(numMat),, drop = FALSE]
#leftBarData = leftBarData[genes,, drop = FALSE]
if(is.null(leftBarLims)){
exprs_bar_lims = round(range(leftBarData$exprn, na.rm = TRUE), digits = 2)
}else{
exprs_bar_lims = leftBarLims
}
if(drawColBar){
par(mar = c(0.25 , 0, 1, 2), xpd = TRUE)
plot(x = NA, y = NA, type = "n", axes = FALSE,
xlim = exprs_bar_lims, ylim = c(0, 1), xaxs = "i")
}
}
#02: Draw top bar plot
if(drawColBar & is.null(topBarData)){
top_bar_data = top_bar_data[,colnames(numMat), drop = FALSE]
if(drawRowBar){
par(mar = c(0.25 , gene_mar, 2, 3), xpd = TRUE)
}else{
par(mar = c(0.25 , gene_mar, 2, 5), xpd = TRUE)
}
if(logColBar){
top_bar_data = apply(top_bar_data, 2, function(x) {
x_fract = x / sum(x)
x_log_total = log10(sum(x))
x_fract * x_log_total
})
top_bar_data[is.infinite(top_bar_data)] = 0
}
if(is.null(topBarLims)){
topBarLims = round(c(0, max(colSums(x = top_bar_data, na.rm = TRUE))), digits = 2)
}
plot(x = NA, y = NA, type = "n", xlim = c(0,ncol(top_bar_data)), ylim = topBarLims,
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i")
axis(side = 2, at = c(0, max(topBarLims, na.rm = TRUE)), las = 2, line = 0.5)
for(i in 1:ncol(top_bar_data)){
x = top_bar_data[,i]
names(x) = rownames(top_bar_data)
x = x[!x == 0]
if(length(x) > 0){
rect(xleft = i-1, ybottom = c(0, cumsum(x)[1:(length(x)-1)]), xright = i-0.1,
ytop = cumsum(x), col = vc_col[names(x)], border = NA, lwd = 0, xpd = FALSE)
}
}
if(logColBar){
mtext(text = "(log10)", side = 2, line = 2, cex = 0.6)
}else{
mtext(text = "TMB", side = 2, line = 2, cex = 0.6)
}
if(!is.null(topBarHline)){
abline(h = topBarHline, lty = 2, col = topBarHlineCol, xpd = FALSE)
mtext(text = topBarHline, side = 2, las = 2, cex = 0.6, line = 0.7)
}
}else if(!is.null(topBarData) & drawColBar){
# Draw extra clinical data in top
if(drawRowBar){
par(mar = c(0.25 , gene_mar, 2, 3), xpd = TRUE)
}else{
par(mar = c(0.25 , gene_mar, 2, 5), xpd = TRUE)
}
if(is.data.frame(topBarData)){
extdata = data.table::copy(x = topBarData)
colnames(extdata)[1] = "Tumor_Sample_Barcode"
rownames(extdata) <- extdata$Tumor_Sample_Barcode
}else{
message("Top barplot data is derived from clinical column: ", topBarData)
extdata <- data.frame(maf@clinical.data)[,c("Tumor_Sample_Barcode", topBarData)]
rownames(extdata) <- extdata$Tumor_Sample_Barcode
}
topbar_title = colnames(extdata)[2]
colnames(extdata)[2] = "barheights"
extdata[, "barheights"] <- as.numeric(as.character(extdata[, "barheights"]))
if(!all(colnames(top_bar_data) %in% rownames(extdata))){
missing_topbars = colnames(top_bar_data)[!colnames(top_bar_data) %in% rownames(extdata)]
warning("Missing column barplot data for ", length(missing_topbars), " samples")
extdata = rbind(extdata, data.frame(row.names = missing_topbars, Tumor_Sample_Barcode = missing_topbars, barheights = 0,stringsAsFactors = FALSE))
}
extdata <- extdata[colnames(top_bar_data),]
if(is.null(topBarLims)){
topBarLims = round(range(extdata[,"barheights"], na.rm = TRUE), digits = 2)
}
plot(x = NA, y = NA, type = "n", xlim = c(0,nrow(extdata)), ylim = topBarLims,
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i")
axis(side = 2, at = topBarLims, las = 2, line = 0.5)
for(i in 1:nrow(extdata)){
graphics::rect(xleft = i-1, xright = i-0.1, ybottom = 0, ytop = as.numeric(extdata[i, "barheights"]),
col = "#535c68", border = NA, lwd = 0, xpd = FALSE)
}
title(ylab = topbar_title, font = 1, cex.lab = legendFontSize, xpd = TRUE)
if(!is.null(topBarHline)){
abline(h = topBarHline, lty = 2, col = topBarHlineCol)
}
}
#03: Draw scale for right barplot
rightBarTitle = NULL
if(drawRowBar){
if(is.null(rightBarData)){
side_bar_lims = c(0, max(unlist(apply(numMat, 1, function(x) cumsum(table(x[x!=0])))), na.rm = TRUE))
rightBarTitle = "No. of samples"
}else{
rightBarTitle = colnames(rightBarData)[2]
colnames(rightBarData) = c('genes', 'exprn')
data.table::setDF(x = rightBarData, rownames = as.character(rightBarData$genes))
missing_samps = rownames(numMat)[!rownames(numMat) %in% rownames(rightBarData)]
if(length(missing_samps) > 0){
temp_data = data.frame(
row.names = missing_samps,
genes = missing_samps,
exprn = rep(0, length(missing_samps)),
stringsAsFactors = FALSE
)
rightBarData = rbind(rightBarData, temp_data)
}
rightBarData = rightBarData[rownames(numMat),, drop = FALSE]
if(is.null(rightBarLims)){
side_bar_lims = round(range(rightBarData$exprn, na.rm = TRUE), digits = 2)
}else{
side_bar_lims = rightBarLims
}
}
if(drawColBar){
par(mar = c(0.25 , 0, 1, 2), xpd = TRUE)
plot(x = NA, y = NA, type = "n", axes = FALSE,
xlim = side_bar_lims, ylim = c(0, 1), xaxs = "i")
}
}
#Draw expression barplot
if(!is.null(leftBarData)){
if(showTumorSampleBarcodes){
if(!drawRowBar & !drawColBar){
par(mar = c(barcode_mar, 1, 2.5, 0), xpd = TRUE)
}else if(!drawRowBar & drawColBar){
par(mar = c(barcode_mar, 1, 0, 5), xpd = TRUE)
}else if(drawRowBar & !drawColBar){
par(mar = c(barcode_mar, 1, 2.5, 0), xpd = TRUE)
} else{
par(mar = c(barcode_mar, 1, 0, 0), xpd = TRUE)
}
}else{
if(!drawRowBar & !drawColBar){
par(mar = c(0.5, 1, 2.5, 0), xpd = TRUE)
}else if(!drawRowBar & drawColBar){
par(mar = c(0.5, 1, 0, 5), xpd = TRUE)
}else if(drawRowBar & !drawColBar){
par(mar = c(0.5, 1, 2.5, 0), xpd = TRUE)
} else{
par(mar = c(0.5, 1, 0, 0), xpd = TRUE)
}
}
plot(x = NA, y = NA, type = "n", xlim = rev(-exprs_bar_lims), ylim = c(0, nrow(leftBarData)),
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i", yaxs = "i") #
for(i in 1:nrow(leftBarData)){
x = rev(leftBarData$exprn)[i]
rect(ybottom = i-1, xleft = -x, ytop = i-0.1,
xright = 0, col = "#535c68", border = NA, lwd = 0)
}
axis(side = 3, at = rev(-exprs_bar_lims), outer = FALSE, line = 0.25, labels = rev(exprs_bar_lims))
mtext(text = leftBarTitle, side = 3, line = 0.50, cex = 0.6)
if(!is.null(leftBarVline)){
abline(h = leftBarVline, lty = 2, col = leftBarVlineCol, xpd = FALSE)
}
}
#04: Draw the main matrix
if(showTumorSampleBarcodes){
if(!drawRowBar & !drawColBar){
par(mar = c(barcode_mar, gene_mar, 2.5, 5), xpd = TRUE)
}else if(!drawRowBar & drawColBar){
par(mar = c(barcode_mar, gene_mar, 0, 5), xpd = TRUE)
}else if(drawRowBar & !drawColBar){
par(mar = c(barcode_mar, gene_mar, 2.5, 3), xpd = TRUE)
} else{
par(mar = c(barcode_mar, gene_mar, 0, 3), xpd = TRUE)
}
}else{
if(!drawRowBar & !drawColBar){
par(mar = c(0.5, gene_mar, 2.5, 5), xpd = TRUE)
}else if(!drawRowBar & drawColBar){
par(mar = c(0.5, gene_mar, 0, 5), xpd = TRUE)
}else if(drawRowBar & !drawColBar){
par(mar = c(0.5, gene_mar, 2.5, 3), xpd = TRUE)
} else{
par(mar = c(0.5, gene_mar, 0, 3), xpd = TRUE)
}
}
nm = t(apply(numMat, 2, rev))
nm[nm == 0] = NA
image(x = 1:nrow(nm), y = 1:ncol(nm), z = nm, axes = FALSE, xaxt="n", yaxt="n",
xlab="", ylab="", col = "white") #col = "#FC8D62"
#Plot for all variant classifications
vc_codes_temp = vc_codes[!vc_codes %in% om$cnvc]
for(i in 2:length(names(vc_codes_temp))){
vc_code = vc_codes_temp[i]
col = vc_col[vc_code]
nm = t(apply(numMat, 2, rev))
#print(names(vc_code))
nm[nm != names(vc_code)] = NA
#print(paste0(vc_code, " ", col))
#Suppress warning due to min/max applied to a vector of NAs; Issue: #286
#This is an harmless warning as matrix is looped over all VC's and missing VC's form NA's (which are plotted in gray)
suppressWarnings(image(x = 1:nrow(nm), y = 1:ncol(nm), z = nm, axes = FALSE, xaxt="n", yaxt="n",
xlab="", ylab="", col = col, add = TRUE))
}
#Add blanks
nm = t(apply(numMat, 2, rev))
nm[nm != 0] = NA
image(x = 1:nrow(nm), y = 1:ncol(nm), z = nm, axes = FALSE, xaxt="n", yaxt="n", xlab="", ylab="", col = bgCol, add = TRUE)
#Add CNVs if any
mat_origin = mat_origin[rownames(numMat), colnames(numMat), drop = FALSE]
if(writeMatrix){
write.table(mat_origin, "onco_matrix.txt", sep = "\t", quote = FALSE)
}
mo = t(apply(mat_origin, 2, rev))
##Complex events (mutated as well as CN altered)
complex_events = unique(grep(pattern = ";", x = mo, value = TRUE))
if(length(complex_events) > 0){
for(i in 1:length(complex_events)){
ce = complex_events[i]
#mo = t(apply(mat_origin, 2, rev))
ce_idx = which(mo == ce, arr.ind = TRUE)
ce = unlist(strsplit(x = ce, split = ";", fixed = TRUE))
nm_temp = matrix(NA, nrow = nrow(nm), ncol = ncol(nm))
nm_temp[ce_idx] = 0
image(x = 1:nrow(nm_temp), y = 1:ncol(nm_temp), z = nm_temp, axes = FALSE, xaxt="n",
yaxt="n", xlab="", ylab="", col = vc_col[ce[2]], add = TRUE)
ce_idx = which(t(nm_temp) == 0, arr.ind = TRUE)
for(i in seq_len(nrow(ce_idx))){
rowi = ce_idx[i,1]
coli = ce_idx[i,2]
rect(xleft = coli-0.5, ybottom = rowi-0.25, xright = coli+0.5, ytop = rowi+0.25, col = vc_col[ce[1]], border = NA, lwd = 0)
}
}
}
for(cnevent in om$cnvc){
cn_idx = which(mo == cnevent, arr.ind = TRUE)
if(nrow(cn_idx) > 0){
nm_temp = matrix(NA, nrow = nrow(nm), ncol = ncol(nm))
nm_temp[cn_idx] = 0
image(x = 1:nrow(nm_temp), y = 1:ncol(nm_temp), z = nm_temp, axes = FALSE, xaxt="n",
yaxt="n", xlab="", ylab="", col = bgCol, add = TRUE)
cn_idx = which(t(nm_temp) == 0, arr.ind = TRUE)
for(i in seq_len(nrow(cn_idx))){
rowi = cn_idx[i,1]
coli = cn_idx[i,2]
rect(xleft = coli-0.5, ybottom = rowi-0.25, xright = coli+0.5, ytop = rowi+0.25, col = vc_col[cnevent], border = NA, lwd = 0)
}
}
}
#Add grids
abline(h = (0:ncol(nm)) + 0.5, col = borderCol, lwd = sepwd_genes, xpd = FALSE)
abline(v = (0:nrow(nm)) + 0.5, col = borderCol, lwd = sepwd_samples, xpd = FALSE)
#Add boxes if pathways are opted
if(!is.null(pathways)){
temp_dat = lapply(seq_along(temp_dat), function(i){
data.table::rbindlist(list(temp_dat[[i]], data.table::data.table(Gene = names(temp_dat)[i])), use.names = TRUE, fill = TRUE)
})
temp_dat = data.table::rbindlist(l = temp_dat, use.names = TRUE, fill = TRUE)
#return(temp_dat)
temp_dat[, row_id := nrow(temp_dat):1]
temp_dat = split(temp_dat, as.factor(temp_dat$Pathway))
if(!collapsePathway){
lapply(temp_dat, function(td){
rect(xleft = 0.5, ybottom = min(td$row_id)-0.499, xright = nrow(nm)+0.5, ytop = max(td$row_id)+0.499, border = pwLineCol, lwd = pwLineWd)
})
}
#New percent alt with pathways
percent_alt = rev(paste0(apply(nm, 2, function(x){
round(length(x[is.na(x)])/totSamps * 100)
}), "%"))
}
#Draw if any additional features are requested
additionalFeature_legend = FALSE
if(!is.null(additionalFeature)){
#If its a list, multi features to be mapped
if(is.list(additionalFeature)){
for(af_idx in seq_along(additionalFeature)){
af = additionalFeature[[af_idx]]
if(length(af) < 2){
stop("additionalFeature must be of length two. See ?oncoplot for details.")
}
af_dat = subsetMaf(maf = maf, genes = rownames(numMat), tsb = colnames(numMat), fields = af[1], includeSyn = FALSE, mafObj = FALSE)
if(length(which(colnames(af_dat) == af[1])) == 0){
message(paste0("Column ", af[1], " not found in maf. Here are available fields.."))
print(getFields(maf))
stop()
}
colnames(af_dat)[which(colnames(af_dat) == af[1])] = 'temp_af'
af_dat = af_dat[temp_af %in% af[2]]
if(nrow(af_dat) == 0){
warning(paste0("No samples are enriched for ", af[2], " in ", af[1]))
}else{
af_mat = data.table::dcast(data = af_dat, Tumor_Sample_Barcode ~ Hugo_Symbol, value.var = "temp_af", fun.aggregate = length)
af_mat = as.matrix(af_mat, rownames = "Tumor_Sample_Barcode")
nm = t(apply(numMat, 2, rev))
if(length(additionalFeaturePch) != length(additionalFeature)){
warning("Provided pch for additional features are recycled")
additionalFeaturePch = rep(additionalFeaturePch, length(additionalFeature))
}
if(length(additionalFeatureCol) != length(additionalFeature)){
warning("Provided colors for additional features are recycled")
additionalFeatureCol = rep(additionalFeatureCol, length(additionalFeature))
}
lapply(seq_len(nrow(af_mat)), function(i){
af_i = af_mat[i,, drop = FALSE]
af_i_genes = colnames(af_i)[which(af_i > 0)]
af_i_sample = rownames(af_i)
lapply(af_i_genes, function(ig){
af_i_mat = matrix(c(which(rownames(nm) == af_i_sample),
which(colnames(nm) == ig)),
nrow = 1)
points(af_i_mat, pch = additionalFeaturePch[af_idx], col= additionalFeatureCol[af_idx], cex = additionalFeatureCex)
})
})
additionalFeature_legend = TRUE
}
}
}else{
if(length(additionalFeature) < 2){
stop("additionalFeature must be of length two. See ?oncoplot for details.")
}
af_dat = subsetMaf(maf = maf, genes = rownames(numMat), tsb = colnames(numMat), fields = additionalFeature[1], includeSyn = FALSE, mafObj = FALSE)
if(length(which(colnames(af_dat) == additionalFeature[1])) == 0){
message(paste0("Column ", additionalFeature[1], " not found in maf. Here are available fields.."))
print(getFields(maf))
stop()
}
colnames(af_dat)[which(colnames(af_dat) == additionalFeature[1])] = 'temp_af'
af_dat = af_dat[temp_af %in% additionalFeature[2]]
if(nrow(af_dat) == 0){
warning(paste0("No samples are enriched for ", additionalFeature[2], " in ", additionalFeature[1]))
}else{
af_mat = data.table::dcast(data = af_dat, Tumor_Sample_Barcode ~ Hugo_Symbol, value.var = "temp_af", fun.aggregate = length)
af_mat = as.matrix(af_mat, rownames = "Tumor_Sample_Barcode")
nm = t(apply(numMat, 2, rev))
lapply(seq_len(nrow(af_mat)), function(i){
af_i = af_mat[i,, drop = FALSE]
af_i_genes = colnames(af_i)[which(af_i > 0)]
af_i_sample = rownames(af_i)
lapply(af_i_genes, function(ig){
af_i_mat = matrix(c(which(rownames(nm) == af_i_sample),
which(colnames(nm) == ig)),
nrow = 1)
points(af_i_mat, pch = additionalFeaturePch, col= additionalFeatureCol, cex = additionalFeatureCex)
})
})
additionalFeature_legend = TRUE
}
}
}
# #Draw grids for the samples
# if(!is.null(gridFeature)){
# To add: https://github.com/PoisonAlien/maftools/issues/528
# }
#Add box border
if (drawBox) {
box(lty = 'solid', col = '#535c68', lwd = 2)
}
mtext(text = colnames(nm), side = 2, at = 1:ncol(nm),
font = 3, line = 0.4, cex = fontSize, las = 2)
if(showPct){
mtext(text = rev(percent_alt), side = 4, at = 1:ncol(nm),
font = 1, line = 0.4, cex = fontSize, las = 2, adj = 0.15)
}
if(showTumorSampleBarcodes){
if(!is.null(tsbToPIDs)){
if(is.data.frame(x = tsbToPIDs)){
colnames(x = tsbToPIDs) = c("Tumor_Sample_Barcode", "PID")
rownames(tsbToPIDs) = tsbToPIDs$Tumor_Sample_Barcode
}else{
clin = getClinicalData(x = maf)
if(!tsbToPIDs %in% colnames(clin)){
stop(tsbToPIDs, " does not exist in the clinical data!")
}
data.table::setDF(x = clin, rownames = clin$Tumor_Sample_Barcode)
tsbToPIDs = clin[,c("Tumor_Sample_Barcode", tsbToPIDs)]
colnames(x = tsbToPIDs) = c("Tumor_Sample_Barcode", "PID")
}
if(all(rownames(nm) %in% rownames(tsbToPIDs))){
text(x =1:nrow(nm), y = par("usr")[3] - 0.2,
labels = tsbToPIDs[rownames(nm), 2], srt = barcodeSrt, font = 1, cex = SampleNamefontSize, adj = 1)
}else{
missing_tsb = setdiff(rownames(nm), rownames(tsbToPIDs))
warning("Follwing Tumor_Sample_Barcodes are missing corresponding PIDs: ", paste(missing_tsb, collapse = ", "))
text(x =1:nrow(nm), y = par("usr")[3] - 0.2,
labels = tsbToPIDs[rownames(nm), 2], srt = barcodeSrt, font = 1, cex = SampleNamefontSize, adj = 1)
}
}else{
text(x =1:nrow(nm), y = par("usr")[3] - 0.2,
labels = rownames(nm), srt = barcodeSrt, font = 1, cex = SampleNamefontSize, adj = 1)
}
}
#05: Draw right side barplot
if(drawRowBar){
if(showTumorSampleBarcodes){
if(!drawRowBar || !drawColBar){
par(mar = c(barcode_mar, 0, 2.5, 1), xpd = TRUE)
}else{
par(mar = c(barcode_mar, 0, 0, 1), xpd = TRUE)
}
}else{
if(!drawRowBar || !drawColBar){
par(mar = c(0.5, 0, 2.5, 1), xpd = TRUE)
}else{
par(mar = c(0.5, 0, 0, 1), xpd = TRUE)
}
}
if(is.null(rightBarData)){
#side_bar_data = apply(numMat, 1, function(x) table(x[x!=0]))
side_bar_data = lapply(seq_len(nrow(numMat)), function(i) {
xi = numMat[i, ]
table(xi[!xi == 0])
})
names(side_bar_data) = rownames(numMat)
plot(x = NA, y = NA, type = "n", xlim = side_bar_lims, ylim = c(0, length(side_bar_data)),
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i", yaxs = "i") #
for(i in 1:length(side_bar_data)){
x = rev(side_bar_data)[[i]]
if(length(x) > 0){
rect(ybottom = i-1, xleft = c(0, cumsum(x)[1:(length(x)-1)]), ytop = i-0.1,
xright = cumsum(x), col = vc_col[vc_codes[names(x)]], border = NA, lwd = 0)
}
}
if(!is.null(rightBarVline)){
abline(v = rightBarVline, lty = 2, col = rightBarVlineCol, xpd = FALSE)
}
axis(side = 3, at = side_bar_lims, outer = FALSE, line = 0.25)
mtext(text = rightBarTitle, side = 3, line = 0.5, cex = 0.6)
}else{
plot(x = NA, y = NA, type = "n", xlim = side_bar_lims, ylim = c(0, nrow(rightBarData)),
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i", yaxs = "i") #
for(i in 1:nrow(rightBarData)){
x = rev(rightBarData$exprn)[i]
rect(ybottom = i-1, xleft = x, ytop = i-0.1,
xright = 0, col = "#535c68", border = NA, lwd = 0)
}
if(!is.null(rightBarVline)){
abline(v = rightBarVline, lty = 2, col = rightBarVlineCol, xpd = FALSE)
}
axis(side = 3, at = side_bar_lims, outer = FALSE, line = 0.25)
mtext(text = rightBarTitle, side = 3, line = 0.5, cex = 0.6)
}
}
#06: Plot annotations if any
if(!is.null(clinicalFeatures)){
annotation = annotation[colnames(numMat), ncol(annotation):1, drop = FALSE]
if(is.null(annotationColor)){
#In case no-color codes are provided
annotationColor = get_anno_cols(ann = annotation)
}else{
#If provided - generate a continuous palette for continuous variables
annotationColor = annotationColor[colnames(annotation)]
annotationColor_temp = lapply(seq_along(annotationColor), function(idx){
if(is.numeric(annotation[,idx])){
#print(annotationColor[[idx]])
seq_col_pal = c("Blues", "BuGn", "BuPu", "GnBu", "Greens", "Greys", "Oranges",
"OrRd", "PuBu", "PuBuGn", "PuRd", "Purples", "RdPu", "Reds",
"YlGn", "YlGnBu", "YlOrBr", "YlOrRd")
col_match = annotationColor[idx] %in% seq_col_pal
if(!col_match){
stop("numeric annotation color for ", names(annotationColor)[idx] , " must be a sequential color palette!\n", paste(seq_col_pal, collapse = " "))
}
x = annotation[,idx]
x_unique = unique(x)
numericAnnoCol = RColorBrewer::brewer.pal(n = 9, name = annotationColor[[idx]])
ann_lvls_cols = colorRampPalette(numericAnnoCol)(length(x_unique))
names(ann_lvls_cols) = x_unique[order(x_unique, na.last = TRUE)]
ann_lvls_cols = ann_lvls_cols[as.character(x)]
}else{
ann_lvls_cols = annotationColor[[idx]]
ann_lvls_cols = ann_lvls_cols[as.character(annotation[,idx])]
}
ann_lvls_cols
})
names(annotationColor_temp) = names(annotationColor)
#dirty way to check which columns are of numeric
temp_col_class = c()
for(i in 1:ncol(annotation)){
if(is.numeric(annotation[,i])){
temp_col_class = c(temp_col_class, "numeric")
}else{
temp_col_class = c(temp_col_class, "char")
}
}
annotationColor = annotationColor_temp
attr(annotationColor, "class") = temp_col_class
}
ann_classes = attributes(annotationColor)$class
#return(list(annotation, annotationColor))
if(!is.null(leftBarData)){
plot.new()
}
if(!drawRowBar){
par(mar = c(0, gene_mar, 0, 5), xpd = TRUE)
}else{
par(mar = c(0, gene_mar, 0, 3), xpd = TRUE)
}
image(x = 1:nrow(annotation), y = 1:ncol(annotation), z = matrix(data = 0, nrow = nrow(annotation), ncol = ncol(annotation)),
axes = FALSE, xaxt="n", yaxt="n", bty = "n",
xlab="", ylab="", col = "gray90") #col = "#FC8D62"
for(i in 1:ncol(annotation)){
temp_anno = annotation
if(ncol(temp_anno) > 1){
rm_idx = setdiff(y = i, x = 1:ncol(temp_anno))
temp_anno[,rm_idx] = NA
}
if(ann_classes[i] == "numeric"){
temp_anno = apply(temp_anno, 2, as.numeric)
suppressWarnings(image(x = 1:nrow(temp_anno), y = 1:ncol(temp_anno), z = temp_anno,
axes = FALSE, xaxt="n", yaxt="n", xlab="", ylab="", col = annotationColor[[i]], add = TRUE))
}else{
temp_lvls = unique(names(annotationColor[[i]])) #unique(temp_anno[,i])
temp_lvls = temp_lvls[complete.cases(temp_lvls)] #Remove NAs
temp_lvls_cols = annotationColor[[i]][temp_lvls]
temp_anno[,i] = match(x = temp_anno[,i], table = unique(temp_lvls), nomatch = NA)
temp_anno = apply(temp_anno, 2, as.numeric)
suppressWarnings(image(x = 1:nrow(temp_anno), y = 1:ncol(temp_anno), z = temp_anno,
axes = FALSE, xaxt="n", yaxt="n", xlab="", ylab="", col = temp_lvls_cols, add = TRUE))
}
}
#Add grids
rect(xleft = 0.5, ybottom = (0:(ncol(annotation))) + 0.5, xright = nrow(annotation)+0.5, ytop = (0:(ncol(annotation))) + 0.5, border = annoBorderCol, lwd = 0.7)
abline(v = (0:(nrow(nm))) + 0.5, col = annoBorderCol, lwd = sepwd_samples, xpd = FALSE)
mtext(text = colnames(annotation), side = 4,
font = 1, line = 0.4, cex = fontSize, las = 2, at = 1:ncol(annotation))
if(drawRowBar){
plot.new()
}
}
#07: Draw TiTv plot
if(draw_titv){
titv_dat = titv(maf = maf, useSyn = TRUE, plot = FALSE)
titv_dat = titv_dat$fraction.contribution
data.table::setDF(x = titv_dat, rownames = as.character(titv_dat$Tumor_Sample_Barcode))
titv_dat = titv_dat[,-1]
titv_dat = t(titv_dat[colnames(numMat), ])
missing_samps = colnames(numMat)[!colnames(numMat) %in% colnames(titv_dat)]
if(length(missing_samps) > 0){
temp_data = matrix(data = 0, nrow = 6, ncol = length(missing_samps))
colnames(temp_data) = missing_samps
titv_dat = cbind(titv_dat, temp_data)
titv_dat = titv_dat[,colnames(numMat)]
}
if(!is.null(leftBarData)){
plot.new()
}
if(!drawRowBar){
par(mar = c(0, gene_mar, 0, 5), xpd = TRUE)
}else{
par(mar = c(0, gene_mar, 0, 3), xpd = TRUE)
}
plot(x = NA, y = NA, type = "n", xlim = c(0,ncol(titv_dat)), ylim = c(0, 100),
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i")
if(is.null(titv_col)){
titv_col = get_titvCol(alpha = 1)
}
for(i in 1:ncol(titv_dat)){
x = titv_dat[,i]
x = x[x > 0]
if(length(x) > 0){
rect(xleft = i-1, ybottom = c(0, cumsum(x)[1:(length(x)-1)]), xright = i-0.1,
ytop = cumsum(x), col = titv_col[names(x)], border = NA, lwd = 0)
}else{
rect(xleft = i-1, ybottom = c(0, 100), xright = i-0.1,
ytop = 100, col = "gray70", border = NA, lwd = 0)
}
}
if(drawRowBar){
par(mar = c(0, 0, 1, 6), xpd = TRUE)
plot(NULL,ylab='',xlab='', xlim=0:1, ylim=0:1, axes = FALSE)
lep = legend("topleft", legend = names(titv_col),
col = titv_col, border = NA, bty = "n",
ncol= 2, pch = 15, xpd = TRUE, xjust = 0, yjust = 0,
cex = legendFontSize)
}else{
mtext(text = names(titv_col)[1:3], side = 2, at = c(25, 50, 75),
font = 1, line = 0.4, cex = fontSize, las = 2, col = titv_col[1:3])
mtext(text = names(titv_col)[4:6], side = 4, at = c(25, 50, 75),
font = 1, line = 0.4, cex = fontSize, las = 2, col = titv_col[4:6], adj = 0.15)
}
}
#08: Add legends
par(mar = c(0, 0.5, 0, 0), xpd = TRUE)
plot(NULL,ylab='',xlab='', xlim=0:1, ylim=0:1, axes = FALSE)
leg_classes = vc_col[vc_codes[2:length(vc_codes)]]
leg_classes_pch = rep(15, length(leg_classes))
if(additionalFeature_legend){
if(!is(object = additionalFeature, class2 = "list")){
if(length(additionalFeature) < 2){
stop("additionalFeature must be of length two. See ?oncoplot for details.")
}else{
additionalFeature = list(additionalFeature)
}
}
if(length(additionalFeaturePch) != length(additionalFeature)){
additionalFeaturePch = rep(additionalFeaturePch, length(additionalFeature))
}
if(length(additionalFeatureCol) != length(additionalFeature)){
additionalFeatureCol = rep(additionalFeatureCol, length(additionalFeature))
}
for(af_idx in 1:length(additionalFeature)){
af = additionalFeature[[af_idx]]
leg_classes = c(leg_classes,additionalFeatureCol[af_idx])
names(leg_classes)[length(leg_classes)] = paste(af, collapse = ":")
leg_classes_pch = c(leg_classes_pch, additionalFeaturePch[af_idx])
}
}
if(collapsePathway){
leg_classes = leg_classes['pathway']
}
lep = legend("topleft", legend = names(leg_classes),
col = leg_classes, border = NA, bty = "n",
ncol= 2, pch = leg_classes_pch, xpd = TRUE, xjust = 0, yjust = 0, cex = legendFontSize)
x_axp = 0+lep$rect$w
if(!is.null(clinicalFeatures)){
#Check which of the clinicalFeatures are of numeric
ann_classes = attributes(annotationColor)$class
num_start = 0.1 #numeric legend starts from here
num_width = 0.1
num_incr = 0.2
for(i in seq_along(ann_classes)){
#x = unique(annotation[,i])
temp_lvls = unique(names(annotationColor[[i]]))
temp_lvls = temp_lvls[complete.cases(temp_lvls)]
temp_lvls_cols = annotationColor[[i]][temp_lvls]
x = annotationColor[[i]]
is_num = ann_classes[i] == "numeric"
xt = names(x)
if(is_num){
xt = as.numeric(xt)
xt = xt[!is.na(xt)]
xt = xt[!is.infinite(xt)]
xt = sort(xt)
x = x[as.character(xt)]
x_range = round(range(xt), 2)
image(
y = c(0.1, 0.2),
x = seq(num_start, num_start+num_width, length.out = 200),
z = t(matrix(seq(0,1,length.out = 200), nrow= 1)),
col = x, xlim = c(0, 1), ylim = c(0, 1), axes = FALSE, xlab = NA, ylab = NA, add = TRUE
)
text(x = c(num_start, num_start+num_width), y = 0.2, labels = x_range, pos = 3)
text(x = num_start+(num_width/2), y = 0.1, labels = names(annotationColor)[i], pos = 1, adj = 1)
num_start = num_start + num_incr
}else{
if(length(temp_lvls_cols) <= 4){
n_col = 1
}else{
n_col = (length(temp_lvls_cols) %/% 4)+1
}
lep = legend(x = x_axp, y = 1, legend = names(temp_lvls_cols),
col = temp_lvls_cols, border = NA,
ncol= n_col, pch = 15, xpd = TRUE, xjust = 0, bty = "n",
cex = annotationFontSize, title = names(annotationColor)[i],
title.adj = 0)
x_axp = x_axp + lep$rect$w
}
}
}
n_mut_samps = length(which(colSums(numMat) != 0))
altStat = paste0("Altered in ", n_mut_samps, " (", round(n_mut_samps/totSamps, digits = 4)*100, "%) of ", totSamps, " samples.")
if(showTitle){
if(is.null(titleText)){
title(main = altStat, outer = TRUE, line = -1, cex.main = titleFontSize)
}else{
title(main = titleText, outer = TRUE, line = -1, cex.main = titleFontSize)
}
}
return(invisible(list(oncomatrix = numMat, vc_legend = vc_codes, vc_color = leg_classes)))
}
PoisonAlien/maftools documentation built on Nov. 5, 2024, 4:12 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions oncoplot
Documented in oncoplot
#' draw an oncoplot
#' @description takes output generated by read.maf and draws an oncoplot
#'
#' @details
#' Takes an \code{\link{MAF}} object as an input and plots it as a matrix. Any desired clincal features can be added at the bottom of the oncoplot by providing \code{clinicalFeatures}.
#' Oncoplot can be sorted either by mutations or by clinicalFeatures using arguments \code{sortByMutation} and \code{sortByAnnotation} respectively.
#'
#'By setting `pathways` argument either `sigpw` or `smgbp` - cohort can be summarized by altered pathways. pathways argument also accepts a custom pathway list in the form of a two column tsv file or a data.frame containing gene names and their corresponding pathway.
#'
#' @references Bailey, Matthew H et al. “Comprehensive Characterization of Cancer Driver Genes and Mutations.” Cell vol. 173,2 (2018): 371-385.e18. doi:10.1016/j.cell.2018.02.060
#' Sanchez-Vega, Francisco et al. “Oncogenic Signaling Pathways in The Cancer Genome Atlas.” Cell vol. 173,2 (2018): 321-337.e10. doi:10.1016/j.cell.2018.03.035
#'
#' @param maf an \code{\link{MAF}} object generated by \code{\link{read.maf}}
#' @param top how many top genes to be drawn. defaults to 20.
#' @param minMut draw all genes with `min` number of mutations. Can be an integer or fraction (of samples mutated), Default NULL
#' @param genes Just draw oncoplot for these genes. Default NULL.
#' @param altered Default FALSE. Chooses top genes based on muatation status. If \code{TRUE} chooses top genes based alterations (CNV or mutation).
#' @param drawColBar logical plots top barplot for each sample. Default \code{TRUE}.
#' @param drawRowBar logical. Plots righ barplot for each gene. Default \code{TRUE}.
#' @param leftBarData Data for leftside barplot. Must be a data.frame with two columns containing gene names and values. Default `NULL`
#' @param leftBarLims limits for `leftBarData`. Default `NULL`.
#' @param leftBarVline Draw vertical lines at these values. Default `NULL`.
#' @param leftBarVlineCol Line color for `leftBarVline` Default gray70
#' @param topBarData Default `NULL` which draws absolute number of mutation load for each sample. Can be overridden by choosing one clinical indicator(Numeric) or by providing a two column data.frame containing sample names and values for each sample. This option is applicable when only `drawColBar` is TRUE.
#' @param topBarLims limits for `topBarData`. Default `NULL`.
#' @param topBarHline Draw horizontal lines at these values. Default `NULL`.
#' @param topBarHlineCol Line color for `topBarHline.` Default gray70
#' @param rightBarData Data for rightside barplot. Must be a data.frame with two columns containing to gene names and values. Default `NULL` which draws distribution by variant classification. This option is applicable when only `drawRowBar` is TRUE.
#' @param rightBarLims limits for `rightBarData`. Default `NULL`.
#' @param rightBarVline Draw vertical lines at these values. Default `NULL`.
#' @param rightBarVlineCol Line color for `rightBarVline` Default gray70
#' @param logColBar Plot top bar plot on log10 scale. Default \code{FALSE}.
#' @param includeColBarCN Whether to include CN in column bar plot. Default TRUE
#' @param clinicalFeatures columns names from `clinical.data` slot of \code{MAF} to be drawn in the plot. Default NULL.
#' @param annotationColor Custom colors to use for `clinicalFeatures`. Must be a named list containing a named vector of colors. Default NULL. See example for more info.
#' @param annotationDat If MAF file was read without clinical data, provide a custom \code{data.frame} with a column \code{Tumor_Sample_Barcode} containing sample names along with rest of columns with annotations.
#' You can specify which columns to be drawn using `clinicalFeatures` argument.
#' @param pathways Default `NULL`. Can be `sigpw`, `smgbp`, or a two column data.frame/tsv-file with genes and corresponding pathway mappings.`
#' @param topPathways Top most altered pathways to draw. Default 3. Mutually exclusive with `selectedPathways`
#' @param path_order Default `NULL` Manually specify the order of pathways
#' @param selectedPathways Manually provide the subset of pathway names to be selected from `pathways`. Default NULL. In case `pathways` is `auto` draws top 3 altered pathways.
#' @param collapsePathway Shows only rows corresponding to the pathways. Default FALSE.
#' @param pwLineCol Color for the box around the pathways Default #535c68
#' @param pwLineWd Line width for the box around the pathways Default Default 1
#' @param draw_titv logical Includes TiTv plot. \code{FALSE}
#' @param titv_col named vector of colors for each transition and transversion classes. Should be of length six with the names "C>T" "C>G" "C>A" "T>A" "T>C" "T>G". Default NULL.
#' @param showTumorSampleBarcodes logical to include sample names.
#' @param tsbToPIDs Custom names for Tumor_Sample_Barcodes. Can be a column name in clinicaldata or a 2 column data.frame of Tumor_Sample_Barcodes to patient ID mappings. Applicable only when `showTumorSampleBarcodes = TRUE`. Default NULL.
#' @param barcode_mar Margin width for sample names. Default 4
#' @param barcodeSrt Rotate sample labels. Default 90.
#' @param gene_mar Margin width for gene names. Default 5
#' @param anno_height Height of plotting area for sample annotations. Default 1
#' @param legend_height Height of plotting area for legend. Default 4
#' @param sortByAnnotation logical sort oncomatrix (samples) by provided `clinicalFeatures`. Sorts based on first `clinicalFeatures`. Defaults to FALSE. column-sort
#' @param groupAnnotationBySize Further group `sortByAnnotation` orders by their size. Defaults to TRUE. Largest groups comes first.
#' @param annotationOrder Manually specify order for annotations. Works only for first `clinicalFeatures`. Default NULL.
#' @param sortByMutation Force sort matrix according mutations. Helpful in case of MAF was read along with copy number data. Default FALSE.
#' @param keepGeneOrder logical whether to keep order of given genes. Default FALSE, order according to mutation frequency
#' @param GeneOrderSort logical this is applicable when `keepGeneOrder` is TRUE. Default TRUE
#' @param sampleOrder Manually speify sample names for oncolplot ordering. Default NULL.
#' @param additionalFeature a vector of length two indicating column name in the MAF and the factor level to be highlighted. Provide a list of values for highlighting more than one features
#' @param additionalFeaturePch Default 20
#' @param additionalFeatureCol Default "gray70"
#' @param additionalFeatureCex Default 0.9
#' @param genesToIgnore do not show these genes in Oncoplot. Default NULL.
#' @param removeNonMutated Logical. If \code{TRUE} removes samples with no mutations in the oncoplot for better visualization. Default \code{FALSE}.
#' @param fill Logical. If \code{TRUE} draws genes and samples as blank grids even when they are not altered.
#' @param cohortSize Number of sequenced samples in the cohort. Default all samples from Cohort. You can manually specify the cohort size. Default \code{NULL}
#' @param colors named vector of colors for each Variant_Classification.
#' @param cBioPortal Adds annotations similar to cBioPortals MutationMapper and collapse Variants into Truncating and rest.
#' @param bgCol Background grid color for wild-type (not-mutated) samples. Default "#ecf0f1"
#' @param borderCol border grid color (not-mutated) samples. Default 'white'.
#' @param annoBorderCol border grid color for annotations. Default NA.
#' @param numericAnnoCol color palette used for numeric annotations. Default 'YlOrBr' from RColorBrewer
#' @param drawBox logical whether to draw a box around main matrix. Default FALSE
#' @param fontSize font size for gene names. Default 0.8.
#' @param SampleNamefontSize font size for sample names. Default 1
#' @param titleFontSize font size for title. Default 1.5
#' @param legendFontSize font size for legend. Default 1.2
#' @param annotationFontSize font size for annotations. Default 1.2
#' @param sepwd_genes size of lines seperating genes. Default 0.5
#' @param sepwd_samples size of lines seperating samples. Default 0.25
#' @param writeMatrix writes character coded matrix used to generate the plot to an output file.
#' @param colbar_pathway Draw top column bar with respect to diplayed pathway. Default FALSE.
#' @param showTitle Default TRUE
#' @param titleText Custom title. Default `NULL`
#' @param showPct Default TRUE. Shows percent altered to the right side of the plot.
#' @returns Invisibly returns a list with components
#' 1. `oncomatrix` A matrix used for drawing the oncoplot. Values are numeric coded for each variant classification
#' 2. `vc_legend` A mapping of variant classification to numeric values in the oncomatrix
#' 3. `vc_color` Color coding used for each variant classification
#' @examples
#' laml.maf <- system.file("extdata", "tcga_laml.maf.gz", package = "maftools")
#' laml.clin = system.file('extdata', 'tcga_laml_annot.tsv', package = 'maftools')
#' laml <- read.maf(maf = laml.maf, clinicalData = laml.clin)
#' #Basic onocplot
#' oncoplot(maf = laml, top = 3)
#' #Changing colors for variant classifications (You can use any colors, here in this example we will use a color palette from RColorBrewer)
#' col = RColorBrewer::brewer.pal(n = 8, name = 'Paired')
#' names(col) = c('Frame_Shift_Del','Missense_Mutation', 'Nonsense_Mutation', 'Multi_Hit', 'Frame_Shift_Ins',
#' 'In_Frame_Ins', 'Splice_Site', 'In_Frame_Del')
#' #Color coding for FAB classification; try getAnnotations(x = laml) to see available annotations.
#' fabcolors = RColorBrewer::brewer.pal(n = 8,name = 'Spectral')
#' names(fabcolors) = c("M0", "M1", "M2", "M3", "M4", "M5", "M6", "M7")
#' fabcolors = list(FAB_classification = fabcolors)
#' oncoplot(maf = laml, colors = col, clinicalFeatures = 'FAB_classification', sortByAnnotation = TRUE, annotationColor = fabcolors)
#' @seealso \code{\link{pathways}}
#' @export
oncoplot = oncoplot = function(maf, top = 20, minMut = NULL, genes = NULL, altered = FALSE,
drawRowBar = TRUE, drawColBar = TRUE,
leftBarData = NULL, leftBarLims = NULL, leftBarVline = NULL, leftBarVlineCol = 'gray70',
rightBarData = NULL, rightBarLims = NULL,rightBarVline = NULL, rightBarVlineCol = 'gray70',
topBarData = NULL, topBarLims = NULL, topBarHline = NULL, topBarHlineCol = 'gray70', logColBar = FALSE, includeColBarCN = TRUE,
clinicalFeatures = NULL, annotationColor = NULL, annotationDat = NULL,
pathways = NULL, topPathways = 3,path_order = NULL, selectedPathways = NULL, collapsePathway = FALSE, pwLineCol = "#535c68", pwLineWd = 1, draw_titv = FALSE, titv_col = NULL,
showTumorSampleBarcodes = FALSE, tsbToPIDs = NULL, barcode_mar = 4, barcodeSrt = 90, gene_mar = 5,
anno_height = 1, legend_height = 4,
sortByAnnotation = FALSE, groupAnnotationBySize = TRUE, annotationOrder = NULL,
sortByMutation = FALSE, keepGeneOrder = FALSE, GeneOrderSort = TRUE, sampleOrder = NULL,
additionalFeature = NULL, additionalFeaturePch = 20, additionalFeatureCol = "gray70", additionalFeatureCex = 0.9,
genesToIgnore = NULL, removeNonMutated = FALSE, fill = TRUE, cohortSize = NULL,
colors = NULL, cBioPortal = FALSE, bgCol = "#ecf0f1", borderCol = 'white', annoBorderCol = NA, numericAnnoCol = NULL,
drawBox = FALSE, fontSize = 0.8, SampleNamefontSize = 1, titleFontSize = 1.5, legendFontSize = 1.2, annotationFontSize = 1.2,
sepwd_genes = 0.5, sepwd_samples = 0.25, writeMatrix = FALSE, colbar_pathway = FALSE, showTitle = TRUE, titleText = NULL, showPct = TRUE){
#Total samples
if(is.null(cohortSize)){
totSamps = as.numeric(maf@summary[3,summary])
}else{
totSamps = cohortSize
}
if(!is.null(genes)){ #If user provides a gene list
genes = unique(as.character(genes))
om = createOncoMatrix(m = maf, g = genes, add_missing = fill, cbio = cBioPortal)
numMat = om$numericMatrix
mat_origin = om$oncoMatrix
if(!is.null(pathways)){
stop("Please use genes or pathways.")
}
} else if(!is.null(pathways)){ #If user provides pathway list
if(is(pathways, 'data.frame')){
pathways = data.table::copy(pathways)
colnames(pathways)[1:2] = c('Gene', 'Pathway')
data.table::setDT(x = pathways)
}else if(file.exists(pathways)){
pathways = data.table::fread(file = pathways)
colnames(pathways)[1:2] = c('Gene', 'Pathway')
data.table::setDT(x = pathways)
}else{
pathways = match.arg(arg = pathways, choices = c("sigpw", "smgbp"))
if(pathways == "sigpw"){
pathdb <- system.file("extdata", "oncogenic_sig_patwhays.tsv", package = "maftools")
#message("Summarizing signalling pathways [Sanchez-Vega et al., https://doi.org/10.1016/j.cell.2018.03.035]")
}else{
#message("Summarizing known driver genes [Bailey et al., https://doi.org/10.1016/j.cell.2018.02.060]")
pathdb <- system.file("extdata", "BP_SMGs.txt.gz", package = "maftools")
}
pathways = data.table::fread(file = pathdb, skip = "Gene")
}
#Convert to characters in case of factors
pathways$Gene = as.character(pathways$Gene)
pathways$Pathway = as.character(pathways$Pathway)
if(is.null(selectedPathways)){
pathwayLoad = get_pw_summary(maf = maf, pathways = pathways[,.(Pathway, Gene)])
message("Drawing upto top ", topPathways, " mutated pathways")
if(ncol(pathwayLoad) >= topPathways){
pathways = pathways[Pathway %in% pathwayLoad[1:topPathways, Pathway]]
}else{
pathways = pathways[Pathway %in% pathwayLoad[, Pathway]]
}
}else{
pws_missing = setdiff(selectedPathways, pathways[,.N,Pathway][,Pathway])
if(length(pws_missing) > 0){
warning("Could not find the following pathways: ", paste(pws_missing, collapse = ", "), immediate. = TRUE)
}
pathways = pathways[Pathway %in% selectedPathways]
if(nrow(pathways) == 0){
stop("Zero pathways to plot!")
}
}
#Make sure pathway names are distinct from the gene names
dup_pw_names = intersect(pathways$Pathway, pathways$Gene)
if(length(dup_pw_names) > 0){
warning("Found following pathways with the same name as gene symbols. Renamed them with the suffix _pw\n", paste(dup_pw_names, collapse = ", "), immediate. = TRUE)
unique_pw_names = setdiff(pathways$Pathway, dup_pw_names)
pw_named_list = setNames(nm = c(unique_pw_names, dup_pw_names), object = c(unique_pw_names, paste0(dup_pw_names, "_pw")))
pathways$Pathway = unname(obj = pw_named_list[pathways$Pathway])
}
if(nrow(pathways[duplicated(Gene)]) > 0){
warning("Duplicated genes found across multiple pathways! This might cause issue while plotting.", immediate. = TRUE)
# pathways = pathways[!duplicated(Gene)]
}
genes = as.character(pathways$Gene)
om = createOncoMatrix(m = maf, g = genes, cbio = cBioPortal)
numMat = om$numericMatrix
mat_origin = om$oncoMatrix
#Check for any missing genes and ignore them if necessary
if(length(genes[!genes %in% rownames(numMat)]) > 0){
#warning('Missing following genes from provided pathways ', paste(genes[!genes %in% rownames(numMat)], collapse = ", "))
genes = genes[genes %in% rownames(numMat)]
}
}else if(!is.null(minMut)){
if(altered){
genes = getGeneSummary(x = maf)[order(AlteredSamples, decreasing = TRUE)][,.(Hugo_Symbol, AlteredSamples)]
}else{
genes = getGeneSummary(x = maf)[order(MutatedSamples, decreasing = TRUE)][,.(Hugo_Symbol, MutatedSamples)]
}
colnames(genes)[2] = "mutload"
genes[,fractMutated := mutload/totSamps]
if(minMut > 1){
genes = genes[mutload >= minMut, Hugo_Symbol]
}else{
genes = genes[fractMutated >= minMut, Hugo_Symbol]
}
om = createOncoMatrix(m = maf, g = genes, cbio = cBioPortal)
numMat = om$numericMatrix
mat_origin = om$oncoMatrix
}else { #If user does not provide gene list or MutSig results, draw TOP (default 20) genes
if(altered){
genes = getGeneSummary(x = maf)[order(AlteredSamples, decreasing = TRUE)][1:top, Hugo_Symbol]
}else{
genes = getGeneSummary(x = maf)[1:top, Hugo_Symbol]
}
om = createOncoMatrix(m = maf, g = genes, cbio = cBioPortal)
numMat = om$numericMatrix
mat_origin = om$oncoMatrix
}
#---remove genes from genesToIgnore if any
if(!is.null(genesToIgnore)){
numMat = numMat[!rownames(numMat) %in% genesToIgnore,]
mat_origin = mat_origin[!rownames(mat_origin) %in% genesToIgnore,]
}
tsbs = levels(getSampleSummary(x = maf)[,Tumor_Sample_Barcode])
if(!removeNonMutated){
tsb.include = matrix(data = 0, nrow = nrow(numMat),
ncol = length(tsbs[!tsbs %in% colnames(numMat)]))
colnames(tsb.include) = tsbs[!tsbs %in% colnames(numMat)]
rownames(tsb.include) = rownames(numMat)
numMat = cbind(numMat, tsb.include)
mat_origin = cbind(mat_origin, tsb.include)
}
#If user wannts to keep given gene order
if(keepGeneOrder){
if(fill){
numMat = numMat[genes, , drop = FALSE]
}else{
if(GeneOrderSort){
numMat = sortByGeneOrder(m = numMat, g = genes)
}else{
numMat = numMat[genes, , drop = FALSE]
}
}
mat = mat_origin[rownames(numMat), , drop = FALSE]
mat = mat[,colnames(numMat), drop = FALSE]
}
if(sortByMutation){
numMat_temp = sortByMutation(numMat = numMat, maf = maf)
numMat = numMat[rownames(numMat_temp), colnames(numMat_temp), drop = FALSE]
}
samp_sum = data.table::copy(getSampleSummary(x = maf))
samp_sum[,total := NULL]
if(cBioPortal){
samp_sum = data.table::melt(data = samp_sum, id.vars = "Tumor_Sample_Barcode")
colnames(samp_sum)[2] = "Variant_Classification"
vc = c("Nonstop_Mutation", "Frame_Shift_Del", "Missense_Mutation",
"Nonsense_Mutation", "Splice_Site", "Frame_Shift_Ins", "In_Frame_Del", "In_Frame_Ins")
vc.cbio = c("Truncating", "Truncating", "Missense", "Truncating", "Truncating", "Truncating",
"In-frame", "In-frame")
names(vc.cbio) = vc
samp_sum[,Variant_Classification_temp := vc.cbio[as.character(samp_sum$Variant_Classification)]]
samp_sum$Variant_Classification_temp = ifelse(test = is.na(samp_sum$Variant_Classification_temp), yes = as.character(samp_sum$Variant_Classification), no = samp_sum$Variant_Classification_temp)
samp_sum[,Variant_Classification := as.factor(as.character(Variant_Classification_temp))]
samp_sum[,Variant_Classification_temp := NULL]
samp_sum = samp_sum[,sum(value), .(Tumor_Sample_Barcode, Variant_Classification)]
samp_sum = data.table::dcast(data = samp_sum, formula = Tumor_Sample_Barcode~Variant_Classification, value.var = "V1")
}
if("CNV_total" %in% colnames(samp_sum)){
samp_sum[,CNV_total := NULL]
}
if(!includeColBarCN){
suppressWarnings(samp_sum[,Amp := NULL])
suppressWarnings(samp_sum[,Del := NULL])
}
data.table::setDF(x = samp_sum, rownames = as.character(samp_sum$Tumor_Sample_Barcode))
samp_sum = samp_sum[,-1, drop = FALSE]
#Parse annotations
if(!is.null(clinicalFeatures)){
if(is.null(annotationDat)){
annotation = parse_annotation_dat(annotationDat = maf, clinicalFeatures = clinicalFeatures)
}else{
annotation = parse_annotation_dat(annotationDat = annotationDat, clinicalFeatures = clinicalFeatures)
}
annotation = annotation[colnames(numMat),, drop = FALSE]
if(sortByAnnotation){
numMat = sortByAnnotation(numMat = numMat, maf = maf, anno = annotation, annoOrder = annotationOrder, group = groupAnnotationBySize, isNumeric = FALSE)
}
}
if(!is.null(sampleOrder)){
sampleOrder = as.character(sampleOrder)
sampleOrder = sampleOrder[sampleOrder %in% colnames(numMat)]
if(length(sampleOrder) == 0){
stop("None of the provided samples are present in the input MAF")
}
numMat = numMat[,sampleOrder, drop = FALSE]
}
gene_sum = apply(numMat, 1, function(x) length(x[x!=0]))
percent_alt = paste0(round(100*(apply(numMat, 1, function(x) length(x[x!=0]))/totSamps)), "%")
nonmut_samps = colnames(numMat)[!colnames(numMat) %in% rownames(samp_sum)]
if(length(nonmut_samps) > 0){
samp_sum_missing_samps = matrix(data = 0, nrow = length(nonmut_samps), ncol = ncol(samp_sum))
colnames(samp_sum_missing_samps) = colnames(samp_sum)
rownames(samp_sum_missing_samps) = nonmut_samps
samp_sum = rbind(samp_sum, samp_sum_missing_samps)
}
if(colbar_pathway){
samp_sum = getSampleSummary(subsetMaf(maf = maf, genes = genes))
samp_sum[,total := NULL]
if("CNV_total" %in% colnames(samp_sum)){
samp_sum[,CNV_total := NULL]
}
if(!includeColBarCN){
suppressWarnings(samp_sum[,Amp := NULL])
suppressWarnings(samp_sum[,Del := NULL])
}
data.table::setDF(x = samp_sum, rownames = as.character(samp_sum$Tumor_Sample_Barcode))
samp_sum = samp_sum[,-1]
top_bar_data = t(samp_sum[colnames(numMat),, drop = FALSE])
}else{
top_bar_data = t(samp_sum[colnames(numMat),, drop = FALSE])
}
#VC codes
if(cBioPortal){
if(is.null(colors)){
vc_col = grDevices::adjustcolor(col = c("black", "#33A02C", "brown"), alpha.f = 0.7)
vc_col = c('Truncating' = vc_col[1], 'Missense' = vc_col[2], 'In-frame' = vc_col[3], "Multi_Hit" = "#9b59b6", "pathway" = "#d35400")
}else{
vc_col = colors
if(!"pathway" %in% names(vc_col)){
vc_col = c(vc_col, "pathway" = "#535C68FF")
}
}
}else{
if(is.null(colors)){
vc_col = get_vcColors(websafe = FALSE)
}else{
vc_col = colors
if(!"pathway" %in% names(vc_col)){
vc_col = c(vc_col, "pathway" = "#535C68FF")
}
}
}
vc_codes = update_vc_codes(om_op = om)
vc_col = update_colors(x = vc_codes, y = vc_col)
if(nrow(numMat) == 1){
stop("Oncoplot requires at-least two genes for plottng.")
}
if(!is.null(pathways)){
vc_codes = c(vc_codes, '99' = 'pathway')
#Order genes in matrix by most frequently altered gene, followed by pathway size
temp_dat = data.table::data.table(Gene = rownames(numMat), percent_alt)
temp_dat = merge(temp_dat, pathways, all.x = TRUE)
temp_dat[is.na(temp_dat)] = "Unknown"
temp_dat$pct_alt = as.numeric(gsub(pattern = "%", replacement = "", x = temp_dat$percent_alt))
temp_dat = split(temp_dat, as.factor(as.character(temp_dat$Pathway)))
temp_dat_freq = lapply(temp_dat, function(x){
data.table::data.table(n = nrow(x), max_alt = max(x$pct_alt))
})
temp_dat_freq = data.table::rbindlist(l = temp_dat_freq, use.names = TRUE, fill = TRUE, idcol = "Pathway")
temp_dat_ord = temp_dat_freq[order(max_alt, n, decreasing = TRUE)][, Pathway]
temp_dat = lapply(temp_dat, function(x){x[order(pct_alt, decreasing = TRUE)]})[temp_dat_ord]
if("Unknown" %in% names(temp_dat)){
temp_dat_ord = c(grep(pattern = "Unknown", x = names(temp_dat), invert = TRUE, value = TRUE), "Unknown")
temp_dat = temp_dat[temp_dat_ord]
}
if(!is.null(path_order)){
path_order = intersect(x = path_order, names(temp_dat) )
temp_dat = temp_dat[path_order]
if(length(temp_dat) == 0){
stop("None of the pathways match the provided order!")
}
}
nm = lapply(seq_along(temp_dat), function(i){
x = numMat[temp_dat[[i]]$Gene,, drop = FALSE]
x = rbind(x, apply(x, 2, function(y) ifelse(test = sum(y) == 0, yes = 0, no = 99)))
if(names(temp_dat)[i] %in% rownames(x)){
stop("pathway name ", names(temp_dat)[i], " is same as one of its member genes! Please rename it.\ne.g, ", names(temp_dat)[i], " to ", names(temp_dat)[i], "_pw")
}
rownames(x)[nrow(x)] = names(temp_dat)[i]
x
})
numMat = do.call(rbind, nm)
#Add pathway information to the character matrix
mat_origin_path = rownames(numMat)[!rownames(numMat) %in% rownames(mat_origin)]
mat_origin_path = numMat[mat_origin_path,, drop = FALSE]
mat_origin_path[mat_origin_path == 0] = ""
mat_origin_path[mat_origin_path == "99"] = "pathway"
mat_origin_path = mat_origin_path[,colnames(mat_origin), drop = FALSE]
if(collapsePathway){
mat_origin = mat_origin_path
numMat = numMat[rownames(mat_origin),, drop = FALSE]
row_ord = names(sort(apply(numMat, 1, function(x) length(x[x == 0])), decreasing = FALSE))
numMat = numMat[row_ord,, drop = FALSE]
tnumMat = t(numMat) #transposematrix
numMat = t(tnumMat[do.call(order, c(as.list(as.data.frame(tnumMat)), decreasing = TRUE)), ]) #sort
if(sortByAnnotation){
numMat = sortByAnnotation(numMat = numMat, maf = maf, anno = annotation, annoOrder = annotationOrder, group = groupAnnotationBySize, isNumeric = FALSE)
}
}else{
mat_origin = rbind(mat_origin, mat_origin_path)
}
mat_origin = mat_origin[rownames(numMat), colnames(numMat), drop = FALSE]
}
#Plot layout
plot_layout(clinicalFeatures = clinicalFeatures, drawRowBar = drawRowBar, anno_height = anno_height,
drawColBar = drawColBar, draw_titv = draw_titv, exprsTbl = leftBarData, legend_height = legend_height)
#01: Draw scale axis for left barplot table
leftBarTitle = NULL
if(!is.null(leftBarData)){
leftBarTitle = colnames(leftBarData)[2]
colnames(leftBarData) = c('genes', 'exprn')
data.table::setDF(x = leftBarData, rownames = as.character(leftBarData$genes))
missing_samps = rownames(numMat)[!rownames(numMat) %in% rownames(leftBarData)]
if(length(missing_samps) > 0){
temp_data = data.frame(
row.names = missing_samps,
genes = missing_samps,
exprn = rep(0, length(missing_samps)),
stringsAsFactors = FALSE
)
leftBarData = rbind(leftBarData, temp_data)
}
leftBarData = leftBarData[rownames(numMat),, drop = FALSE]
#leftBarData = leftBarData[genes,, drop = FALSE]
if(is.null(leftBarLims)){
exprs_bar_lims = round(range(leftBarData$exprn, na.rm = TRUE), digits = 2)
}else{
exprs_bar_lims = leftBarLims
}
if(drawColBar){
par(mar = c(0.25 , 0, 1, 2), xpd = TRUE)
plot(x = NA, y = NA, type = "n", axes = FALSE,
xlim = exprs_bar_lims, ylim = c(0, 1), xaxs = "i")
}
}
#02: Draw top bar plot
if(drawColBar & is.null(topBarData)){
top_bar_data = top_bar_data[,colnames(numMat), drop = FALSE]
if(drawRowBar){
par(mar = c(0.25 , gene_mar, 2, 3), xpd = TRUE)
}else{
par(mar = c(0.25 , gene_mar, 2, 5), xpd = TRUE)
}
if(logColBar){
top_bar_data = apply(top_bar_data, 2, function(x) {
x_fract = x / sum(x)
x_log_total = log10(sum(x))
x_fract * x_log_total
})
top_bar_data[is.infinite(top_bar_data)] = 0
}
if(is.null(topBarLims)){
topBarLims = round(c(0, max(colSums(x = top_bar_data, na.rm = TRUE))), digits = 2)
}
plot(x = NA, y = NA, type = "n", xlim = c(0,ncol(top_bar_data)), ylim = topBarLims,
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i")
axis(side = 2, at = c(0, max(topBarLims, na.rm = TRUE)), las = 2, line = 0.5)
for(i in 1:ncol(top_bar_data)){
x = top_bar_data[,i]
names(x) = rownames(top_bar_data)
x = x[!x == 0]
if(length(x) > 0){
rect(xleft = i-1, ybottom = c(0, cumsum(x)[1:(length(x)-1)]), xright = i-0.1,
ytop = cumsum(x), col = vc_col[names(x)], border = NA, lwd = 0, xpd = FALSE)
}
}
if(logColBar){
mtext(text = "(log10)", side = 2, line = 2, cex = 0.6)
}else{
mtext(text = "TMB", side = 2, line = 2, cex = 0.6)
}
if(!is.null(topBarHline)){
abline(h = topBarHline, lty = 2, col = topBarHlineCol, xpd = FALSE)
mtext(text = topBarHline, side = 2, las = 2, cex = 0.6, line = 0.7)
}
}else if(!is.null(topBarData) & drawColBar){
# Draw extra clinical data in top
if(drawRowBar){
par(mar = c(0.25 , gene_mar, 2, 3), xpd = TRUE)
}else{
par(mar = c(0.25 , gene_mar, 2, 5), xpd = TRUE)
}
if(is.data.frame(topBarData)){
extdata = data.table::copy(x = topBarData)
colnames(extdata)[1] = "Tumor_Sample_Barcode"
rownames(extdata) <- extdata$Tumor_Sample_Barcode
}else{
message("Top barplot data is derived from clinical column: ", topBarData)
extdata <- data.frame(maf@clinical.data)[,c("Tumor_Sample_Barcode", topBarData)]
rownames(extdata) <- extdata$Tumor_Sample_Barcode
}
topbar_title = colnames(extdata)[2]
colnames(extdata)[2] = "barheights"
extdata[, "barheights"] <- as.numeric(as.character(extdata[, "barheights"]))
if(!all(colnames(top_bar_data) %in% rownames(extdata))){
missing_topbars = colnames(top_bar_data)[!colnames(top_bar_data) %in% rownames(extdata)]
warning("Missing column barplot data for ", length(missing_topbars), " samples")
extdata = rbind(extdata, data.frame(row.names = missing_topbars, Tumor_Sample_Barcode = missing_topbars, barheights = 0,stringsAsFactors = FALSE))
}
extdata <- extdata[colnames(top_bar_data),]
if(is.null(topBarLims)){
topBarLims = round(range(extdata[,"barheights"], na.rm = TRUE), digits = 2)
}
plot(x = NA, y = NA, type = "n", xlim = c(0,nrow(extdata)), ylim = topBarLims,
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i")
axis(side = 2, at = topBarLims, las = 2, line = 0.5)
for(i in 1:nrow(extdata)){
graphics::rect(xleft = i-1, xright = i-0.1, ybottom = 0, ytop = as.numeric(extdata[i, "barheights"]),
col = "#535c68", border = NA, lwd = 0, xpd = FALSE)
}
title(ylab = topbar_title, font = 1, cex.lab = legendFontSize, xpd = TRUE)
if(!is.null(topBarHline)){
abline(h = topBarHline, lty = 2, col = topBarHlineCol)
}
}
#03: Draw scale for right barplot
rightBarTitle = NULL
if(drawRowBar){
if(is.null(rightBarData)){
side_bar_lims = c(0, max(unlist(apply(numMat, 1, function(x) cumsum(table(x[x!=0])))), na.rm = TRUE))
rightBarTitle = "No. of samples"
}else{
rightBarTitle = colnames(rightBarData)[2]
colnames(rightBarData) = c('genes', 'exprn')
data.table::setDF(x = rightBarData, rownames = as.character(rightBarData$genes))
missing_samps = rownames(numMat)[!rownames(numMat) %in% rownames(rightBarData)]
if(length(missing_samps) > 0){
temp_data = data.frame(
row.names = missing_samps,
genes = missing_samps,
exprn = rep(0, length(missing_samps)),
stringsAsFactors = FALSE
)
rightBarData = rbind(rightBarData, temp_data)
}
rightBarData = rightBarData[rownames(numMat),, drop = FALSE]
if(is.null(rightBarLims)){
side_bar_lims = round(range(rightBarData$exprn, na.rm = TRUE), digits = 2)
}else{
side_bar_lims = rightBarLims
}
}
if(drawColBar){
par(mar = c(0.25 , 0, 1, 2), xpd = TRUE)
plot(x = NA, y = NA, type = "n", axes = FALSE,
xlim = side_bar_lims, ylim = c(0, 1), xaxs = "i")
}
}
#Draw expression barplot
if(!is.null(leftBarData)){
if(showTumorSampleBarcodes){
if(!drawRowBar & !drawColBar){
par(mar = c(barcode_mar, 1, 2.5, 0), xpd = TRUE)
}else if(!drawRowBar & drawColBar){
par(mar = c(barcode_mar, 1, 0, 5), xpd = TRUE)
}else if(drawRowBar & !drawColBar){
par(mar = c(barcode_mar, 1, 2.5, 0), xpd = TRUE)
} else{
par(mar = c(barcode_mar, 1, 0, 0), xpd = TRUE)
}
}else{
if(!drawRowBar & !drawColBar){
par(mar = c(0.5, 1, 2.5, 0), xpd = TRUE)
}else if(!drawRowBar & drawColBar){
par(mar = c(0.5, 1, 0, 5), xpd = TRUE)
}else if(drawRowBar & !drawColBar){
par(mar = c(0.5, 1, 2.5, 0), xpd = TRUE)
} else{
par(mar = c(0.5, 1, 0, 0), xpd = TRUE)
}
}
plot(x = NA, y = NA, type = "n", xlim = rev(-exprs_bar_lims), ylim = c(0, nrow(leftBarData)),
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i", yaxs = "i") #
for(i in 1:nrow(leftBarData)){
x = rev(leftBarData$exprn)[i]
rect(ybottom = i-1, xleft = -x, ytop = i-0.1,
xright = 0, col = "#535c68", border = NA, lwd = 0)
}
axis(side = 3, at = rev(-exprs_bar_lims), outer = FALSE, line = 0.25, labels = rev(exprs_bar_lims))
mtext(text = leftBarTitle, side = 3, line = 0.50, cex = 0.6)
if(!is.null(leftBarVline)){
abline(h = leftBarVline, lty = 2, col = leftBarVlineCol, xpd = FALSE)
}
}
#04: Draw the main matrix
if(showTumorSampleBarcodes){
if(!drawRowBar & !drawColBar){
par(mar = c(barcode_mar, gene_mar, 2.5, 5), xpd = TRUE)
}else if(!drawRowBar & drawColBar){
par(mar = c(barcode_mar, gene_mar, 0, 5), xpd = TRUE)
}else if(drawRowBar & !drawColBar){
par(mar = c(barcode_mar, gene_mar, 2.5, 3), xpd = TRUE)
} else{
par(mar = c(barcode_mar, gene_mar, 0, 3), xpd = TRUE)
}
}else{
if(!drawRowBar & !drawColBar){
par(mar = c(0.5, gene_mar, 2.5, 5), xpd = TRUE)
}else if(!drawRowBar & drawColBar){
par(mar = c(0.5, gene_mar, 0, 5), xpd = TRUE)
}else if(drawRowBar & !drawColBar){
par(mar = c(0.5, gene_mar, 2.5, 3), xpd = TRUE)
} else{
par(mar = c(0.5, gene_mar, 0, 3), xpd = TRUE)
}
}
nm = t(apply(numMat, 2, rev))
nm[nm == 0] = NA
image(x = 1:nrow(nm), y = 1:ncol(nm), z = nm, axes = FALSE, xaxt="n", yaxt="n",
xlab="", ylab="", col = "white") #col = "#FC8D62"
#Plot for all variant classifications
vc_codes_temp = vc_codes[!vc_codes %in% om$cnvc]
for(i in 2:length(names(vc_codes_temp))){
vc_code = vc_codes_temp[i]
col = vc_col[vc_code]
nm = t(apply(numMat, 2, rev))
#print(names(vc_code))
nm[nm != names(vc_code)] = NA
#print(paste0(vc_code, " ", col))
#Suppress warning due to min/max applied to a vector of NAs; Issue: #286
#This is an harmless warning as matrix is looped over all VC's and missing VC's form NA's (which are plotted in gray)
suppressWarnings(image(x = 1:nrow(nm), y = 1:ncol(nm), z = nm, axes = FALSE, xaxt="n", yaxt="n",
xlab="", ylab="", col = col, add = TRUE))
}
#Add blanks
nm = t(apply(numMat, 2, rev))
nm[nm != 0] = NA
image(x = 1:nrow(nm), y = 1:ncol(nm), z = nm, axes = FALSE, xaxt="n", yaxt="n", xlab="", ylab="", col = bgCol, add = TRUE)
#Add CNVs if any
mat_origin = mat_origin[rownames(numMat), colnames(numMat), drop = FALSE]
if(writeMatrix){
write.table(mat_origin, "onco_matrix.txt", sep = "\t", quote = FALSE)
}
mo = t(apply(mat_origin, 2, rev))
##Complex events (mutated as well as CN altered)
complex_events = unique(grep(pattern = ";", x = mo, value = TRUE))
if(length(complex_events) > 0){
for(i in 1:length(complex_events)){
ce = complex_events[i]
#mo = t(apply(mat_origin, 2, rev))
ce_idx = which(mo == ce, arr.ind = TRUE)
ce = unlist(strsplit(x = ce, split = ";", fixed = TRUE))
nm_temp = matrix(NA, nrow = nrow(nm), ncol = ncol(nm))
nm_temp[ce_idx] = 0
image(x = 1:nrow(nm_temp), y = 1:ncol(nm_temp), z = nm_temp, axes = FALSE, xaxt="n",
yaxt="n", xlab="", ylab="", col = vc_col[ce[2]], add = TRUE)
ce_idx = which(t(nm_temp) == 0, arr.ind = TRUE)
for(i in seq_len(nrow(ce_idx))){
rowi = ce_idx[i,1]
coli = ce_idx[i,2]
rect(xleft = coli-0.5, ybottom = rowi-0.25, xright = coli+0.5, ytop = rowi+0.25, col = vc_col[ce[1]], border = NA, lwd = 0)
}
}
}
for(cnevent in om$cnvc){
cn_idx = which(mo == cnevent, arr.ind = TRUE)
if(nrow(cn_idx) > 0){
nm_temp = matrix(NA, nrow = nrow(nm), ncol = ncol(nm))
nm_temp[cn_idx] = 0
image(x = 1:nrow(nm_temp), y = 1:ncol(nm_temp), z = nm_temp, axes = FALSE, xaxt="n",
yaxt="n", xlab="", ylab="", col = bgCol, add = TRUE)
cn_idx = which(t(nm_temp) == 0, arr.ind = TRUE)
for(i in seq_len(nrow(cn_idx))){
rowi = cn_idx[i,1]
coli = cn_idx[i,2]
rect(xleft = coli-0.5, ybottom = rowi-0.25, xright = coli+0.5, ytop = rowi+0.25, col = vc_col[cnevent], border = NA, lwd = 0)
}
}
}
#Add grids
abline(h = (0:ncol(nm)) + 0.5, col = borderCol, lwd = sepwd_genes, xpd = FALSE)
abline(v = (0:nrow(nm)) + 0.5, col = borderCol, lwd = sepwd_samples, xpd = FALSE)
#Add boxes if pathways are opted
if(!is.null(pathways)){
temp_dat = lapply(seq_along(temp_dat), function(i){
data.table::rbindlist(list(temp_dat[[i]], data.table::data.table(Gene = names(temp_dat)[i])), use.names = TRUE, fill = TRUE)
})
temp_dat = data.table::rbindlist(l = temp_dat, use.names = TRUE, fill = TRUE)
#return(temp_dat)
temp_dat[, row_id := nrow(temp_dat):1]
temp_dat = split(temp_dat, as.factor(temp_dat$Pathway))
if(!collapsePathway){
lapply(temp_dat, function(td){
rect(xleft = 0.5, ybottom = min(td$row_id)-0.499, xright = nrow(nm)+0.5, ytop = max(td$row_id)+0.499, border = pwLineCol, lwd = pwLineWd)
})
}
#New percent alt with pathways
percent_alt = rev(paste0(apply(nm, 2, function(x){
round(length(x[is.na(x)])/totSamps * 100)
}), "%"))
}
#Draw if any additional features are requested
additionalFeature_legend = FALSE
if(!is.null(additionalFeature)){
#If its a list, multi features to be mapped
if(is.list(additionalFeature)){
for(af_idx in seq_along(additionalFeature)){
af = additionalFeature[[af_idx]]
if(length(af) < 2){
stop("additionalFeature must be of length two. See ?oncoplot for details.")
}
af_dat = subsetMaf(maf = maf, genes = rownames(numMat), tsb = colnames(numMat), fields = af[1], includeSyn = FALSE, mafObj = FALSE)
if(length(which(colnames(af_dat) == af[1])) == 0){
message(paste0("Column ", af[1], " not found in maf. Here are available fields.."))
print(getFields(maf))
stop()
}
colnames(af_dat)[which(colnames(af_dat) == af[1])] = 'temp_af'
af_dat = af_dat[temp_af %in% af[2]]
if(nrow(af_dat) == 0){
warning(paste0("No samples are enriched for ", af[2], " in ", af[1]))
}else{
af_mat = data.table::dcast(data = af_dat, Tumor_Sample_Barcode ~ Hugo_Symbol, value.var = "temp_af", fun.aggregate = length)
af_mat = as.matrix(af_mat, rownames = "Tumor_Sample_Barcode")
nm = t(apply(numMat, 2, rev))
if(length(additionalFeaturePch) != length(additionalFeature)){
warning("Provided pch for additional features are recycled")
additionalFeaturePch = rep(additionalFeaturePch, length(additionalFeature))
}
if(length(additionalFeatureCol) != length(additionalFeature)){
warning("Provided colors for additional features are recycled")
additionalFeatureCol = rep(additionalFeatureCol, length(additionalFeature))
}
lapply(seq_len(nrow(af_mat)), function(i){
af_i = af_mat[i,, drop = FALSE]
af_i_genes = colnames(af_i)[which(af_i > 0)]
af_i_sample = rownames(af_i)
lapply(af_i_genes, function(ig){
af_i_mat = matrix(c(which(rownames(nm) == af_i_sample),
which(colnames(nm) == ig)),
nrow = 1)
points(af_i_mat, pch = additionalFeaturePch[af_idx], col= additionalFeatureCol[af_idx], cex = additionalFeatureCex)
})
})
additionalFeature_legend = TRUE
}
}
}else{
if(length(additionalFeature) < 2){
stop("additionalFeature must be of length two. See ?oncoplot for details.")
}
af_dat = subsetMaf(maf = maf, genes = rownames(numMat), tsb = colnames(numMat), fields = additionalFeature[1], includeSyn = FALSE, mafObj = FALSE)
if(length(which(colnames(af_dat) == additionalFeature[1])) == 0){
message(paste0("Column ", additionalFeature[1], " not found in maf. Here are available fields.."))
print(getFields(maf))
stop()
}
colnames(af_dat)[which(colnames(af_dat) == additionalFeature[1])] = 'temp_af'
af_dat = af_dat[temp_af %in% additionalFeature[2]]
if(nrow(af_dat) == 0){
warning(paste0("No samples are enriched for ", additionalFeature[2], " in ", additionalFeature[1]))
}else{
af_mat = data.table::dcast(data = af_dat, Tumor_Sample_Barcode ~ Hugo_Symbol, value.var = "temp_af", fun.aggregate = length)
af_mat = as.matrix(af_mat, rownames = "Tumor_Sample_Barcode")
nm = t(apply(numMat, 2, rev))
lapply(seq_len(nrow(af_mat)), function(i){
af_i = af_mat[i,, drop = FALSE]
af_i_genes = colnames(af_i)[which(af_i > 0)]
af_i_sample = rownames(af_i)
lapply(af_i_genes, function(ig){
af_i_mat = matrix(c(which(rownames(nm) == af_i_sample),
which(colnames(nm) == ig)),
nrow = 1)
points(af_i_mat, pch = additionalFeaturePch, col= additionalFeatureCol, cex = additionalFeatureCex)
})
})
additionalFeature_legend = TRUE
}
}
}
# #Draw grids for the samples
# if(!is.null(gridFeature)){
# To add: https://github.com/PoisonAlien/maftools/issues/528
# }
#Add box border
if (drawBox) {
box(lty = 'solid', col = '#535c68', lwd = 2)
}
mtext(text = colnames(nm), side = 2, at = 1:ncol(nm),
font = 3, line = 0.4, cex = fontSize, las = 2)
if(showPct){
mtext(text = rev(percent_alt), side = 4, at = 1:ncol(nm),
font = 1, line = 0.4, cex = fontSize, las = 2, adj = 0.15)
}
if(showTumorSampleBarcodes){
if(!is.null(tsbToPIDs)){
if(is.data.frame(x = tsbToPIDs)){
colnames(x = tsbToPIDs) = c("Tumor_Sample_Barcode", "PID")
rownames(tsbToPIDs) = tsbToPIDs$Tumor_Sample_Barcode
}else{
clin = getClinicalData(x = maf)
if(!tsbToPIDs %in% colnames(clin)){
stop(tsbToPIDs, " does not exist in the clinical data!")
}
data.table::setDF(x = clin, rownames = clin$Tumor_Sample_Barcode)
tsbToPIDs = clin[,c("Tumor_Sample_Barcode", tsbToPIDs)]
colnames(x = tsbToPIDs) = c("Tumor_Sample_Barcode", "PID")
}
if(all(rownames(nm) %in% rownames(tsbToPIDs))){
text(x =1:nrow(nm), y = par("usr")[3] - 0.2,
labels = tsbToPIDs[rownames(nm), 2], srt = barcodeSrt, font = 1, cex = SampleNamefontSize, adj = 1)
}else{
missing_tsb = setdiff(rownames(nm), rownames(tsbToPIDs))
warning("Follwing Tumor_Sample_Barcodes are missing corresponding PIDs: ", paste(missing_tsb, collapse = ", "))
text(x =1:nrow(nm), y = par("usr")[3] - 0.2,
labels = tsbToPIDs[rownames(nm), 2], srt = barcodeSrt, font = 1, cex = SampleNamefontSize, adj = 1)
}
}else{
text(x =1:nrow(nm), y = par("usr")[3] - 0.2,
labels = rownames(nm), srt = barcodeSrt, font = 1, cex = SampleNamefontSize, adj = 1)
}
}
#05: Draw right side barplot
if(drawRowBar){
if(showTumorSampleBarcodes){
if(!drawRowBar || !drawColBar){
par(mar = c(barcode_mar, 0, 2.5, 1), xpd = TRUE)
}else{
par(mar = c(barcode_mar, 0, 0, 1), xpd = TRUE)
}
}else{
if(!drawRowBar || !drawColBar){
par(mar = c(0.5, 0, 2.5, 1), xpd = TRUE)
}else{
par(mar = c(0.5, 0, 0, 1), xpd = TRUE)
}
}
if(is.null(rightBarData)){
#side_bar_data = apply(numMat, 1, function(x) table(x[x!=0]))
side_bar_data = lapply(seq_len(nrow(numMat)), function(i) {
xi = numMat[i, ]
table(xi[!xi == 0])
})
names(side_bar_data) = rownames(numMat)
plot(x = NA, y = NA, type = "n", xlim = side_bar_lims, ylim = c(0, length(side_bar_data)),
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i", yaxs = "i") #
for(i in 1:length(side_bar_data)){
x = rev(side_bar_data)[[i]]
if(length(x) > 0){
rect(ybottom = i-1, xleft = c(0, cumsum(x)[1:(length(x)-1)]), ytop = i-0.1,
xright = cumsum(x), col = vc_col[vc_codes[names(x)]], border = NA, lwd = 0)
}
}
if(!is.null(rightBarVline)){
abline(v = rightBarVline, lty = 2, col = rightBarVlineCol, xpd = FALSE)
}
axis(side = 3, at = side_bar_lims, outer = FALSE, line = 0.25)
mtext(text = rightBarTitle, side = 3, line = 0.5, cex = 0.6)
}else{
plot(x = NA, y = NA, type = "n", xlim = side_bar_lims, ylim = c(0, nrow(rightBarData)),
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i", yaxs = "i") #
for(i in 1:nrow(rightBarData)){
x = rev(rightBarData$exprn)[i]
rect(ybottom = i-1, xleft = x, ytop = i-0.1,
xright = 0, col = "#535c68", border = NA, lwd = 0)
}
if(!is.null(rightBarVline)){
abline(v = rightBarVline, lty = 2, col = rightBarVlineCol, xpd = FALSE)
}
axis(side = 3, at = side_bar_lims, outer = FALSE, line = 0.25)
mtext(text = rightBarTitle, side = 3, line = 0.5, cex = 0.6)
}
}
#06: Plot annotations if any
if(!is.null(clinicalFeatures)){
annotation = annotation[colnames(numMat), ncol(annotation):1, drop = FALSE]
if(is.null(annotationColor)){
#In case no-color codes are provided
annotationColor = get_anno_cols(ann = annotation)
}else{
#If provided - generate a continuous palette for continuous variables
annotationColor = annotationColor[colnames(annotation)]
annotationColor_temp = lapply(seq_along(annotationColor), function(idx){
if(is.numeric(annotation[,idx])){
#print(annotationColor[[idx]])
seq_col_pal = c("Blues", "BuGn", "BuPu", "GnBu", "Greens", "Greys", "Oranges",
"OrRd", "PuBu", "PuBuGn", "PuRd", "Purples", "RdPu", "Reds",
"YlGn", "YlGnBu", "YlOrBr", "YlOrRd")
col_match = annotationColor[idx] %in% seq_col_pal
if(!col_match){
stop("numeric annotation color for ", names(annotationColor)[idx] , " must be a sequential color palette!\n", paste(seq_col_pal, collapse = " "))
}
x = annotation[,idx]
x_unique = unique(x)
numericAnnoCol = RColorBrewer::brewer.pal(n = 9, name = annotationColor[[idx]])
ann_lvls_cols = colorRampPalette(numericAnnoCol)(length(x_unique))
names(ann_lvls_cols) = x_unique[order(x_unique, na.last = TRUE)]
ann_lvls_cols = ann_lvls_cols[as.character(x)]
}else{
ann_lvls_cols = annotationColor[[idx]]
ann_lvls_cols = ann_lvls_cols[as.character(annotation[,idx])]
}
ann_lvls_cols
})
names(annotationColor_temp) = names(annotationColor)
#dirty way to check which columns are of numeric
temp_col_class = c()
for(i in 1:ncol(annotation)){
if(is.numeric(annotation[,i])){
temp_col_class = c(temp_col_class, "numeric")
}else{
temp_col_class = c(temp_col_class, "char")
}
}
annotationColor = annotationColor_temp
attr(annotationColor, "class") = temp_col_class
}
ann_classes = attributes(annotationColor)$class
#return(list(annotation, annotationColor))
if(!is.null(leftBarData)){
plot.new()
}
if(!drawRowBar){
par(mar = c(0, gene_mar, 0, 5), xpd = TRUE)
}else{
par(mar = c(0, gene_mar, 0, 3), xpd = TRUE)
}
image(x = 1:nrow(annotation), y = 1:ncol(annotation), z = matrix(data = 0, nrow = nrow(annotation), ncol = ncol(annotation)),
axes = FALSE, xaxt="n", yaxt="n", bty = "n",
xlab="", ylab="", col = "gray90") #col = "#FC8D62"
for(i in 1:ncol(annotation)){
temp_anno = annotation
if(ncol(temp_anno) > 1){
rm_idx = setdiff(y = i, x = 1:ncol(temp_anno))
temp_anno[,rm_idx] = NA
}
if(ann_classes[i] == "numeric"){
temp_anno = apply(temp_anno, 2, as.numeric)
suppressWarnings(image(x = 1:nrow(temp_anno), y = 1:ncol(temp_anno), z = temp_anno,
axes = FALSE, xaxt="n", yaxt="n", xlab="", ylab="", col = annotationColor[[i]], add = TRUE))
}else{
temp_lvls = unique(names(annotationColor[[i]])) #unique(temp_anno[,i])
temp_lvls = temp_lvls[complete.cases(temp_lvls)] #Remove NAs
temp_lvls_cols = annotationColor[[i]][temp_lvls]
temp_anno[,i] = match(x = temp_anno[,i], table = unique(temp_lvls), nomatch = NA)
temp_anno = apply(temp_anno, 2, as.numeric)
suppressWarnings(image(x = 1:nrow(temp_anno), y = 1:ncol(temp_anno), z = temp_anno,
axes = FALSE, xaxt="n", yaxt="n", xlab="", ylab="", col = temp_lvls_cols, add = TRUE))
}
}
#Add grids
rect(xleft = 0.5, ybottom = (0:(ncol(annotation))) + 0.5, xright = nrow(annotation)+0.5, ytop = (0:(ncol(annotation))) + 0.5, border = annoBorderCol, lwd = 0.7)
abline(v = (0:(nrow(nm))) + 0.5, col = annoBorderCol, lwd = sepwd_samples, xpd = FALSE)
mtext(text = colnames(annotation), side = 4,
font = 1, line = 0.4, cex = fontSize, las = 2, at = 1:ncol(annotation))
if(drawRowBar){
plot.new()
}
}
#07: Draw TiTv plot
if(draw_titv){
titv_dat = titv(maf = maf, useSyn = TRUE, plot = FALSE)
titv_dat = titv_dat$fraction.contribution
data.table::setDF(x = titv_dat, rownames = as.character(titv_dat$Tumor_Sample_Barcode))
titv_dat = titv_dat[,-1]
titv_dat = t(titv_dat[colnames(numMat), ])
missing_samps = colnames(numMat)[!colnames(numMat) %in% colnames(titv_dat)]
if(length(missing_samps) > 0){
temp_data = matrix(data = 0, nrow = 6, ncol = length(missing_samps))
colnames(temp_data) = missing_samps
titv_dat = cbind(titv_dat, temp_data)
titv_dat = titv_dat[,colnames(numMat)]
}
if(!is.null(leftBarData)){
plot.new()
}
if(!drawRowBar){
par(mar = c(0, gene_mar, 0, 5), xpd = TRUE)
}else{
par(mar = c(0, gene_mar, 0, 3), xpd = TRUE)
}
plot(x = NA, y = NA, type = "n", xlim = c(0,ncol(titv_dat)), ylim = c(0, 100),
axes = FALSE, frame.plot = FALSE, xlab = NA, ylab = NA, xaxs = "i")
if(is.null(titv_col)){
titv_col = get_titvCol(alpha = 1)
}
for(i in 1:ncol(titv_dat)){
x = titv_dat[,i]
x = x[x > 0]
if(length(x) > 0){
rect(xleft = i-1, ybottom = c(0, cumsum(x)[1:(length(x)-1)]), xright = i-0.1,
ytop = cumsum(x), col = titv_col[names(x)], border = NA, lwd = 0)
}else{
rect(xleft = i-1, ybottom = c(0, 100), xright = i-0.1,
ytop = 100, col = "gray70", border = NA, lwd = 0)
}
}
if(drawRowBar){
par(mar = c(0, 0, 1, 6), xpd = TRUE)
plot(NULL,ylab='',xlab='', xlim=0:1, ylim=0:1, axes = FALSE)
lep = legend("topleft", legend = names(titv_col),
col = titv_col, border = NA, bty = "n",
ncol= 2, pch = 15, xpd = TRUE, xjust = 0, yjust = 0,
cex = legendFontSize)
}else{
mtext(text = names(titv_col)[1:3], side = 2, at = c(25, 50, 75),
font = 1, line = 0.4, cex = fontSize, las = 2, col = titv_col[1:3])
mtext(text = names(titv_col)[4:6], side = 4, at = c(25, 50, 75),
font = 1, line = 0.4, cex = fontSize, las = 2, col = titv_col[4:6], adj = 0.15)
}
}
#08: Add legends
par(mar = c(0, 0.5, 0, 0), xpd = TRUE)
plot(NULL,ylab='',xlab='', xlim=0:1, ylim=0:1, axes = FALSE)
leg_classes = vc_col[vc_codes[2:length(vc_codes)]]
leg_classes_pch = rep(15, length(leg_classes))
if(additionalFeature_legend){
if(!is(object = additionalFeature, class2 = "list")){
if(length(additionalFeature) < 2){
stop("additionalFeature must be of length two. See ?oncoplot for details.")
}else{
additionalFeature = list(additionalFeature)
}
}
if(length(additionalFeaturePch) != length(additionalFeature)){
additionalFeaturePch = rep(additionalFeaturePch, length(additionalFeature))
}
if(length(additionalFeatureCol) != length(additionalFeature)){
additionalFeatureCol = rep(additionalFeatureCol, length(additionalFeature))
}
for(af_idx in 1:length(additionalFeature)){
af = additionalFeature[[af_idx]]
leg_classes = c(leg_classes,additionalFeatureCol[af_idx])
names(leg_classes)[length(leg_classes)] = paste(af, collapse = ":")
leg_classes_pch = c(leg_classes_pch, additionalFeaturePch[af_idx])
}
}
if(collapsePathway){
leg_classes = leg_classes['pathway']
}
lep = legend("topleft", legend = names(leg_classes),
col = leg_classes, border = NA, bty = "n",
ncol= 2, pch = leg_classes_pch, xpd = TRUE, xjust = 0, yjust = 0, cex = legendFontSize)
x_axp = 0+lep$rect$w
if(!is.null(clinicalFeatures)){
#Check which of the clinicalFeatures are of numeric
ann_classes = attributes(annotationColor)$class
num_start = 0.1 #numeric legend starts from here
num_width = 0.1
num_incr = 0.2
for(i in seq_along(ann_classes)){
#x = unique(annotation[,i])
temp_lvls = unique(names(annotationColor[[i]]))
temp_lvls = temp_lvls[complete.cases(temp_lvls)]
temp_lvls_cols = annotationColor[[i]][temp_lvls]
x = annotationColor[[i]]
is_num = ann_classes[i] == "numeric"
xt = names(x)
if(is_num){
xt = as.numeric(xt)
xt = xt[!is.na(xt)]
xt = xt[!is.infinite(xt)]
xt = sort(xt)
x = x[as.character(xt)]
x_range = round(range(xt), 2)
image(
y = c(0.1, 0.2),
x = seq(num_start, num_start+num_width, length.out = 200),
z = t(matrix(seq(0,1,length.out = 200), nrow= 1)),
col = x, xlim = c(0, 1), ylim = c(0, 1), axes = FALSE, xlab = NA, ylab = NA, add = TRUE
)
text(x = c(num_start, num_start+num_width), y = 0.2, labels = x_range, pos = 3)
text(x = num_start+(num_width/2), y = 0.1, labels = names(annotationColor)[i], pos = 1, adj = 1)
num_start = num_start + num_incr
}else{
if(length(temp_lvls_cols) <= 4){
n_col = 1
}else{
n_col = (length(temp_lvls_cols) %/% 4)+1
}
lep = legend(x = x_axp, y = 1, legend = names(temp_lvls_cols),
col = temp_lvls_cols, border = NA,
ncol= n_col, pch = 15, xpd = TRUE, xjust = 0, bty = "n",
cex = annotationFontSize, title = names(annotationColor)[i],
title.adj = 0)
x_axp = x_axp + lep$rect$w
}
}
}
n_mut_samps = length(which(colSums(numMat) != 0))
altStat = paste0("Altered in ", n_mut_samps, " (", round(n_mut_samps/totSamps, digits = 4)*100, "%) of ", totSamps, " samples.")
if(showTitle){
if(is.null(titleText)){
title(main = altStat, outer = TRUE, line = -1, cex.main = titleFontSize)
}else{
title(main = titleText, outer = TRUE, line = -1, cex.main = titleFontSize)
}
}
return(invisible(list(oncomatrix = numMat, vc_legend = vc_codes, vc_color = leg_classes)))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.