other/pipeline/trunk/mutation_signature.R

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

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