suppressPackageStartupMessages(library(optparse))
suppressPackageStartupMessages(library(futile.logger))

### Parsing command line ------------------------------------------------------

option_list <- list(
    make_option(c("--rds"), action="store", type="character", default=NULL,
        help="PureCN output RDS file"),
    make_option(c("--callable"), action="store", type="character", default=NULL,
        help="TMB: File parsable by rtracklayer specifying all callable regions, e.g. generated by GATK CallableLoci"),
    make_option(c("--exclude"), action="store", type="character", default=NULL, 
        help=paste("TMB: File parsable by rtracklayer specifying regions to exclude",
         "from mutation burden calculation, e.g. intronic regions")),
    make_option(c("--maxpriorsomatic"), action = "store", type = "double",
        default = formals(PureCN::callMutationBurden)$max.prior.somatic,
        help = "TMB: Can be used to exclude hotspot mutations with high somatic prior probability [default %default]"),
    make_option(c("--keepindels"), action="store_true", default=FALSE, 
        help="TMB: count indels"),
    make_option(c("--signatures"), action="store_true", default=FALSE, 
        help="Attempt the deconstruction of COSMIC signatures (requires deconstructSigs package)"),
    make_option(c("--out"), action="store", type="character", default=NULL,
        help="File name prefix to which results should be written"),
    make_option(c("-v", "--version"), action="store_true", default=FALSE, 
        help="Print PureCN version"),
    make_option(c("-f", "--force"), action="store_true", default=FALSE, 
        help="Overwrite existing files")
)

opt <- parse_args(OptionParser(option_list=option_list))

if (opt$version) {
    message(as.character(packageVersion("PureCN")))
    q(status=1)
}    

# Parse input rds
infileRds <- opt$rds
if (is.null(infileRds)) stop("Need --rds")
infileRds <- normalizePath(infileRds, mustWork=TRUE)

# Parse outdir

# Parse both BED files restricting covered region
callableFile <- opt$callable
callable <- NULL
if (!is.null(callableFile)) {
    suppressPackageStartupMessages(library(rtracklayer))
    callableFile <- normalizePath(callableFile, mustWork=TRUE)
    flog.info("Reading %s...", callableFile)
    callable <- GRanges(import(callableFile))
}

excludeFile <- opt$exclude
exclude <- NULL
if (!is.null(excludeFile)) {
    suppressPackageStartupMessages(library(rtracklayer))
    excludeFile <- normalizePath(excludeFile, mustWork=TRUE)
    flog.info("Reading %s...", excludeFile)
    exclude <- GRanges(import(excludeFile))
}


### Run PureCN ----------------------------------------------------------------

flog.info("Loading PureCN...")
suppressPackageStartupMessages(library(PureCN))

res <- readCurationFile(infileRds)
sampleid <- res$input$sampleid

.getOutPrefix <- function(opt, infile, sampleid) {
    out <- opt[["out"]]
    if (is.null(out)) {
        if (!is.null(infile) && file.exists(infile)) {
            outdir <- dirname(infile)
            prefix <- sampleid
        } else {
            stop("Need --out")
        }    
    } else {
        outdir <- dirname(out)
        prefix <- basename(out)
    }    
    outdir <- normalizePath(outdir, mustWork=TRUE)
    outPrefix <- file.path(outdir, prefix)
}
outPrefix <- .getOutPrefix(opt, infileRds, sampleid)
     
outfileMb <- paste0(outPrefix, "_mutation_burden.csv")


if (!opt$force && file.exists(outfileMb)) {
    stop(outfileMb, " exists. Use --force to overwrite.")
}    

flog.info("Calling mutation burden...")

fun.countMutation <- eval(formals(callMutationBurden)$fun.countMutation)
if (opt$keepindels) fun.countMutation <- function(vcf) width(vcf) >= 1
    
mb <- callMutationBurden(res, callable=callable, exclude=exclude,
        max.prior.somatic = opt$maxpriorsomatic,
        fun.countMutation = fun.countMutation)

write.csv(cbind(Sampleid=sampleid, mb), file=outfileMb, row.names=FALSE, quote=FALSE)

if (opt$signatures && require(deconstructSigs)) {
    x <- predictSomatic(res)
    x$Sampleid <- res$input$sampleid

    flog.info("deconstructSigs package found.")
    s <- x[ x$ML.SOMATIC & x$prior.somatic > 0.1 & !x$FLAGGED ,]

    genome <- genome(res$input$vcf)[1]
    bsg <- NULL
    if (genome != "hg19") { 
        bsgPkg <- paste0("BSgenome.Hsapiens.UCSC.", genome)
        if (!require(bsgPkg, character.only = TRUE)) { 
            flog.fatal("%s not found.", bsgPkg)
            q(status=0)
        }
        bsg <- get(bsgPkg)
    }
    if (nrow(s) >= 10) {
        sigs.input <- mut.to.sigs.input(s, 
                        sample.id = "Sampleid", chr = "chr", pos = "start", 
                        ref = "REF", alt = "ALT", bsg = bsg)

        sigs <- whichSignatures(tumor.ref = sigs.input, 
                                signatures.ref = signatures.cosmic, 
                                contexts.needed = TRUE, 
                                tri.counts.method = "default")

        outfile <- paste0(outPrefix, "_signatures.pdf")
        pdf(outfile)
        plotSignatures(sigs)
        makePie(sigs)
        dev.off()
        outfile <- paste0(outPrefix, "_signatures.csv")
        write.csv(sigs$weights, file = outfile)
    } else {
        flog.warn("Not enough somatic calls to deconstruct signatures.")
    }
}
