extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/grid/scripts/mutation_signature.R

Code
Comments
Other
Rev Date Author Line
6656 29 Mar 22 nicklas 1 #!/usr/bin/env Rscript
6656 29 Mar 22 nicklas 2
6656 29 Mar 22 nicklas 3 # $Id: mutation_signature.R 5807 2020-01-10 07:18:28Z nicklas $
6656 29 Mar 22 nicklas 4 # A simple helper script for calculating mutation signature scores
6656 29 Mar 22 nicklas 5 # based on COSMIC_Cancer_signatures_probabilities.RData
6656 29 Mar 22 nicklas 6 #
6656 29 Mar 22 nicklas 7 # The script need 3 or 4 parameters
6656 29 Mar 22 nicklas 8 # 1: Name of sample (used as table header and in the plot)
6656 29 Mar 22 nicklas 9 # 2: Path to VCF file with mutations
6656 29 Mar 22 nicklas 10 # 3: Path to COSMIC_Cancer_signatures_probabilities.RData
6656 29 Mar 22 nicklas 11 # 4: (Optional) Path to output pdf file
6656 29 Mar 22 nicklas 12 #
6656 29 Mar 22 nicklas 13 # The script will output scores for each of the 30 signatures as a
6656 29 Mar 22 nicklas 14 # tab-separated table to stdout
6656 29 Mar 22 nicklas 15 #
6656 29 Mar 22 nicklas 16
6656 29 Mar 22 nicklas 17 ## Get command line arguments
6656 29 Mar 22 nicklas 18 args = commandArgs(trailingOnly=TRUE)
6656 29 Mar 22 nicklas 19 if (length(args) < 3) stop('At least 3 parameters are required')
6656 29 Mar 22 nicklas 20
6656 29 Mar 22 nicklas 21 sampleName <- args[1]
6656 29 Mar 22 nicklas 22 vcfPath <- args[2]
6656 29 Mar 22 nicklas 23 cosmicSignature <- args[3]
6656 29 Mar 22 nicklas 24
6656 29 Mar 22 nicklas 25 if (!file.exists(vcfPath)) stop(paste('Can\'t find VCF file: ', vcfPath, sep=''))
6656 29 Mar 22 nicklas 26 if (!file.exists(cosmicSignature)) stop(paste('Can\'t find COSMIC signature file: ', cosmicSignature, sep=''))
6656 29 Mar 22 nicklas 27
6656 29 Mar 22 nicklas 28 plot <- length(args) >= 4
6656 29 Mar 22 nicklas 29 if (plot)
6656 29 Mar 22 nicklas 30 {
6656 29 Mar 22 nicklas 31   pdfPath = args[4]
6656 29 Mar 22 nicklas 32   if (!dir.exists(dirname(pdfPath))) stop(paste('Can`t find directory for PDF file: ', pdfPath, sep=''))
6656 29 Mar 22 nicklas 33 }
6656 29 Mar 22 nicklas 34
6656 29 Mar 22 nicklas 35 ## Other configuration options
6656 29 Mar 22 nicklas 36 ref_genome <- "BSgenome.Hsapiens.UCSC.hg38"
6656 29 Mar 22 nicklas 37
6656 29 Mar 22 nicklas 38 ## Load libraries
6656 29 Mar 22 nicklas 39 library(GenomeInfoDbData)
6656 29 Mar 22 nicklas 40 library(MutationalPatterns)
6656 29 Mar 22 nicklas 41 library(BSgenome)
6656 29 Mar 22 nicklas 42 library(ref_genome, character.only = TRUE)
6656 29 Mar 22 nicklas 43
6656 29 Mar 22 nicklas 44 ## Load the mutation signatures
6656 29 Mar 22 nicklas 45 load(cosmicSignature)
6656 29 Mar 22 nicklas 46 ## reduce object matrix to signature columns only
6656 29 Mar 22 nicklas 47 cancer_signatures = as.matrix(cancer_signatures[,4:33])
6656 29 Mar 22 nicklas 48
6656 29 Mar 22 nicklas 49 ## read vcf 
6656 29 Mar 22 nicklas 50 vcfs <- read_vcfs_as_granges(vcfPath, sampleName, ref_genome)
6656 29 Mar 22 nicklas 51
6656 29 Mar 22 nicklas 52 ## extract mutations
6656 29 Mar 22 nicklas 53 mut_mat <- mut_matrix(vcf_list = vcfs, ref_genome = ref_genome)
6656 29 Mar 22 nicklas 54
6656 29 Mar 22 nicklas 55 ## fit mutations to cancer signatures
6656 29 Mar 22 nicklas 56 fit_res <- fit_to_signatures(mut_mat, cancer_signatures)
6656 29 Mar 22 nicklas 57
6656 29 Mar 22 nicklas 58 ## get contribution
6656 29 Mar 22 nicklas 59 x <- fit_res$contribution
6656 29 Mar 22 nicklas 60
6656 29 Mar 22 nicklas 61 ## plot as pdf?
6656 29 Mar 22 nicklas 62 if (plot)
6656 29 Mar 22 nicklas 63 {
6656 29 Mar 22 nicklas 64   signature.color<-c(
6656 29 Mar 22 nicklas 65       rgb(85/255,46/255,141/255), #s1
6656 29 Mar 22 nicklas 66       rgb(240/255,105/255,160/255), #s2
6656 29 Mar 22 nicklas 67       rgb(7/255,83/255,157/255),#s3
6656 29 Mar 22 nicklas 68       "lightgray",#s4
6656 29 Mar 22 nicklas 69       rgb(140/255,108/255,175/255),#s5
6656 29 Mar 22 nicklas 70       rgb(10/255,110/255,56/255),#s6
6656 29 Mar 22 nicklas 71       "darkgray",#s7
6656 29 Mar 22 nicklas 72       rgb(68/255,147/255,198/255),#s8
6656 29 Mar 22 nicklas 73       "beige",#s9
6656 29 Mar 22 nicklas 74       "darkslategray",#s10
6656 29 Mar 22 nicklas 75       "goldenrod2",#s11
6656 29 Mar 22 nicklas 76       "firebrick4",#s12
6656 29 Mar 22 nicklas 77       rgb(250/255,197/255,192/255),#s13
6656 29 Mar 22 nicklas 78       "azure2",#s14
6656 29 Mar 22 nicklas 79       "burlywood3",#s15
6656 29 Mar 22 nicklas 80       "darksalmon",#s16
6656 29 Mar 22 nicklas 81       rgb(218/255,218/255,234/255),#s17
6656 29 Mar 22 nicklas 82       rgb(158/255,154/255,200/255),#s18
6656 29 Mar 22 nicklas 83       "darkorchid4",#s19
6656 29 Mar 22 nicklas 84       rgb(66/255,171/255,93/255),#s20
6656 29 Mar 22 nicklas 85       "orange",#s21
6656 29 Mar 22 nicklas 86       "darkred",#s22
6656 29 Mar 22 nicklas 87       "black",#s23
6656 29 Mar 22 nicklas 88       "chocolate",#s24
6656 29 Mar 22 nicklas 89       "cornflowerblue",#s25
6656 29 Mar 22 nicklas 90       "cadetblue4",#s26
6656 29 Mar 22 nicklas 91       "aquamarine4",#s27
6656 29 Mar 22 nicklas 92       "aquamarine",#s28
6656 29 Mar 22 nicklas 93       "antiquewhite",#s29
6656 29 Mar 22 nicklas 94       "yellow" #s30
6656 29 Mar 22 nicklas 95   )
6656 29 Mar 22 nicklas 96   
6656 29 Mar 22 nicklas 97   names(signature.color) <- paste("Signature.",1:30,sep="")
6656 29 Mar 22 nicklas 98   
6656 29 Mar 22 nicklas 99   main.label.file <- basename(vcfPath)
6656 29 Mar 22 nicklas 100   
6656 29 Mar 22 nicklas 101   ####### Plot and calculate proportions
6656 29 Mar 22 nicklas 102   pdf(pdfPath)
6656 29 Mar 22 nicklas 103   
6656 29 Mar 22 nicklas 104   prop.x <- prop.table(x, 2)
6656 29 Mar 22 nicklas 105   ii <- match(rownames(prop.x), names(signature.color))
6656 29 Mar 22 nicklas 106   signature.color <- signature.color[ii]
6656 29 Mar 22 nicklas 107   
6656 29 Mar 22 nicklas 108   #select signatures with some contribution
6656 29 Mar 22 nicklas 109   select <- which(rowSums(prop.x) > 0)
6656 29 Mar 22 nicklas 110   
6656 29 Mar 22 nicklas 111   ## plot
6656 29 Mar 22 nicklas 112   my.xlim <- c(0, ncol(prop.x))
6656 29 Mar 22 nicklas 113   if (ncol(prop.x)==1)
6656 29 Mar 22 nicklas 114   {
6656 29 Mar 22 nicklas 115     my.xlim <- c(0, 2) #increases margins on right side
6656 29 Mar 22 nicklas 116     prop.x <- as.matrix(prop.x[select,]) 
6656 29 Mar 22 nicklas 117     signature.color <- signature.color[select]
6656 29 Mar 22 nicklas 118   }
6656 29 Mar 22 nicklas 119   barplot(prop.x*100,
6656 29 Mar 22 nicklas 120       xlim=my.xlim,
6656 29 Mar 22 nicklas 121       ylab="%",
6656 29 Mar 22 nicklas 122       main=sampleName,
6656 29 Mar 22 nicklas 123       cex.main=0.7,
6656 29 Mar 22 nicklas 124       las=2,
6656 29 Mar 22 nicklas 125       cex.names = 0.6,
6656 29 Mar 22 nicklas 126       col=signature.color,
6656 29 Mar 22 nicklas 127       legend.text=rownames(prop.x),
6656 29 Mar 22 nicklas 128       args.legend=list("top", bg="white", ncol=1, cex=0.7)
6656 29 Mar 22 nicklas 129   )
6656 29 Mar 22 nicklas 130   
6656 29 Mar 22 nicklas 131   my.at <- as.numeric(filter(x=c(0,cumsum(round(prop.x[,1], 3)*100)), filter=rep(1, 2)) / 2)[1:length(prop.x[,1])]
6656 29 Mar 22 nicklas 132   my.sign <- paste("s.", sapply(row.names(prop.x), function(x) strsplit(x, "\\.")[[1]][2]), sep="")
6656 29 Mar 22 nicklas 133   my.pcr <- paste(round(prop.x[,1], 3)*100, "%", sep="")
6656 29 Mar 22 nicklas 134   my.labels <- paste(my.pcr, " (", my.sign,")", sep="")
6656 29 Mar 22 nicklas 135   axis(4, at=my.at, labels=my.labels, las=2, cex.axis=0.7, line=-12, lwd=0)
6656 29 Mar 22 nicklas 136   
6656 29 Mar 22 nicklas 137   mtext(paste("signature proportions:", main.label.file), line=0.5, cex=0.6)
6656 29 Mar 22 nicklas 138   foo <- dev.off()
6656 29 Mar 22 nicklas 139   
6656 29 Mar 22 nicklas 140 }# end if plot
6656 29 Mar 22 nicklas 141
6656 29 Mar 22 nicklas 142 ## print table to stdout
6656 29 Mar 22 nicklas 143 write.table(x, file=stdout(), row.names=TRUE, col.names=FALSE, quote=FALSE, sep="\t")  

Generated on Tue Nov 21 03:33:42 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)