SomaticSignatures.R

#!/usr/bin/env Rscript
args <- commandArgs()
scriptPath = strsplit(args[4], '=', fixed = T)[[1]][2]
scriptName = basename(scriptPath)
scriptDir = dirname(scriptPath)
args = args[-(1:5)]
source(paste0(scriptDir, '/common.R'))

usage = function(){
    cat(paste0("Usage: ", scriptName) )
    cat(" -p=outputName.pdf input1.vcf input2.vcf [input3.vcf]... >motifMatrix.tsv
Option:
  Input:
    -r|ref          FILE       The reference genome in FASTA or 2bit format
    -m|mutect                  The input vcf is generated by MuTect1
  Output:
    -p|pdf          PDF        The output figure in pdf[figure.pdf]
    -c|compare      TSV        Comparison to knwon signature21
    -w|width        INT        The pdf width
       height       INT        The pdf height
  Efficiency:
  Tool:
       pca                     Use PCA instead of NMF
    -s|sigN         INT        The signature number[5]
    -a|access       INTs       The comma-separated signature numbers for assessment[2,3,4,5,6,7,8]
  Other:
       sample       STRs       The comma-separated sample names
    -h|help                    Show help
")
    q(save = 'no')
}

myPdf = 'figure.pdf'
access = 2:8
sigN = 5

if(length(args) >= 1){
    for(i in 1:length(args)){
        arg = args[i]
        tmp = parseArg(arg, 'r(ef)?', 'r')
        if(!is.null(tmp)){
            ref = tmp
            args[i] = NA
            next
        }
        if(arg == '-m' || arg == '-mutect'){
            mutect = TRUE
            args[i] = NA
            next
        }
        tmp = parseArg(arg, 'p(df)?', 'p')
        if(!is.null(tmp)){
            myPdf = tmp
            args[i] = NA
            next
        }
        tmp = parseArg(arg, 'c(omparison)?', 'c')
        if(!is.null(tmp)){
            cmp = tmp
            args[i] = NA
            next
        }
        tmp = parseArgAsNum(arg, 'w(idth)?', 'w')
        if(!is.null(tmp)){
            myWidth = tmp
            args[i] = NA
            next
        }
        tmp = parseArgAsNum(arg, 'height', 'height')
        if(!is.null(tmp)){
            height = tmp
            args[i] = NA
            next
        }
        if(arg == '-pca'){
            PCA = TRUE
            args[i] = NA
            next
        }
        tmp = parseArgAsNum(arg, 's(igN)?', 's')
        if(!is.null(tmp)){
            sigN = tmp
            args[i] = NA
            next
        }
        tmp = parseArgNums(arg, 'a(ccess)?', 'a')
        if(!is.null(tmp)){
            access = tmp
            args[i] = NA
            next
        }
        tmp = parseArgStrs(arg, 'sample', 'sample')
        if(!is.null(tmp)){
            mySampleNames = tmp
            args[i] = NA
            next
        }
        if(arg == '-h' || arg == '-help') usage()
    }
}
args = args[!is.na(args)]


suppressPackageStartupMessages(library(SomaticSignatures))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(ggdendro))


myCmd = 'pdf(myPdf'
if(exists('myWidth')) myCmd = paste0(myCmd, ', width = myWidth')
if(exists('height')) myCmd = paste0(myCmd, ', height = height')
myCmd = paste0(myCmd, ')')
eval(parse(text = myCmd))


if(grepl('.2bit$', ref)){
    refGenome = TwoBitFile(ref)
}else{
    refGenome = FaFile(ref)
}

