cnv_facets.R

#!/usr/bin/env Rscript
#
# Copyright (C) 2018-2019 University of Glasgow
#
# Author: Dario Beraldi <dario.beraldi@glasgow.ac.uk>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
# DEALINGS IN THE SOFTWARE.

suppressMessages(library(parallel))
suppressMessages(library(argparse))
suppressMessages(library(facets))
suppressMessages(library(data.table))
suppressMessages(library(Rsamtools))
suppressMessages(library(ggplot2))
suppressMessages(library(gridExtra))

# -----------------------------------------------------------------------------

VERSION= sprintf('0.16.1; facets=%s', packageVersion('facets'))

docstring<- sprintf('DESCRIPTION \\n\\
Detect somatic copy number variants (CNVs) and estimate purity and ploidy in a\\n\\
tumour sample given a matched normal sample. The core of the analysis is based\\n\\
on the package FACETS (https://github.com/mskcc/facets).\\n\\
\\n\\
INPUT\\n\\
\\n\\
* A bam file for the tumour sample, a bam file for the matched normal sample, and a vcf\\n\\
  file of SNP positions. Bam and vcf files must be sorted and indexed\\n\\
\\n\\
Or:\\n\\
\\n\\
* A csv file of counts (pileup) generated by a previous run of cnv_facets.R\\n\\
  or by the accompanying program snp-pileup\\n\\
\\n\\
EXAMPLE\\n\\
\\n\\
    cnv_facets.R -n normal.bam -t tumour.bam -vcf snps.vcf.gz -o my_cnv\\n\\
\\n\\
See the online documentation for more details and usage.\\n\\
Version %s', VERSION)

epilog<- 'Options --snp-* are used only to generate the pileup file. They are ignored with option --pileup'
                        
parser<- ArgumentParser(description= docstring, formatter_class= 'argparse.RawTextHelpFormatter', epilog= epilog)

parser$add_argument('--out', '-o', help= 'Output prefix for the output files', required= TRUE)

parser$add_argument('--snp-tumour', '-t', help= 'BAM file for tumour sample', required= FALSE)
parser$add_argument('--snp-normal', '-n', help= 'BAM file for normal sample', required= FALSE)

parser$add_argument('--snp-vcf', '-vcf', help= 'VCF file of SNPs where pileup is to be computed', required= FALSE)

def<- 5
parser$add_argument('--snp-mapq', '-mq', help= sprintf(
'Sets the minimum threshold for mapping quality. Default %s', def), required= FALSE, default= def, type= 'integer')

def<- 10
parser$add_argument('--snp-baq', '-bq', help= sprintf(
'Sets the minimum threshold for base quality. Default %s', def), required= FALSE, default= def, type= 'integer')

parser$add_argument('--snp-count-orphans', '-A', help= 
'Do not discard anomalous read pairs', action= 'store_true')

def<- 1
parser$add_argument('--snp-nprocs', '-N', help= sprintf(
'Number of parallel processes to run to prepare the read pileup file.\\n\\
Each chromsome is assigned to a process. Default %s', def), required= FALSE, default= def, type= 'integer')

parser$add_argument('--pileup', '-p', help= 
'Pileup for matched normal (first sample) and tumour\\n\\
(second sample). Not needed if using BAM input. This file\\n\\
is the <prefix>.cvs.gz output of of a previous run of cnv_facet.R', required= FALSE)

def<- c(25, 4000)
parser$add_argument('--depth', '-d', help= sprintf(
'Minimum and maximum depth in normal sample for a position\\n\\
to be considered. Default %s', 
    paste(def, collapse= ' ')), type= 'integer', default= def, nargs= 2)

parser$add_argument('--targets', '-T', help= sprintf(
'BED file of target regions to scan. It may be the target regions\\n\\
from WEX or panel sequencing protocols. It is not required, even\\n\\
for targeted sequencing, but it may improve the results'))

def<- c(25, 150)
parser$add_argument('--cval', '-cv', help= sprintf(
"Critical values for segmentation in pre-processing and\\n\\
processing. Larger values reduce segmentation. [25 150] is\\n\\
facets default based on exome data. For whole genome consider\\n\\
increasing to [25 400] and for targeted sequencing consider\\n\\
reducing them. Default %s", paste(def, collapse= ' ')), type= 'double', default= def, nargs= 2)

def_nbhd<- list('auto', 250)
parser$add_argument('--nbhd-snp', '-snp', help= sprintf(
'If an interval of size nbhd-snp contains more than one SNP,\\n\\
sample a random one.  This sampling reduces the SNP serial\\n\\
correlation. This value should be similar to the median insert\\n\\
size of the libraries. If "auto" and if using paired end BAM\\n\\
input, use the estimated insert size from the normal bam file.\\n\\
Otherwise use 250. Default %s', def_nbhd[[1]]), type= 'character', default= def_nbhd[[1]])

parser$add_argument('--annotation', '-a', help= sprintf(
'Optional annotation file in BED format where the 4th column\\n\\
contains the feature name (e.g. gene name). CNVs will be\\n\\
annotated with an additional INFO/TAG reporting all the\\n\\
overalapping features'), type= 'character', required= FALSE)

def<- 'hg38'
parser$add_argument('--gbuild', '-g', help= sprintf(
'String indicating the reference genome build. Default %s.', def),
    type= 'character', required= FALSE, default= def, 
    choices= c('hg19', 'hg38', 'mm9', 'mm10'))

parser$add_argument('--unmatched', '-u', help= 
'Normal sample is unmatched. If set, heterozygote SNPs are\\n\\
called using tumor reads only and logOR calculations are different', action= 'store_true')

parser$add_argument('--no-cov-plot', '-np', help= 
'Do not produce coverage plots (they can be memory intensive).\\n\\
Other plots remain unaffected', action= 'store_true')

def_rnd<- 'The name of the input file'
parser$add_argument('--rnd-seed', '-s', help= sprintf(
'Seed for random number generator. Default: %s', def_rnd), type= 'character', default= def_rnd)

# NB: argparse v1.1.1+ required for -v option to work.
parser$add_argument("-v", "--version", action= 'version', version= VERSION)

# -----------------------------------------------------------------------------
genomes<- list(
    HG19= c(chr1= 249250621, 
            chr2= 243199373, 
            chr3= 198022430, 
            chr4= 191154276, 
            chr5= 180915260, 
            chr6= 171115067, 
            chr7= 159138663, 
            chr8= 146364022, 
            chr9= 141213431, 
            chr10= 135534747, 
            chr11= 135006516, 
            chr12= 133851895, 
            chr13= 115169878, 
            chr14= 107349540, 
            chr15= 102531392, 
            chr16= 90354753, 
            chr17= 81195210, 
            chr18= 78077248, 
            chr19= 59128983, 
            chr20= 63025520,  
            chr21= 48129895, 
            chr22= 51304566, 
            chrM= 16571, 
            chrX= 155270560, 
            chrY= 59373566),

    HG38= c(chr1= 248956422,
            chr2= 242193529,
            chr3= 198295559,
            chr4= 190214555,
            chr5= 181538259,
            chr6= 170805979,
            chr7= 159345973,
            chr8= 145138636,
            chr9= 138394717,
            chr11= 135086622,
            chr10= 133797422,
            chr12= 133275309,
            chr13= 114364328,
            chr14= 107043718,
            chr15= 101991189,
            chr16= 90338345,
            chr17= 83257441,
            chr18= 80373285,
            chr19= 58617616,
            chr20= 64444167,
            chr21= 46709983,
            chr22= 50818468,
            chrM=  16569,
            chrX= 156040895,
            chrY= 57227415),

    MM9= c(chr1= 197195432,
           chr2= 181748087,
           chr3= 159599783,
           chr4= 155630120,
           chr5= 152537259,
           chr6= 149517037,
           chr7= 152524553,
           chr8= 131738871,
           chr9= 124076172,
           chr10= 129993255,
           chr11= 121843856,
           chr12= 121257530,
           chr13= 120284312,
           chr14= 125194864,
           chr15= 103494974,
           chr16= 98319150,
           chr17= 95272651,
           chr18= 90772031,
           chr19= 61342430,
           chrM= 16299,
           chrX= 166650296,
           chrY= 15902555),

    MM10= c(chr1= 195471971,
            chr2= 182113224,
            chr3= 160039680,
            chr4= 156508116,
            chr5= 151834684,
            chr6= 149736546,
            chr7= 145441459,
            chr8= 129401213,
            chr9= 124595110,
            chr10= 130694993,
            chr11= 122082543,
            chr12= 120129022,
            chr13= 120421639,
            chr14= 124902244,
            chr15= 104043685,
            chr16= 98207768,
            chr17= 94987271,
            chr18= 90702639,
            chr19= 61431566,
            chrM= 16299,
            chrX= 171031299,
            chrY= 91744698)
)
# ----------------------------------------------------------------------------

avg_insert_size<- function(bam, default){
    idx<- data.table(idxstatsBam(bam))[mapped > 0][order(-mapped)]
    if(nrow(idx) > 10){
        idx<- idx[1:10,]
    }
    ins_size<- c()
    n_reads<- c()
    hdr_len<- length(scanBamHeader(bam)[[1]][['text']])
    for(chrom in idx$seqnames){
        cmd<- sprintf('samtools view -q 3 -f 3 -F 3840 -h %s %s | head -n %s | samtools stats | grep ^SN | cut -f 2-',
                      bam, chrom, 200000 + hdr_len)
        stats<- system(cmd, intern= TRUE)
        size<- grep('insert size average', stats, value= TRUE)
        size<- as.numeric(unlist(strsplit(size, '\t'))[2])
        n<- grep('reads properly paired', stats, value= TRUE)
        n<- as.numeric(unlist(strsplit(n, '\t'))[2])
        if(n == 0){
            next 
        } 
        ins_size<- c(ins_size, size)
        n_reads<- c(n_reads, n)
    }
    if(length(n_reads) == 0){
        return(as.numeric(default))
    }
    return(weighted.mean(ins_size, n_reads))
}

exec_snp_pileup<- function(chrom, snp_vcf, output, normal_bam, tumour_bam, mapq, baq, pseudo_snp, keep_orphans){
    # Execute snp-pileup on chromosome `chrom` 
    # Send tmp output to the same directory of the final output so we are sure
    # we can write there. 
    d<- dirname(output)
    chrom_vcf<- file.path(d, paste0(sub('\\.vcf\\.gz$|\\.vcf\\.bgz$', '', basename(snp_vcf)), '.', chrom, '.vcf'))
    chrom_nbam<- file.path(d, paste0(sub('\\.bam', '', basename(normal_bam)), '.', chrom, '.bam'))
    chrom_tbam<- file.path(d, paste0(sub('\\.bam', '', basename(tumour_bam)), '.', chrom, '.bam'))

    if(keep_orphans) {
        orphans<- '--count-orphans'
    } else {
        orphans<- ''
    }
       
    cmd<- c('#!/bin/bash', '\n',
            '\n',
            'set -eo pipefail', '\n',
            '\n',
            'mkfifo', chrom_vcf, '\n',
            'mkfifo', chrom_nbam, '\n',
            'mkfifo', chrom_tbam, '\n',
            '\n',
            'bcftools view --output-type u', snp_vcf, chrom, '>', chrom_vcf, '&\n',
            'pid_bcf=$!', '\n',
            '\n',
            'samtools view -u', normal_bam, chrom, '>', chrom_nbam, '&\n',
            'pid_nbam=$!', '\n',
            '\n',
            'samtools view -u', tumour_bam, chrom, '>', chrom_tbam, '&\n',
            'pid_tbam=$!', '\n',
            '\n',
            'snp-pileup', 
            '--gzip',
            '--pseudo-snps', pseudo_snp,
            '--min-map-quality', mapq,
            '--min-base-quality', baq,
            '--max-depth 10000000',
            '--min-read-counts', '0,0',
            orphans,
            chrom_vcf, output, chrom_nbam, chrom_tbam,
            '\n',
            'set +e', '\n',
            '\n',
            'wait $pid_bcf; exit_code=$?', '\n',
            'if [[ $exit_code != 0 && $exit_code != 141 ]]; then exit $exit_code; fi', '\n',
            '\n',
            'wait $pid_nbam; exit_code=$?', '\n',
            'if [[ $exit_code != 0 && $exit_code != 141 ]]; then exit $exit_code; fi', '\n',
            '\n',
            'wait $pid_tbam; exit_code=$?', '\n',
            'if [[ $exit_code != 0 && $exit_code != 141 ]]; then exit $exit_code; fi', '\n'
       )
    cmd<- paste(cmd, collapse= ' ')
    script<- file.path(d, sprintf('snp_pileup.%s.sh', chrom))
    write(cmd, script)
    status<- system2("/bin/bash", args= script)
    if(status != 0){
        stop(sprintf('\nError in computing snp pileup. Exit code %s from execution of %s\n' , status, script))
    }
    unlink(c(chrom_vcf, chrom_nbam, chrom_tbam))
    return(cmd)
}

exec_snp_pileup_parallel<- function(snp_vcf, output, normal_bam, tumour_bam, mapq, baq, pseudo_snp, nprocs, keep_orphans){
    
    dtm<- format(Sys.time(), "%y%m%d-%H%M%S")
    tmpdir<- tempfile(pattern = paste0(basename(sub('\\.cvs\\.gz$', '', output)), '_', dtm, '_'), tmpdir = dirname(output))
    dir.create(tmpdir)
    
    chroms<- headerTabix(snp_vcf)$seqnames

    cl<- makeCluster(nprocs, type= 'FORK')
    chrom_csv<- parLapply(cl, chroms, function(chrom){
        chrom_csv<- file.path(tmpdir, paste0(chrom, '.csv.gz'))
        cmd<- exec_snp_pileup(chrom= chrom, 
                              snp_vcf= snp_vcf,
                              output= chrom_csv, 
                              normal_bam= normal_bam, 
                              tumour_bam= tumour_bam, 
                              mapq= mapq, 
                              baq= baq, 
                              pseudo_snp= pseudo_snp,
                              keep_orphans= keep_orphans)
        return(chrom_csv)
    })
    stopCluster(cl)
    concat_csv(chrom_csv, output, tmpdir= tmpdir)
    unlink(tmpdir, recursive= TRUE)
}

concat_csv<- function(csv_list, xfile, tmpdir){
    isGzip<- ifelse(grepl('\\.gz$', xfile), TRUE, FALSE)
    if(isGzip == TRUE){
        conn<- file.path(tmpdir, sub('\\.gz', '', basename(xfile))) 
    } else {
        conn<- xfile
    }
    col.names<- TRUE
    for(csv in csv_list){
        options(datatable.fread.input.cmd.message= FALSE)
        fwrite(x= fread(sprintf('gzip -d -c %s', csv)), file= conn, sep= ',', col.names= col.names, row.name= FALSE, quote= FALSE, append= !isTRUE(col.names))
        options(datatable.fread.input.cmd.message= TRUE)
        col.names<- FALSE
    }
    if(isGzip == TRUE){
        system2(c('gzip', '-c', conn), stdout= xfile)
        unlink(conn)
    }
}

readSnpMatrix2<- function(pileup, gbuild){
    xf<- file(pileup, open= 'r')
    if(summary(xf)$class == 'gzfile'){
        conn<- sprintf('gzip -d -c %s', pileup)
    } else {
        conn<- pileup
    }
    close(xf)
    
    options(datatable.fread.input.cmd.message= FALSE)
    rcmat<- fread(conn, select= c('Chromosome', 'Position', 'File1R', 'File1A', 'File2R', 'File2A'))
    options(datatable.fread.input.cmd.message= TRUE)

    setnames(rcmat, c('File1R', 'File1A', 'File2R', 'File2A'),
                    c('NOR.RD', 'NOR.DP', 'TUM.RD', 'TUM.DP'))
    rcmat[, NOR.DP := NOR.DP + NOR.RD]
    rcmat[, TUM.DP := TUM.DP + TUM.RD]

    chr_prefix<- any(rcmat$Chromosome %in% c(paste0('chr', 1:22), 'chrX'))

    # Unfortunately, facets needs numeric chromsomes. X will be converted later by facets
    rcmat[, Chromosome := sub("^chr", "", Chromosome)]
    setcolorder(rcmat, c('Chromosome', 'Position', 'NOR.DP', 'NOR.RD', 'TUM.DP', 'TUM.RD'))

    # We only keep the major chromosomes 
    if(gbuild %in% c("hg19", "hg38")){
        rcmat<- rcmat[Chromosome %in% c(1:22, 'X')]
    } else if(gbuild %in% c("mm9", "mm10")){
        rcmat<- rcmat[Chromosome %in% c(1:19, 'X')]
    } else {
        write(sprintf('Invalid genome build: %s', gbuild), stderr())
        quit(status= 1)
    }
    return(list(pileup= rcmat, chr_prefix= chr_prefix))
}

facetsRecordToVcf<- function(x){
    stopifnot(is.data.table(x))
    stopifnot(nrow(x) == 1)
    # Convert the annotated facets record to a VCF record.
    vcf<- vector(length= 8)
    vcf[1]<- x$chrom
    vcf[2]<- x$start+1 # TODO: Check +1 is correct
    vcf[3]<- x$seg
    vcf[4]<- 'N'
    vcf[5]<- '<CNV>'
    vcf[6]<- '.'
    vcf[7]<- ifelse(x$type == "NEUTR", 'neutral', 'PASS')

    # INFO field: Keep consistent with header
    vcf[8]<- paste0(
        'SVTYPE=', x$type,
        ';SVLEN=', x$end-x$start,
        ';END=', x$end,
        ';NUM_MARK=', x$num.mark,
        ';NHET=', x$nhet,
        ';CNLR_MEDIAN=', ifelse(is.na(x$cnlr.median), '.', round(x$cnlr.median, 3)), 
        ';MAF_R=', ifelse(is.na(x$mafR), '.', round(x$mafR, 3)), 
        ';SEGCLUST=', x$segclust, 
        ';CNLR_MEDIAN_CLUST=', ifelse(is.na(x$cnlr.median.clust), '.', round(x$cnlr.median.clust, 3)), 
        ';MAF_R_CLUST=', ifelse(is.na(x$mafR.clust), '.', round(x$mafR.clust, 3)), 
        ';CF_EM=', ifelse(is.na(x$cf.em), '.', round(x$cf.em, 3)),
        ';TCN_EM=', ifelse(is.na(x$tcn.em), '.', x$tcn.em),
        ';LCN_EM=', ifelse(is.na(x$lcn.em), '.', x$lcn.em),
        ';CNV_ANN=', ifelse(is.na(x$annotation) || is.null(x$annotation) || x$annotation == '', '.', x$annotation))
    return(vcf)
}

getScriptName<- function(){
    opt<- grep('^--file=', commandArgs(trailingOnly = FALSE), value= TRUE)
    name<- basename(sub('^--file=', '', opt))
    return(name)
}

annotate<- function(cnv, bed_file){
    # Annotate data.table cnv with the features in bed_file

    ann<- data.table(read.table(bed_file, comment.char= '#', header= FALSE, sep= '\t', stringsAsFactors= FALSE, na.strings= ""))
    stopifnot(ncol(ann) >= 4)
    ann<- ann[, 1:4]
    setnames(ann, names(ann), c('chrom', 'start', 'end', 'name'))
    ann<- ann[!is.na(name)]

    # Make feature names VCF compliant. From https://samtools.github.io/hts-specs/VCFv4.1.pdf:
    #
    #     String, no white-space, semi-colons, or equals-signs permitted; commas are
    #     permitted only as delimiters for lists of values
    #
    # Convert special characters using URL encoding
    ann[, name := gsub('%', '%25', name, fixed= TRUE)]
    ann[, name := gsub(',', '%2C', name, fixed= TRUE)]
    ann[, name := gsub('=', '%3D', name, fixed= TRUE)]
    ann[, name := gsub(';', '%3B', name, fixed= TRUE)]
    ann[, name := gsub('|', '%7C', name, fixed= TRUE)]
    ann[, name := gsub(' ', '_', name, fixed= TRUE)] # NB: Not URL encoding for spaces!

    ann<- makeGRangesFromDataFrame(ann, seqnames.field= 'chrom', 
        start.field= 'start', end.field= 'end', keep.extra.columns= TRUE,
        starts.in.df.are.0based= TRUE)

    gcnv<- makeGRangesFromDataFrame(cnv, seqnames.field= 'chrom', start.field= 'start', end.field= 'end', keep.extra.columns= TRUE)
    suppressWarnings({
        ovl<- findOverlaps(query= gcnv, subject= ann, ignore.strand= TRUE)
    })
    hits<- data.table(queryHits= ovl@from, subjectHits= ovl@to, feature= ann$name[ovl@to])
    hits<- hits[, list(feature= paste(feature,  collapse= ',')), by= queryHits]
    gcnv<- as.data.table(gcnv)
    gcnv[, annotation := NA]
    gcnv$annotation[hits$queryHits]<- hits$feature
    gcnv$annotation[gcnv$type == 'NEUTR']<- NA
    setnames(gcnv, 'seqnames', 'chrom')
    gcnv[, chrom := as.character(chrom)]

    stopifnot(cnv$chrom == gcnv$chrom)
    stopifnot(cnv$start == gcnv$start)
    stopifnot(cnv$end == gcnv$end)
    cnv[, annotation := gcnv$annotation]
    return(cnv)
}

make_header<- function(gbuild, genomes, is_chrom_prefixed, cmd, extra){
    # extra: Named vector of additional information. E.g. c(purity=0.5, ploidy= 2.1) 
    header<- c(
        '##fileformat=VCFv4.2',
        sprintf('##reference=%s', gbuild)
    )
    # Contigs
    # --------------------------
    chrom_size<- NULL
    if(gbuild == 'hg19') {chrom_size<- genomes$HG19}
    if(gbuild == 'hg38') {chrom_size<- genomes$HG38}
    if(gbuild == 'mm9') {chrom_size<- genomes$MM9}
    if(gbuild == 'mm10') {chrom_size<- genomes$MM10}

    stopifnot( ! is.null(chrom_size) )

    for(i in 1:length(chrom_size)){
        size<- chrom_size[i]
        name<- names(chrom_size)[i]
        if(is_chrom_prefixed != TRUE){
            name<- sub('^chr', '', name)
        }
        header<- c(header, sprintf('##contig=<ID=%s,length=%s>', name, size))
    }
    # --------------------------

    header<- c(header, '##FILTER=<ID=PASS,Description="All filters passed">')
    header<- c(header, '##FILTER=<ID=neutral,Description="Copy number neutral">')
    header<- c(header, '##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">')
    header<- c(header, '##INFO=<ID=SVLEN,Number=1,Type=Integer,Description="Difference in length between REF and ALT alleles">')
    header<- c(header, '##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">')
    header<- c(header, '##INFO=<ID=NUM_MARK,Number=1,Type=Integer,Description="Number of SNPs in the segment">')
    header<- c(header, '##INFO=<ID=NHET,Number=1,Type=Integer,Description="Number of SNPs that are deemed heterozygous">')
    header<- c(header, '##INFO=<ID=CNLR_MEDIAN,Number=1,Type=Float,Description="Median log-ratio (logR) of the segment. logR is defined by the log-ratio of total read depth in the tumor versus that in the normal">')
    header<- c(header, '##INFO=<ID=CNLR_MEDIAN_CLUST,Number=1,Type=Float,Description="Median log-ratio (logR) of the segment cluster. logR is defined by the log-ratio of total read depth in the tumor versus that in the normal">')
    header<- c(header, '##INFO=<ID=MAF_R,Number=1,Type=Float,Description="Log-odds-ratio (logOR) summary for the segment. logOR is defined by the log-odds ratio of the variant allele count in the tumor versus in the normal">')
    header<- c(header, '##INFO=<ID=MAF_R_CLUST,Number=1,Type=Float,Description="Log-odds-ratio (logOR) summary for the segment cluster. logOR is defined by the log-odds ratio of the variant allele count in the tumor versus that in the normal">')
    header<- c(header, '##INFO=<ID=SEGCLUST,Number=1,Type=Integer,Description="Segment cluster to which the segment belongs">')
    header<- c(header, '##INFO=<ID=CF_EM,Number=1,Type=Float,Description="Cellular fraction, fraction of DNA associated with the aberrant genotype. Set to 1 for normal diploid">')
    header<- c(header, '##INFO=<ID=TCN_EM,Number=1,Type=Integer,Description="Total copy number. 2 for normal diploid">')
    header<- c(header, '##INFO=<ID=LCN_EM,Number=1,Type=Integer,Description="Lesser (minor) copy number. 1 for normal diploid">')
    header<- c(header, '##INFO=<ID=CNV_ANN,Number=.,Type=String,Description="Annotation features assigned to this CNV">')
    header<- c(header, '##ALT=<ID=CNV,Description="Copy number variable region">')

    header<- c(header, cmd)
    
    stopifnot(! is.null(names(extra)) )
    for(i in 1:length(extra)){
        key<- names(extra)[i]
        value<- extra[i]
        stopifnot( key != "" )
        header<- c(header, sprintf('##%s=%s', key, value))
    }
    header<- c(header, '#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO')
    return(paste(header, collapse= '\n'))
}

classify_cnv<- function(dat){
    # Classify CNV. See also https://github.com/mskcc/facets/issues/62
    # dat is a data.table modifed in-place
    dat[, type := NA]
    dat[, type := ifelse((tcn.em == 2 & (lcn.em == 1 | is.na(lcn.em))), 'NEUTR', type)]
    dat[, type := ifelse(is.na(type) & tcn.em == 2 & lcn.em == 2, 'DUP', type)]
    dat[, type := ifelse(is.na(type) & tcn.em == 0, 'DEL', type)]
    dat[, type := ifelse(is.na(type) & tcn.em > 2 & (lcn.em > 0 | is.na(lcn.em)), 'DUP', type)]
    dat[, type := ifelse(is.na(type) & tcn.em == 1, 'HEMIZYG', type)]
    dat[, type := ifelse(is.na(type) & tcn.em == 2 & lcn.em == 0, 'LOH', type)]
    dat[, type := ifelse(is.na(type) & tcn.em > 2 & lcn.em == 0, 'DUP-LOH', type)]
    stopifnot(all(!is.na(dat$type))) # Everything has been classified

    return(dat)
}

prep_coverage_data<- function(rcmat){
    # Reformat the rcmat data.table to make it suitable for plotting histogram
    dhist<- melt(data= rcmat, id= 'Position', measure.vars= c('NOR.DP', 'TUM.DP'),
        variable.name= 'sample', value.name= 'depth')
    dhist[, Position := NULL]
    dhist[, sample := ifelse(sample == 'NOR.DP', 'Normal', 'Tumour')]
    dhist<- dhist[depth > 0]
    dhist<- dhist[, list(depth= ifelse(.SD$depth > quantile(.SD$depth, 0.99), quantile(.SD$depth, 0.99), .SD$depth)), by= sample]
    nsites<- dhist[, .N, by= sample]
    dhist<- merge(dhist, nsites)
    dhist[, label := paste(sample, 'N=', N)]
    dhist[, sample := NULL]
    return(dhist)
}

plot_coverage<- function(rcmat, rcmat_flt, fname, title){
    
    xall<- ggplot(data= prep_coverage_data(rcmat), aes(x= depth)) +
        geom_histogram(bins= 20, colour= 'white', fill= '#F8766D') +
        facet_wrap(~ label) +
        xlab('All positions [capped]') + 
        theme(axis.title.x = element_text(colour = "#F8766D"))
        
    xflt<- ggplot(data= prep_coverage_data(rcmat_flt), aes(x= depth)) +
        geom_histogram(bins= 20, colour= 'white', fill= '#00BFC4') +
        facet_wrap(~ label) + 
        xlab('Filtered positions [capped]') + 
        theme(axis.title.x = element_text(colour = "#00BFC4"))

    pdf(NULL) # Prevent Rplots.pdf to be generated
    gg<- arrangeGrob(xall, xflt, top= title)
    ggsave(fname, gg, width= 16, height= 18, units= 'cm')
    rm(xall, xflt, gg)
    x_<- gc(verbose= FALSE) 
}

filter_rcmat<- function(rcmat, min_ndepth, max_ndepth, target_bed){
    rcmat_flt<- rcmat[NOR.DP >= min_ndepth & NOR.DP < max_ndepth]
    if(is.null(target_bed) || nrow(rcmat_flt) == 0){
        return(rcmat_flt)
    }
    targets<- makeGRangesFromDataFrame(target_bed, seqnames.field= 'V1', 
        start.field= 'V2', end.field= 'V3', starts.in.df.are.0based= TRUE)

    # Convert rcmat to GRanges
    rcmat_flt<- makeGRangesFromDataFrame(rcmat_flt, seqnames.field= 'Chromosome', 
        start.field= 'Position', end.field= 'Position', keep.extra.columns= TRUE)
    
    hits<- findOverlaps(query= rcmat_flt, subject= targets, ignore.strand= TRUE)
    
    rcmat_flt<- as.data.table(rcmat_flt[unique(hits@from)])
    rcmat_flt[, start := NULL]
    rcmat_flt[, width := NULL]
    rcmat_flt[, strand := NULL]
    setnames(rcmat_flt, c('seqnames', 'end'), c('Chromosome', 'Position'))
    # This is to print (i.e., for debugging) the data.table obj without calling
    # print twice. See 2.23 at
    # https://cran.r-project.org/web/packages/data.table/vignettes/datatable-faq.html#why-do-i-have-to-type-dt-sometimes-twice-after-using-to-print-the-result-to-console
    rcmat_flt[]
    return(rcmat_flt)
}

run_facets<- function(
            pre_rcmat,
            pre_gbuild,
            pre_snp.nbhd,
            pre_het.thresh,
            pre_cval,
            pre_deltaCN,
            pre_unmatched,
            pre_ndepth,
            pre_ndepthmax,

            proc_cval,
            proc_min.nhet,
            proc_dipLogR,
            
            emcncf_unif,
            emcncf_min.nhet,
            emcncf_maxiter,
            emcncf_eps
    ){
    # Run the core functions of facets for segmentation, purity etc.
    # Here is where the actual CNV discovery happen.
    # Param prefix matches the facets function they go to.
    write(sprintf('[%s] Preprocessing sample...', Sys.time()), stderr())
    
    xx<- preProcSample(
           rcmat=      pre_rcmat,
           gbuild=     pre_gbuild, 
           snp.nbhd=   pre_snp.nbhd, 
           het.thresh= pre_het.thresh, 
           cval=       pre_cval, 
           deltaCN=    pre_deltaCN, 
           unmatched=  pre_unmatched, 
           ndepth=     pre_ndepth,    
           ndepthmax=  pre_ndepthmax
    )
    rm(pre_rcmat)
    x_ <- gc(verbose= FALSE)

    write(sprintf('[%s] Processing sample...', Sys.time()), stderr())
    proc_out<- procSample(xx, 
                          cval=     proc_cval, 
                          min.nhet= proc_min.nhet, 
                          dipLogR=  proc_dipLogR)
    proc_out$jointseg<- data.table(proc_out$jointseg)
    proc_out$out<- data.table(proc_out$out)

    write(sprintf('[%s] Fitting model...', Sys.time()), stderr())
    emcncf_fit<- emcncf(x=        proc_out, 
                        unif=     emcncf_unif, 
                        min.nhet= emcncf_min.nhet, 
                        maxiter=  emcncf_maxiter, 
                        eps=      emcncf_eps)
    names(emcncf_fit$purity)<- NULL
    names(emcncf_fit$ploidy)<- NULL
    emcncf_fit[['cncf']]<- data.table(emcncf_fit[['cncf']])[order(chrom, start)]
    return(list(proc_out= proc_out, emcncf_fit= emcncf_fit))
}

reset_chroms<- function(cncf, gbuild, chr_prefix){
    # Reset chomosome names *in-place*
    # Reset chrom X
    if(gbuild %in% c('hg19', 'hg38')){
        stopifnot(length(unique(cncf$chrom)) <= 23)
        cncf[, chrom := ifelse(chrom == 23, 'X', chrom)]
    } else if(gbuild %in% c('mm9', 'mm10')){
        stopifnot(length(unique(cncf$chrom)) <= 20)
        cncf[, chrom := ifelse(chrom == 20, 'X', chrom)]
    } else {
        stop(sprintf('Invalid genome build: %s', gbuild))
    }
    # Reset chrom names
    stopifnot(is.logical(chr_prefix))
    if(chr_prefix == TRUE){
        cncf[, chrom := paste0('chr', chrom)]
    }
}

if(sys.nframe() == 0){
    # Script is being executed from the command line

    # ---------------- [Validate arguments] -------------

    xargs<- parser$parse_args()

    if(is.null(xargs$pileup) && (is.null(xargs$snp_tumour) || is.null(xargs$snp_normal))){
        write('Please provide a pileup file or both a tumour and a normal bam file', stderr()) 
        quit(status= 1)
    }
    if(! is.null(xargs$snp_tumour) && xargs$snp_tumour == xargs$snp_normal){
        write('Input bam files for tumour and normal cannot be the same', stderr())
        quit(status= 1)
    }

    if(is.null(xargs$pileup) && is.null(xargs$snp_vcf)){
        write('Please use the --snp_vcf option to provide a VCF file of SNPs', stderr()) 
        quit(status= 1)
    }
    if(xargs$depth[1] > xargs$depth[2]){
        write('Error in argument --depth: min > max', stderr()) 
        quit(status= 1)
    }
    if(xargs$cval[1] > xargs$cval[2]){
        write('Error in argument --cval. Value for pre-processing (first argument)\nmust be lower than the value for processing (second argument)', stderr()) 
        quit(status= 1)
    }
    
    est_insert_size<- NA
    if(xargs$nbhd_snp == "auto"){
        if(!is.null(xargs$pileup)){
            nbhd_snp<- def_nbhd[[2]]
        } else {
            nbhd_snp<- avg_insert_size(xargs$snp_normal, default= def_nbhd[[2]])
            est_insert_size<- round(nbhd_snp, 2)
        }
    } else if(is.na(as.integer(xargs$nbhd_snp))){
        write('Argument to --nbhd-snp must be "auto" or an integer', stderr())
        quit(status= 1)
    } else {
        nbhd_snp<- as.integer(xargs$nbhd_snp)  
    }
    nbhd_snp<- as.integer(nbhd_snp)

    # --------------- [Start processing] ----------------

    dir.create(dirname(xargs$out), recursive= TRUE, showWarnings= FALSE)

    if(xargs$rnd_seed == def_rnd){
        if( ! is.null(xargs$pileup) ){
            seed<- sum(utf8ToInt(xargs$pileup))
        } else {
            seed<- sum(utf8ToInt(paste(xargs$snp_normal, xargs$snp_tumour)))
        }
    } else if( ! is.na(suppressWarnings(as.numeric(xargs$rnd_seed)))){
        seed<- as.numeric(xargs$rnd_seed)
    } else {
        seed<- sum(utf8ToInt(xargs$rnd_seed))
    }
    set.seed(seed)

    if(is.null(xargs$pileup)){
        pileup<- paste0(xargs$out, '.csv.gz')
        basename(xargs$snp_tumour)
        exec_snp_pileup_parallel(snp_vcf= xargs$snp_vcf, 
                                 output= pileup, 
                                 normal_bam= xargs$snp_normal, 
                                 tumour_bam= xargs$snp_tumour, 
                                 mapq= xargs$snp_mapq, 
                                 baq= xargs$snp_baq, 
                                 pseudo_snp= nbhd_snp,
                                 nprocs= ceiling(xargs$snp_nprocs / 2),
                                 keep_orphans= xargs$snp_count_orphans
                                 )
    } else {
        pileup<- xargs$pileup
    }

    write(sprintf('[%s] Loading file %s...', Sys.time(), pileup), stderr())
    rcmat<- readSnpMatrix2(pileup, xargs$gbuild)
    if(nrow(rcmat$pileup) == 0){
        write(sprintf('[%s] Analysis stopped: pileup file %s has no records', Sys.time(), pileup), stderr())
        quit(status= 1)
    }

    # Read targets file
    if(!is.null(xargs$targets)){
        targets<- data.table(read.table(xargs$targets, comment.char= '#', header= FALSE, sep= '\t'))
        if(rcmat[['chr_prefix']]){
            # Original input has chr prefix, so we need to strip it from targets DF
            n<- length(unique(targets$V1))
            targets[, V1 := sub("^chr", "", V1)]
            stopifnot(n == length(unique(targets$V1)))
        }
    } else {
        targets<- NULL
    }
    rcmat_flt<- filter_rcmat(rcmat= rcmat[['pileup']], min_ndepth= xargs$depth[1], 
        max_ndepth= xargs$depth[2], target_bed= targets)
   
    if(xargs$no_cov_plot == FALSE){
        write(sprintf('[%s] Plotting histogram of coverage...', Sys.time()), stderr())
        plot_coverage(rcmat= rcmat[['pileup']], rcmat_flt= rcmat_flt, 
            fname= paste0(xargs$out, '.cov.pdf'), title= paste0('Depth of coverage\n', xargs$out))
    }

    rcmat[['pileup']]<- NULL
    
    x_ <- gc(verbose= FALSE)
    
    # ------------- [Run FACETS] --------------
    facets<- run_facets(
           pre_rcmat=       rcmat_flt,
           pre_gbuild=      xargs$gbuild, 
           pre_snp.nbhd=    nbhd_snp, 
           pre_het.thresh=  ifelse(xargs$unmatched, 0.1, 0.25), 
           pre_cval=        xargs$cval[1], 
           pre_deltaCN=     0, 
           pre_unmatched=   xargs$unmatched,
           pre_ndepth=      1,   # We subset the input matrix instead of using ndepth*
           pre_ndepthmax=   1e8, # options
           proc_cval=       xargs$cval[2], 
           proc_min.nhet=   15, 
           proc_dipLogR=    NULL,
           emcncf_unif=     FALSE, 
           emcncf_min.nhet= 15,
           emcncf_maxiter=  20,
           emcncf_eps=      1e-3
    ) 
    # -----------------------------------------

    write(sprintf('[%s] Writing output', Sys.time()), stderr())
    cncf<- copy(facets$emcncf_fit$cncf)

    reset_chroms(cncf= cncf, gbuild= xargs$gbuild, chr_prefix= rcmat$chr_prefix)
    classify_cnv(cncf)

    if(is.null(xargs$annotation) == FALSE){
        annotate(cncf, xargs$annotation)
    }

    vcf_tmp<- tempfile(pattern= paste0(basename(xargs$out), '.'), tmpdir= dirname(xargs$out), fileext= '.vcf')
    cmd<- sprintf('##%sCommand=%s; Version=%s; Date=%s', getScriptName(), paste(commandArgs(), collapse= ' '), VERSION, Sys.time())
    
    extra<- c(purity= facets$emcncf_fit$purity, 
              ploidy= facets$emcncf_fit$ploidy, 
              dipLogR= facets$emcncf_fit$dipLogR, 
              emflags= facets$emcncf_fit$emflags,
              est_insert_size= est_insert_size)

    header<- make_header(gbuild= xargs$gbuild, 
                         genomes= genomes, 
                         is_chrom_prefixed= rcmat$chr_prefix, 
                         cmd= cmd, 
                         extra= extra)
    write(header, vcf_tmp)

    for(i in 1:nrow(cncf)){
        write(paste(facetsRecordToVcf(cncf[i]), collapse= '\t'), vcf_tmp, append= TRUE)
    }

    x_<- bgzip(vcf_tmp, dest= paste0(xargs$out, '.vcf.gz'), overwrite= TRUE)
    x_<- indexTabix(paste0(xargs$out, '.vcf.gz'), format= 'vcf4')
    unlink(vcf_tmp)

    write(sprintf('[%s] Plotting genome...', Sys.time()), stderr())
    png(paste0(xargs$out, '.cnv.png'), units="px", width=1600, height=1600, res=300)
    sname<- sprintf('%s; ploidy= %.2f; purity= %.2f', basename(xargs$out), facets$emcncf_fit$ploidy, facets$emcncf_fit$purity)
     plotSample(x= facets$proc_out, emfit= facets$emcncf_fit, sname= sname)
    x_ <- dev.off()

    write(sprintf('[%s] Plotting spider...', Sys.time()), stderr())
    pdf(paste0(xargs$out, '.spider.pdf'), width= 16/2.54, height= 14/2.54)
    par(las= 1, mar= c(3, 3, 1, 1), mgp= c(1.5, 0.5, 0), tcl= -0.3)
    logRlogORspider(facets$proc_out$out, facets$proc_out$dipLogR)
    x_ <- dev.off()
    
    si<- capture.output(sessionInfo())
    write(si, stderr())
}