sink(stderr())
cat(paste0('[DEBUG] ', Sys.time(), ' Check if the following variables are correct as expected:\n'))
cat(paste0('\nref\t', ref))
cat('\nmuTect\t'); if(exists('muTect')) cat(muTect)
cat('\nmyPdf\t'); if(exists('myPdf')) cat(myPdf)
cat('\nwidth\t'); if(exists('myWidth')) cat(myWidth)
cat('\nheight\t'); if(exists('height')) cat(height)
cat('\npca\t'); if(exists('PCA')) cat(PCA)
cat('\nsigN\t'); if(exists('sigN')) cat(sigN)
cat('\naccess\t'); if(exists('access')) cat(access)
cat('\nsample\t'); if(exists('mySampleNames')) cat(mySampleNames)
cat('\nInput vcf\t'); cat(args)
cat('\n')

cat(paste0('[INFO] ', Sys.time(), ' Start fetching context...\n'))
for(i in 1:length(args)){
    if(exists('mutect')){
        VR = readMutect(args[i])
    }else{
        cat(paste0('Reading file ', args[i]), '\n')
        VR = readVcfAsVRanges(args[i], 'myGenome')
        softFilterMatrix(VR) = FilterMatrix(matrix(nrow=length(VR),ncol=0),FilterRules())
    }
    if(exists('mySampleNames')){
        VR$sampleNames = rep(mySampleNames[i], length(VR))
    }
    cat(paste0('Get motif for ', args[i]), '\n')
    motif = mutationContext(VR, refGenome)
    if(exists('motifs')){
        motifs = c(motifs, motif)
    }else{
        motifs = motif
    }
}
cat(paste0('[INFO] ', Sys.time(), ' Finish fetching context...\n'))


motifMatrix = motifMatrix(motifs, normalize = TRUE)

if(exists('cmp')){
    data(signatures21)
    library(reshape2)
    cmpResult = melt(cor(motifMatrix, signatures21))
    dist = apply(signatures21, 2, function(x){
        apply(motifMatrix, 2, function(v1, v2 = x){
            dist(matrix(c(v1, v2), byrow = T, nrow = 2))
        })
    })
    dist = melt(dist)
    colnames(cmpResult) = c('Sample', 'Signature', 'Pearson.Correlation')
    cmpResult=cbind(cmpResult, Euclidean.Distance=dist$value)
    write.table(cmpResult, cmp, quote = F, sep = '\t', row.names = F)
}

cat(paste0('[INFO] ', Sys.time(), ' Start assessment...\n'))
if(exists('PCA')){
    gof = assessNumberSignatures(motifMatrix, access, pcaDecomposition)
    cat(paste0('[INFO] ', Sys.time(), ' Finish assessment...\n'))
    cat(paste0('[INFO] ', Sys.time(), ' Start identify signatures...\n'))
    sigs = identifySignatures(motifMatrix, sigN, pcaDecomposition)
    method = 'PCA'
}else{
    gof = assessNumberSignatures(motifMatrix, access, nReplicates = 5)
    cat(paste0('[INFO] ', Sys.time(), ' Finish assessment...\n'))
    cat(paste0('[INFO] ', Sys.time(), ' Start identify signatures...\n'))
    sigs = identifySignatures(motifMatrix, sigN, nmfDecomposition)
    method = 'NMF'
}
cat(paste0('[INFO] ', Sys.time(), ' Finish identify signatures...\n'))


plotNumberSignatures(gof)
plotObservedSpectrum(sigs) + ggtitle("Observed Spectrum")
plotFittedSpectrum(sigs) + ggtitle("Fitted Spectrum")
plotSignatureMap(sigs) + ggtitle(paste0("Somatic Signatures: ", method, " - Heatmap"))
plotSignatures(sigs) + ggtitle(paste0("Somatic Signatures: ", method, " - Barchart"))
plotSampleMap(sigs) + ggtitle(paste0("Sample Map: ", method, " - Heatmap"))
plotSamples(sigs) + ggtitle(paste0("Signature Contribution: ", method, " - Barchart"))


clu_motif = clusterSpectrum(motifMatrix, "motif")
p = ggdendrogram(clu_motif, rotate = TRUE)
p

sink()
write.table(cbind(Motif=row.names(motifMatrix),motifMatrix), stdout(), quote = F, sep = '\t', row.names = F)