metrics.Rmd

---
title: "Metrics"
output: html_document
---

# in silico mixture results

Perform benchmarking on in silico mixture results.

Contains code to regenerate figures 2b-d, 3, 4 and supplementary figures 8-9 from the SVclone paper.

```{r setup, include=FALSE}
knitr::opts_chunk$set(autodep        = TRUE,
                      cache          = TRUE,
                      cache.path     = "cache/",
                      cache.comments = TRUE,
                      echo           = FALSE,
                      error          = FALSE,
                      fig.path       = "figures/",
                      fig.width      = 18,
                      fig.height     = 12,
                      dev            = c('png', 'pdf'),
                      message        = FALSE,
                      warning        = FALSE)
```


```{r libraries}
library(ccube)
library(dplyr)
library(colorspace)
library(ggplot2)
library(pROC)
```


```{r options}
options(stringsAsFactors = F)
```


```{r source}
source('R/metrics_funcs.R')
source('R/util.R')
source('R/load.R')
```


```{r funcs_and_vars}
gg_color_hue <- function(n) {
  hues = seq(15, 375, length=n+1)
  hcl(h=hues, l=65, c=100)[1:n]
}
myColors <- gg_color_hue(10)

cutoff_snvs <- cutoff_svs <- subclonal_cutoff <- 0.7
cutoff_cna <-  0.943

basedir <- '001_svclone_results/'
wd_truth <- '001_truth/'
resultFolder <- '001_svclone_results/001_mixes/'
clonal_basedir <- snvFolder <- '001_svclone_results/001_mixes_clonal/'
bb_out <- 'battenberg_out/'
base_pyc_results <- '001_pyclone_results/'
pycInputFolder <-  '001_pyclone_results/001_q20_mq20_pyclone/'
```


```{r load_svclone_mix_results}
sv_data <- NULL
snv_data <- NULL
for(i in seq(0.1, 0.9, 0.1)) {
    mix <- paste(i, 1-i, sep='-')
    basename <- paste('001bM_p', i, '_001gM_p', 1-i, sep='')

    ################################################################################################
    # SVs
    ################################################################################################
    ccubeResultRDataFile = paste0(resultFolder, '/', basename, '/ccube_out/ccube_sv_results.RData')
    load(ccubeResultRDataFile)
    svs <- data.frame(data='svs', doubleBreakPtsRes$ssm, mix=mix, p1=i, p2=1-i)

    svin <- read.delim(paste0(resultFolder, '/', basename, '/', basename, '_filtered_svs.tsv'))
    svin$sv <- paste(svin$chr1, svin$pos1, svin$dir1, svin$chr2, svin$pos2, svin$dir2, sep=':')
    svin$mutation_id <- paste(paste(svin$chr1, svin$pos1, svin$dir1, sep=':'),
                              paste(svin$chr2, svin$pos2, svin$dir2, sep=':'), sep='_')
    svs <- inner_join(svs, svin, by='mutation_id')
    svs <- attach_001_ground_truth_svs(wd_truth, svs)

    pp <- read.delim(paste0(resultFolder, basename, '/purity_ploidy.txt'))
    pur <- pp$purity; pl <- pp$ploidy
    af <- 1 - (pur / pl)
    gains <- svs$classification%in%c('DUP','INTDUP')

    svs$norm1_adjusted <- svs$norm1; svs$norm2_adjusted <- svs$norm2;
    svs$norm1_adjusted[gains] <- svs$norm1_adjusted[gains] * af
    svs$norm2_adjusted[gains] <- svs$norm2_adjusted[gains] * af
    svs$norm_cn <- 2; svs$norm_cn[svs$chr1%in%c('X','Y','MT')] <- 1
    svs$true_ccf <- 1
    svs[svs$sample%in%'bM_only', 'true_ccf'] <- svs$p1[svs$sample%in%'bM_only']
    svs[svs$sample%in%'gM_only', 'true_ccf'] <- svs$p2[svs$sample%in%'gM_only']
    rownames(svs) <- 1:nrow(svs)

    x <- t(data.frame(apply(svs, 1, get_true_sc_cn_by_side, pur=pur, side=1)))
    svs$true_cn1 <- x[,1]
    svs$wtotal_cn1 <- x[,2]
    svs$pv1 <- x[,3]
    svs$best_ccf1 <- x[,4]

    x <- t(data.frame(apply(svs, 1, get_true_sc_cn_by_side, pur=pur, side=2)))
    svs$true_cn2 <- x[,1]
    svs$wtotal_cn2 <- x[,2]
    svs$pv2 <- x[,3]
    svs$best_ccf2 <- x[,4]

    sv_data <- rbind(sv_data, svs)

    ################################################################################################
    # SNVs
    ################################################################################################
    ccubeResultRDataFile = paste0(clonal_basedir, basename, '/ccube_out/snvs/', basename, '_ccube_snv_results.RData')
    load(ccubeResultRDataFile)

    snvs <- data.frame(data='snvs', snvRes$ssm, sv=snvRes$ssm$mutation_id, mix=mix, p1=i, p2=1-i)
    snvs <- attach_001_ground_truth_snvs(wd_truth, snvs)

    snvin <- read.delim(paste0(resultFolder, '/', basename, '/', basename, '_filtered_snvs.tsv'))
    snvin$mutation_id <- paste(snvin$chrom, snvin$pos, sep='_')
    snvs <- inner_join(snvs, snvin[,c('mutation_id', 'gtype')], by='mutation_id')

    snvs$true_ccf <- 1
    snvs[snvs$sample%in%'bM_only','true_ccf'] <- i
    snvs[snvs$sample%in%'gM_only','true_ccf'] <- 1-i

    snvs$adjusted_support <- snvs$var_counts
    snvs$adjusted_depth <- snvs$var_counts + snvs$ref_counts

    snvs$norm_cn <- snvs$normal_cn
    x <- t(data.frame(apply(snvs, 1, get_true_sc_cn)))
    snvs$true_cn <- x[,1]; snvs$wtotal_cn <- x[,2]
    snvs$pv <- x[,3]; snvs$best_ccf <- x[,4]

    snv_data <- rbind(snv_data, snvs)
}
snv_clonal <- snv_data
```


```{r load_pyclone_results}
pyc_all <- NULL
for(p1 in seq(0.1, 0.9, 0.1)) {
    pycloneFolder <- paste0(base_pyc_results, '001bM_p', p1)
    p2 <- 1 - p1; mix <- paste(p1, p2, sep='-')

    pp <- read.delim(paste0(resultFolder, '/001bM_p', p1, '_001gM_p', p2, '/purity_ploidy.txt'))
    pyc <- read.delim(file = paste0(pycloneFolder, '/pyclone_results.tsv'), sep = "\t")
    pyc$purity <- pp$purity

    sampleTsvFile <- paste0(pycInputFolder, "001bM_p", p1, "_pyclone_input.txt")
    sampleTsv <- read.delim(sampleTsvFile, stringsAsFactors = F)
    sampleTsv <- mutate(sampleTsv, vaf = var_counts/(ref_counts+var_counts))
    allData <- left_join(pyc, sampleTsv, by = "mutation_id")
    rm(sampleTsv)
    allData <- dplyr::filter(allData, !is.na(average_ccf))

    allData <- dplyr::mutate(allData,
                    multiplicity = GetMultFromCcf(bn = var_counts,
                                                  dn = var_counts + ref_counts,
                                                  ccf = ccf,
                                                  major_cn = major_cn,
                                                  minor_cn = minor_cn,
                                                  purity = purity))
    pyc <- allData;

    tmp1 <- as.data.frame(table(pyc$cluster_id), stringsAsFactors = F)
    tmp2 <- as.data.frame(table(pyc$average_ccf), stringsAsFactors = F)
    tmp <- left_join(tmp1, tmp2, by ="Freq")
    colnames(tmp) <- c('cluster_id', 'n_ssms', 'cluster_ccf')
    pyc <- merge(pyc, tmp, by='cluster_id')

    pyc$sv <- pyc$mutation_id; pyc$true_ccf <- 1
    pyc <- attach_001_ground_truth_snvs(wd_truth, pyc)
    pyc[pyc$sample%in%'bM_only', 'true_ccf'] <- p1
    pyc[pyc$sample%in%'gM_only', 'true_ccf'] <- p2

    pyc$norm_cn <- pyc$normal_cn
    pyc$gtype <- paste(pyc$major_cn, pyc$minor_cn, 1, sep=',')
    pyc$adjusted_support <- pyc$var_counts
    pyc$adjusted_depth <- pyc$var_counts + pyc$ref_counts

    x <- t(data.frame(apply(pyc, 1, get_true_sc_cn)))
    pyc$true_cn <- x[,1]
    pyc$wtotal_cn <- x[,2]
    pyc$pv <- x[,3]
    pyc$best_ccf <- x[,4]

    pyc$best_guess_subclonal <- pyc$best_ccf < subclonal_cutoff
    pyc$pyc_guess_subclonal <- pyc$ccf < subclonal_cutoff

    pyc$mix <- mix
    pyc_all <- rbind(pyc_all, pyc)
}

for(nclus in c('4clus', '5clus')) {
    pycloneFolder <- paste0(base_pyc_results, nclus)
    mix <- nclus

    pp <- read.delim(paste0(basedir, nclus, '_out/purity_ploidy.txt'))
    pyc_file <- paste0(pycloneFolder, '/pyclone_results.tsv')
    if(!file.exists(pyc_file)) {
        proc_pyclone_results(pycloneFolder, pp$purity)
    }
    pyc <- read.delim(file = pyc_file, sep = "\t")
    pyc$purity <- pp$purity

    sampleTsvFile <- paste0(pycInputFolder, nclus, '_pyclone_input.txt')
    sampleTsv <- read.delim(sampleTsvFile, stringsAsFactors = F)
    sampleTsv <- mutate(sampleTsv, vaf = var_counts/(ref_counts+var_counts))
    allData <- left_join(pyc, sampleTsv, by = "mutation_id")
    allData <- dplyr::filter(allData, !is.na(average_ccf))

    allData <- dplyr::mutate(allData,
                    multiplicity = GetMultFromCcf(bn = var_counts,
                                                  dn = var_counts + ref_counts,
                                                  ccf = ccf,
                                                  major_cn = major_cn,
                                                  minor_cn = minor_cn,
                                                  purity = purity))
    pyc <- allData
    pyc$chr <- as.numeric(sapply(pyc$mutation_id, function(x){strsplit(x, '_')[[1]][1]}))

    tmp1 <- as.data.frame(table(pyc$cluster_id), stringsAsFactors = F)
    tmp2 <- as.data.frame(table(pyc$average_ccf), stringsAsFactors = F)
    tmp <- left_join(tmp1, tmp2, by ="Freq")
    colnames(tmp) <- c('cluster_id', 'n_ssms', 'cluster_ccf')
    pyc <- merge(pyc, tmp, by='cluster_id')

    pyc$sv <- pyc$mutation_id; pyc$true_ccf <- 1
    pyc <- attach_001_ground_truth_snvs(wd_truth, pyc)

    if(nclus == '4clus') {
        pyc[pyc$sample%in%'bM_only' & pyc$chr %% 2 == 1, 'true_ccf'] <- 0.2
        pyc[pyc$sample%in%'bM_only' & pyc$chr %% 2 == 0, 'true_ccf'] <- 0.6
        pyc[pyc$sample%in%'gM_only' & pyc$chr %% 2 == 1, 'true_ccf'] <- 0.4
    } else {
        pyc[pyc$sample%in%'bM_only' & pyc$chr %% 2 == 1, 'true_ccf'] <- 0.8
        pyc[pyc$sample%in%'bM_only' & pyc$chr %% 2 == 0, 'true_ccf'] <- 0.6
        pyc[pyc$sample%in%'gM_only' & pyc$chr %% 2 == 1, 'true_ccf'] <- 0.2
        pyc[pyc$sample%in%'gM_only' & pyc$chr %% 2 == 0, 'true_ccf'] <- 0.4
    }

    pyc$norm_cn <- pyc$normal_cn
    pyc$gtype <- paste(pyc$major_cn, pyc$minor_cn, 1, sep=',')
    pyc$adjusted_support <- pyc$var_counts
    pyc$adjusted_depth <- pyc$var_counts + pyc$ref_counts

    x <- t(data.frame(apply(pyc, 1, get_true_sc_cn)))
    pyc$true_cn <- x[,1]
    pyc$wtotal_cn <- x[,2]
    pyc$pv <- x[,3]
    pyc$best_ccf <- x[,4]

    pyc$best_guess_subclonal <- pyc$best_ccf < subclonal_cutoff
    pyc$pyc_guess_subclonal <- pyc$ccf < subclonal_cutoff
    pyc$chr <- NULL

    pyc$mix <- mix
    pyc_all <- rbind(pyc_all, pyc)
}
```


```{r load_cna_results}
bcn <- read.delim('001_truth/001bM_subclones.txt')[,2:14]
gcn <- read.delim('001_truth/001gM_subclones.txt')[,2:14]

# filter out CN neutral
bcn <- bcn[!((bcn$frac1 == 1 & bcn$nMaj1_A == 1 & bcn$nMin1_A == 1) |
           (bcn$chr %in% c('X','Y') &  bcn$frac1 == 1 & bcn$nMaj1_A == 1)), ]
gcn <- gcn[!((gcn$frac1 == 1 & gcn$nMaj1_A == 1 & gcn$nMin1_A == 1) |
           (gcn$chr %in% c('X','Y') &  gcn$frac1 == 1 & gcn$nMaj1_A == 1)), ]

boundary <- 5000
bx <- GRanges(seqnames = bcn$chr, ranges=IRanges(start = bcn$startpos,
                                                 end = bcn$endpos))
gx <- GRanges(seqnames = gcn$chr, ranges=IRanges(start = gcn$startpos,
                                                 end = gcn$endpos))

olaps_both <- findOverlaps(bx, gx, maxgap=boundary, type='equal')
olaps_starts <- findOverlaps(bx, gx, maxgap=boundary, type='start')
olaps_ends <- findOverlaps(bx, gx, maxgap=boundary, type='end')

hits <- data.frame(olaps_both)
match <- ((bcn[hits$queryHits,]$nMaj1_A == gcn[hits$subjectHits,]$nMaj1_A) &
         (bcn[hits$queryHits,]$nMin1_A == gcn[hits$subjectHits,]$nMin1_A)) &
        ((bcn[hits$queryHits,]$nMaj2_A == gcn[hits$subjectHits,]$nMaj2_A) &
         (bcn[hits$queryHits,]$nMin2_A == gcn[hits$subjectHits,]$nMin2_A)) |
        (is.na(bcn[hits$queryHits,]$nMaj2_A) == is.na(gcn[hits$subjectHits,]$nMaj2_A))
match[is.na(match)] <- FALSE

# shared CNs lookup
p1cn <- bcn[queryHits(olaps_both),]
p1cn$sample <- 'shared'

part1 <- union(queryHits(olaps_starts), queryHits(olaps_ends))
part2 <- union(subjectHits(olaps_starts), subjectHits(olaps_ends))

p1cn <- rbind(p1cn, data.frame(bcn[-c(part1),], sample='bM_only'))
p1cn <- rbind(p1cn, data.frame(gcn[-c(part2),], sample='gM_only'))
rownames(p1cn) <- 1:nrow(p1cn)

# get the dominant copynumber
p1cn$major <- p1cn$nMaj1_A
p1cn$minor <- p1cn$nMin1_A

dom_cn <- t(apply(p1cn[grep('only', p1cn$sample),], 1,
                function(x){
                    if(x['frac1_A'] > 0.5){
                        return(c(x['nMaj1_A'], x['nMin1_A']))
                    } else {
                        return(c(x['nMaj2_A'], x['nMin2_A']))
                    }}))
dom_cn <- apply(dom_cn, 2, as.numeric)

p1cn[grep('only', p1cn$sample),'major'] <- dom_cn[,1]
p1cn[grep('only', p1cn$sample),'minor'] <- dom_cn[,2]

p1x <- GRanges(seqnames = p1cn$chr, ranges = IRanges(start = p1cn$startpos, end = p1cn$endpos))

# now compare with 10-90 sample
all_comps <- NULL
for(i in seq(0.1, 0.9, 0.1)) {
    basename <- paste('001bM_p', i, '_001gM_p', 1-i, sep='')
    cn <- read.delim(paste0(bb_out, basename, '_subclones.txt'))[,2:14]
    cn <- cn[!((cn$frac1 == 1 & cn$nMaj1_A == 1 & cn$nMin1_A == 1) |
               (cn$chr %in% c('X','Y') &  cn$frac1 == 1 & cn$nMaj1_A == 1)), ]
    cx <- GRanges(seqnames = cn$chr, ranges = IRanges(start = cn$startpos, end = cn$endpos))
    olaps <- findOverlaps(cx, p1x, maxgap=boundary, type='equal')
    tmp <- cn[queryHits(olaps),]
    tmp$sample <- p1cn[subjectHits(olaps),]$sample
    tmp$subjectHits <- subjectHits(olaps)

    tmp1 <- tmp[grep('only',tmp$sample),]
    fracs <- apply(tmp1, 1,
                function(x){
                    maj <- p1cn[x['subjectHits'],'major']
                    min <- p1cn[x['subjectHits'],'minor']
                    cmaj1 <- as.numeric(x['nMaj1_A'])
                    cmin1 <- as.numeric(x['nMin1_A'])
                    cmaj2 <- as.numeric(x['nMaj2_A'])
                    cmin2 <- as.numeric(x['nMin2_A'])
                    if(maj == cmaj1 & min == cmin1){
                        return(x['frac1_A'])
                    } else if(!is.na(cmaj2) & (maj == cmaj2 & min == cmin2)) {
                        return(x['frac2_A'])
                    } else {return(NA)}
                })
    tmp$true_ccf <- 1
    tmp[tmp$sample%in%'bM_only', 'true_ccf'] <- i
    tmp[tmp$sample%in%'gM_only', 'true_ccf'] <- 1-i
    tmp$ccf <- tmp$frac1_A
    tmp[grep('only',tmp$sample), 'ccf'] <- fracs
    all_comps <- rbind(all_comps, data.frame(tmp, mix=paste0(i, '-', 1-i)))
}

res <- all_comps[!is.na(all_comps$ccf),]
res$ccf <- as.numeric(res$ccf)
res <- data.table(res)[,mean(true_ccf - ccf), by=c('mix')]

p1 <- ggplot(res, aes(mix, V1, group=1)) +
        geom_point() + geom_line() + ylim(-0.5, 0.5) + theme_bw() + ggtitle('mean_ccf_error')

res <- all_comps[!is.na(all_comps$ccf),]
res$is_subclonal_truth <- res$true_ccf<subclonal_cutoff
res$is_subclonal_bb <- res$ccf<subclonal_cutoff

cna_data <- all_comps[!is.na(all_comps$ccf),]
```

## Figure 2b

Optimal versus ground truth CCF distributions. From left to right, by column: i) SV theoretical mixture truth CCFs, ii) 'optimal' SV CCFs based on sample membership, iii) SNV theoretical mixture truth CCF, iv) 'optimal' SNV CCFs based on sample membership.

This figure is an extended version of the figure that appears in the paper.

```{r Figure2b, fig.height=20, fig.width=10}
cols <- brewer.pal(10, 'RdBu')
names(cols) <- seq(0.1, 1, 0.1)

ccf_met <- NULL
opts <- NULL
for(p1 in seq(0.1, 0.9, 0.1)) {
    p2 <- 1 - p1
    mix <- paste(p1, p2, sep='-')
    x <- sv_data[sv_data$mix%in%mix,]
    x$best_ccf <- as.numeric(apply(x[,c('best_ccf1', 'best_ccf2')], 1, mean))
    opts[[paste0(mix, '_true_sv')]] <-
            ggplot(x, aes(true_ccf, fill=as.factor(true_ccf))) +
                geom_histogram(colour='grey', alpha=0.6, position='identity', binwidth=0.05) +
                ylim(0,80) + geom_vline(xintercept=c(1, p1, p2), colour='blue', size=1, alpha=0.2) +
                scale_x_continuous(breaks = seq(0, 2, 0.5), lim=c(0,2)) + theme_bw() +
                theme(legend.position = 'none', axis.text = element_text(size = 12)) +
                scale_fill_manual(values = cols) +
                ylab('') + xlab('')
    opts[[paste0(mix, '_best_sv')]] <-
        ggplot(x, aes(best_ccf, fill=as.factor(true_ccf))) +
                    geom_histogram(colour='grey', alpha=0.6, position='identity', binwidth=0.05) +
                    ylim(0,20) + geom_vline(xintercept=c(1, p1, p2), colour='blue', size=1, alpha=0.2) +
                    scale_x_continuous(breaks = seq(0, 2, 0.5), lim=c(0,2)) + theme_bw() +
                    theme(legend.position = 'none', axis.text = element_text(size = 12)) +
                scale_fill_manual(values = cols) +
                ylab('') + xlab('')

    ccf_err <- mean(x$true_ccf - x$best_ccf)
    mult_err <- mean(c(x$true_cn1 - x$ccube_mult1, x$true_cn2 - x$ccube_mult2))
    ccf_met <- rbind(ccf_met, data.frame(data='SVs', mix, ccf_err))

    x <- snv_data[snv_data$mix%in%mix,]
    opts[[paste0(mix, '_true_snv')]] <-
            ggplot(x, aes(true_ccf, fill=as.factor(true_ccf))) +
                geom_histogram(colour='grey', alpha=0.6, position='identity', binwidth=0.05) +
                geom_vline(xintercept=c(1, p1, p2), colour='blue', size=1, alpha=0.2) +
                scale_x_continuous(breaks = seq(0, 2, 0.5), lim=c(0,2)) + theme_bw() + ylim(0,4000) +
                theme(legend.position = 'none', axis.text = element_text(size = 12)) +
                scale_fill_manual(values = cols) +
                ylab('') + xlab('')
    opts[[paste0(mix, '_best_snv')]] <-
        ggplot(x, aes(best_ccf, fill=as.factor(true_ccf))) +
                    geom_histogram(colour='grey', alpha=0.6, position='identity', binwidth=0.05) +
                    geom_vline(xintercept=c(1, p1, p2), colour='blue', size=1, alpha=0.2) +
                    scale_x_continuous(breaks = seq(0, 2, 0.5), lim=c(0,2)) + theme_bw() + ylim(0,500) +
                    theme(legend.position = 'none', axis.text = element_text(size = 12)) +
                scale_fill_manual(values = cols) +
                ylab('') + xlab('')

    ccf_err <- mean(x$true_ccf - x$best_ccf)
    ccf_met <- rbind(ccf_met, data.frame(data='SNVs', mix, ccf_err))
}

# add hiclus mixes
for(p1 in c('4clus', '5clus')) {
    x <- load_hiclus_svs(p1, basedir)
    x$best_ccf <- as.numeric(apply(x[,c('best_ccf1', 'best_ccf2')], 1, mean))
    clones <- as.numeric(names(table(x$true_ccf)))
    opts[[paste0(p1, '_true_sv')]] <-
            ggplot(x, aes(true_ccf, fill=as.factor(true_ccf))) +
                    geom_histogram(colour='grey', alpha=0.6, position='identity', binwidth=0.05) +
                    ylim(0,80) + geom_vline(xintercept=clones, colour='blue', size=1, alpha=0.2) +
                    scale_x_continuous(breaks = seq(0, 2, 0.5), lim=c(0,2)) + theme_bw() +
                    theme(legend.position = 'none', axis.text = element_text(size = 12)) +
                scale_fill_manual(values = cols) +
                ylab('') + xlab('')
    opts[[paste0(p1, '_best_sv')]] <-
            ggplot(x, aes(best_ccf, fill=as.factor(true_ccf))) +
                    geom_histogram(colour='grey', alpha=0.6, position='identity', binwidth=0.05) +
                    ylim(0,20) + geom_vline(xintercept=clones, colour='blue', size=1, alpha=0.2) +
                    scale_x_continuous(breaks = seq(0, 2, 0.5), lim=c(0,2)) + theme_bw() +
                    theme(legend.position = 'none', axis.text = element_text(size = 12)) +
                scale_fill_manual(values = cols) +
                ylab('') + xlab('')

    ccf_err <- mean(x$true_ccf - x$best_ccf)
    mult_err <- mean(c(x$true_cn1 - x$ccube_mult1, x$true_cn2 - x$ccube_mult2))
    ccf_met <- rbind(ccf_met, data.frame(data='SVs', mix=p1, ccf_err))

    x <- load_hiclus_snvs(p1, basedir)
    opts[[paste0(p1, '_true_snv')]] <-
        ggplot(x, aes(true_ccf, fill=as.factor(true_ccf))) +
            geom_histogram(colour='grey', alpha=0.6, position='identity', binwidth=0.05) +
            geom_vline(xintercept=clones, colour='blue', size=1, alpha=0.2) +
            scale_x_continuous(breaks = seq(0, 2, 0.5), lim=c(0,2)) + theme_bw() + ylim(0,4000) +
            theme(legend.position = 'none', axis.text = element_text(size = 12)) +
            scale_fill_manual(values = cols) +  ylab('') + xlab('')
    opts[[paste0(p1, '_best_snv')]] <-
        ggplot(x, aes(best_ccf, fill=as.factor(true_ccf))) +
                    geom_histogram(colour='grey', alpha=0.6, position='identity', binwidth=0.05) +
                    geom_vline(xintercept=clones, colour='blue', size=1, alpha=0.2) +
                    scale_x_continuous(breaks = seq(0, 2, 0.5), lim=c(0,2)) + theme_bw() + ylim(0,500) +
                    theme(legend.position = 'none', axis.text = element_text(size = 12)) +
                    scale_fill_manual(values = cols) +  ylab('') + xlab('')

    ccf_err <- mean(x$true_ccf - x$best_ccf)
    ccf_met <- rbind(ccf_met, data.frame(data='SNVs', mix=p1, ccf_err))
}
do.call("grid.arrange", c(opts, ncol=4))
```

CCF mean error

```{r}
print(data.table(ccf_met)[,mean(ccf_err),by=c('data')])
```


## Figure 2c

True minus optimal CCF averages for SNVs (red) and SVs (blue).

```{r Figure2c, fig.height=5, fig.width=6}
cm <- melt(ccf_met)
mixes <- sapply(cm$mix[grep('-',cm$mix)],
                function(x){mix <- strsplit(x, '-')[[1]];
                   paste(as.numeric(mix[1])*100, as.numeric(mix[2])*100, sep='-') })
cm$mix[grep('-',cm$mix)] <- as.character(mixes)
cm$mix <- factor(cm$mix, levels=c('10-90', '20-80', '30-70', '40-60', '50-50',
                                  '60-40', '70-30', '80-20', '90-10', '4clus', '5clus'))

ggplot(cm, aes(mix, value, group=data, colour=data)) +
        geom_line() + geom_point() + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('optimal - estimated CCF') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=10),
              axis.text.y = element_text(size=14), legend.title = element_blank(),
              legend.text = element_text(size=14)) +
        scale_colour_brewer(palette = 'Dark2')
```

## Figure 2d

ROC curves showing sensitivity and specificity of classifying variants as subclonal:

- red: SNVs
- gold: SV max
- green: SV mean
- purple: SV min

```{r Figure2d, fig.height=5, fig.width=7}
source('https://gist.githubusercontent.com/charly06/91578196fc615c5a79c7174318be4349/raw/d96a98c2933f5af141eac91af83c7895062d68a5/ggrocs.R')

x <- sv_data
x$best_ccf_mean <- as.numeric(apply(x[,c('best_ccf1', 'best_ccf2')], 1, mean))
x$best_ccf_min <- as.numeric(apply(x[,c('best_ccf1', 'best_ccf2')], 1, min))
x$best_ccf_max <- as.numeric(apply(x[,c('best_ccf1', 'best_ccf2')], 1, max))
x$is_subclonal_truth <- !x$sample%in%'shared'

y <- snv_data
y$is_subclonal_truth <- !y$sample%in%'shared'

rocs <- NULL
rocs[['SV_ccf_mean']] <- roc(x$is_subclonal_truth, x$best_ccf_mean)
rocs[['SV_ccf_min']] <- roc(x$is_subclonal_truth, x$best_ccf_min)
rocs[['SV_ccf_max']] <- roc(x$is_subclonal_truth, x$best_ccf_max)
rocs[['SNV_ccf']] <- roc(y$is_subclonal_truth, y$best_ccf)

tsize = 15
ggrocs(rocs) + theme(axis.text = element_text(size = tsize),
                     legend.text = element_text(size = tsize)) +
               geom_step(size=1) +
               scale_colour_brewer(palette = 'RdBu')
```

Respective AUCs

```{r}
sapply(rocs, auc)
```

SV CCF mean optimal cutoff

```{r}
pROC::coords(rocs[['SV_ccf_mean']], "b", ret="t", best.method="youden")
```

SV CCF min optimal cutoff

```{r}
pROC::coords(rocs[['SV_ccf_min']], "b", ret="t", best.method="youden")
```

SV CCF min optimal cutoff

```{r}
pROC::coords(rocs[['SV_ccf_max']], "b", ret="t", best.method="youden")
```

SNV CCF optimal cutoff

```{r}
pROC::coords(rocs[['SNV_ccf']], "b", ret="t", best.method="youden")
```

## Supplementary Table 3

Percent of SVs with heterogenous copy-number states across both SV sides.

```{r SuppTable3}
# add 001bM and gM data
bm <- read.delim('001_svclone_results/001bM/001bM_filtered_svs.tsv', sep='\t')
gm <- read.delim('001_svclone_results/001gM/001gM_filtered_svs.tsv', sep='\t')

hetCN <- data.table(mix='001bM', heterogenousCN=bm$gtype1 != bm$gtype2)
hetCN <- rbind(hetCN,
               data.table(mix='001gM', heterogenousCN=gm$gtype1 != gm$gtype2))
hetCN <- rbind(hetCN,
               data.table(mix=sv_data$mix, heterogenousCN=sv_data$gtype1 != sv_data$gtype2))
hetCN <- hetCN[, sum(heterogenousCN)/length(heterogenousCN), by=c('mix')]
hetCN$V1 <- round(hetCN$V1, 2)
colnames(hetCN)[2] <- 'heterogenous_CN_percent'

print(hetCN)
```



```{r calc_se_metrics}
source('R/CN_perturbation_experiments.R')

se_basedir <- '001_svclone_results/001_mixes_single_end/'
# uncomment below to generate data
# SVclone/cluster_with_ccube.R will then
# need to be run on output
# generate_se_svs(basedir, se_basedir)

svse_data <- NULL
for(i in seq(0.1, 0.9, 0.1)) {
    mix <- paste(i, 1-i, sep='-')
    basename <- paste('001bM_p', i, '_001gM_p', 1-i, sep='')
    for(side in c('side1', 'side2')) {
        ccubeResultRDataFile = paste0(se_basedir, basename,
                                      paste('/ccube', side, 'out/', sep='_'),
                                      'ccube_out/snvs/', basename, '_ccube_snv_results.RData')
        load(ccubeResultRDataFile)

        svs <- data.frame(data='snvs', side=side, snvRes$ssm, 
                          sv=snvRes$ssm$mutation_id, mix=mix, p1=i, p2=1-i)

        svin <- read.delim(paste0(resultFolder, basename, '/', basename, '_filtered_svs.tsv'))
        svin$mutation_id <- paste(svin$chr1, svin$pos1, svin$dir1, svin$chr2, svin$pos2, svin$dir2, sep=':')
        svs <- inner_join(svs, svin, by='mutation_id')
        svs$sv <- svs$mutation_id; svs$side <- side
        svs <- attach_001_ground_truth_svs(wd_truth, svs)

        svs$true_ccf <- 1
        svs[svs$sample%in%'bM_only','true_ccf'] <- i
        svs[svs$sample%in%'gM_only','true_ccf'] <- 1-i

        svs$adjusted_support <- svs$var_counts
        svs$adjusted_depth <- svs$var_counts + svs$ref_counts

        svs$norm_cn <- svs$normal_cn
        svs$gtype <- ifelse(side == 'side1', svs$gtype1, svs$gtype2)
        x <- t(data.frame(apply(svs, 1, get_true_sc_cn)))
        svs$true_cn <- x[,1]; svs$wtotal_cn <- x[,2]
        svs$pv <- x[,3]; svs$best_ccf <- x[,4]

        svse_data <- rbind(svse_data, svs)
    }
}

metrics <- NULL
for(p1 in seq(0.1, 0.9, 0.1)) {
    p2 <- 1 - p1
    mix <- paste(p1, p2, sep='-')
    mixname <- paste(p1*100, p2*100, sep='-')

    x <- sv_data[sv_data$mix%in%mix, ]
    metrics <- rbind(metrics, calc_metrics_3clus(x, mix, 'SVs_dual_end'))

    x <- svse_data[svse_data$mix%in%mix & svse_data$side%in%'side1', ]
    metrics <- rbind(metrics, calc_metrics_3clus(x, mix, 'SVs_single_end1', type='snv'))

    x <- svse_data[svse_data$mix%in%mix & svse_data$side%in%'side2', ]
    metrics <- rbind(metrics, calc_metrics_3clus(x, mix, 'SVs_single_end2', type='snv'))
}
```

## Figure 4

Performance comparison of the dual-end SV model versus single-end.

```{r Figure4, fig.height=8, fig.width=12}
mm <- melt(metrics)
mm$variable <- as.character(mm$variable)
mm$mix <- sapply(mm$mix, function(x){tmp <- as.numeric(strsplit(x, '-')[[1]]); paste(tmp[1]*100, tmp[2]*100, sep='-')})
mm[mm$variable%like%'mean_mult_error',]$variable <- 'mean multiplicity error'
mm[mm$variable%like%'mean_ccf',]$variable <- 'mean variant CCF error'
mm[mm$variable%like%'sens',]$variable <- 'is subclonal sensitivity'
mm[mm$variable%like%'spec',]$variable <- 'is subclonal specificity'
mm[mm$variable%like%'clus_ccf',]$variable <- 'mean cluster CCF error'
mm[mm$variable%like%'clus_num',]$variable <- 'cluster number error'

colpal <- 'Dark2'
g1 <- ggplot(mm[mm$variable%like%'cluster num',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-3, 3) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated cluster number') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g2 <- ggplot(mm[mm$variable%like%'cluster CCF',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-0.6, 0.6) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated cluster CCF') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g3 <- ggplot(mm[mm$variable%like%'mean var',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-0.6, 0.6) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated CCF') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g4 <- ggplot(mm[mm$variable%like%'mult',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-0.6, 0.6) +
        geom_hline(yintercept=0, colour='grey') + ylab('optimal - estimated multiplicity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g5 <- ggplot(mm[mm$variable%like%'sens',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(0,1) +
        ylab('sensitivity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g6 <- ggplot(mm[mm$variable%like%'spec',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(0,1) +
        ylab('specificity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=7)) +
        scale_colour_brewer(palette = colpal)

grid.arrange(g1, g3, g5,
             g2, g4, g6, ncol=3)
```

## Supplementary Figure 10

Comparing ROC curves of dual-end and single-end models.

```{r SuppFigure10, fig.height=7, fig.width=6}
x <- sv_data
x$ccf_mean <- apply(x[,c('ccube_ccf1', 'ccube_ccf2')], 1, mean)
x$is_subclonal_truth <- !x$sample%in%'shared'

y <- svse_data
y$is_subclonal_truth <- !y$sample%in%'shared'

rocs <- NULL

rocs[['SVs_dual_end']] <- roc(x$is_subclonal_truth, x$ccf_mean)
rocs[['SVs_single_end']] <- roc(y$is_subclonal_truth, y$ccube_ccf)

tsize = 15
ggrocs(rocs) + theme(axis.text = element_text(size = tsize),
                     legend.text = element_text(size = tsize),
                     legend.title = element_blank(),
                     legend.position = 'bottom') + geom_step(size=1) +
                scale_colour_brewer(palette = colpal)
```
Respective AUCs

```{r de_vs_se_AUC}
sapply(rocs, auc)
```


```{r calc_metrics_clonal_vs_subcl}
svc_data <- NULL
for(i in seq(0.1, 0.9, 0.1)) {
    mix <- paste(i, 1-i, sep='-')
    basename <- paste('001bM_p', i, '_001gM_p', 1-i, sep='')

    ccubeResultRDataFile = paste0(clonal_basedir, basename, '/ccube_out/',
                                  basename, '_ccube_sv_results.RData')
    load(ccubeResultRDataFile)
    svs <- data.frame(data='svs', doubleBreakPtsRes$ssm, mix=mix, p1=i, p2=1-i)

    svin <- read.delim(paste0(clonal_basedir, basename, '/', basename, '_filtered_svs.tsv'))
    svin$sv <- paste(svin$chr1, svin$pos1, svin$dir1, svin$chr2, svin$pos2, svin$dir2, sep=':')
    svin$mutation_id <- paste(paste(svin$chr1, svin$pos1, svin$dir1, sep=':'),
                              paste(svin$chr2, svin$pos2, svin$dir2, sep=':'), sep='_')
    svs <- inner_join(svs, svin, by='mutation_id')
    svs <- attach_001_ground_truth_svs(wd_truth, svs)

    pp <- read.delim(paste0(resultFolder, basename, '/purity_ploidy.txt'))
    pur <- pp$purity; pl <- pp$ploidy
    af <- 1 - (pur / pl)
    gains <- svs$classification%in%c('DUP','INTDUP')

    svs$norm1_adjusted <- svs$norm1; svs$norm2_adjusted <- svs$norm2;
    svs$norm1_adjusted[gains] <- svs$norm1_adjusted[gains] * af
    svs$norm2_adjusted[gains] <- svs$norm2_adjusted[gains] * af
    svs$vaf1 <- svs$adjusted_support / (svs$norm1_adjusted + svs$adjusted_support)
    svs$vaf2 <- svs$adjusted_support / (svs$norm2_adjusted + svs$adjusted_support)
    svs$norm_cn <- 2; svs$norm_cn[svs$chr1%in%c('X','Y','MT')] <- 1
    svs$true_ccf <- 1
    svs[svs$sample%in%'bM_only', 'true_ccf'] <- svs$p1[svs$sample%in%'bM_only']
    svs[svs$sample%in%'gM_only', 'true_ccf'] <- svs$p2[svs$sample%in%'gM_only']
    svs$ccf_true <- svs$true_ccf
    rownames(svs) <- 1:nrow(svs)

    x <- t(data.frame(apply(svs, 1, get_true_sc_cn_by_side, pur=pur, side=1)))
    svs$true_cn1 <- x[,1]
    svs$wtotal_cn1 <- x[,2]
    svs$pv1 <- x[,3]
    svs$best_ccf1 <- x[,4]

    x <- t(data.frame(apply(svs, 1, get_true_sc_cn_by_side, pur=pur, side=2)))
    svs$true_cn2 <- x[,1]
    svs$wtotal_cn2 <- x[,2]
    svs$pv2 <- x[,3]
    svs$best_ccf2 <- x[,4]

    svc_data <- rbind(svc_data, svs)
}

metrics <- NULL
for(p1 in seq(0.1, 0.9, 0.1)) {
    p2 <- 1 - p1
    minor <- min(p1, p2)
    major <- max(p1, p2)
    mix <- paste(p1, p2, sep='-')
    mixname <- paste(p1*100, p2*100, sep='-')

    x <- sv_data[sv_data$mix%in%mix, ]
    metrics <- rbind(metrics, calc_metrics_3clus(x, mix, 'SVs_all_CN'))

    x <- svc_data[svc_data$mix%in%mix, ]
    metrics <- rbind(metrics, calc_metrics_3clus(x, mix, 'SVs_clonal_CN'))
}
```

## Figure 5

Performance for SVs that have clonal background copy-numbers vs. all SVs (clonal + subclonal background copy-numbers).

```{r Figure5, fig.height=8, fig.width=12}
mm <- melt(metrics)
mm$variable <- as.character(mm$variable)
mm$mix <- sapply(mm$mix, function(x){tmp <- as.numeric(strsplit(x, '-')[[1]]); paste(tmp[1]*100, tmp[2]*100, sep='-')})
mm[mm$variable%like%'mean_mult_error',]$variable <- 'mean multiplicity error'
mm[mm$variable%like%'mean_ccf',]$variable <- 'mean variant CCF error'
mm[mm$variable%like%'sens',]$variable <- 'is subclonal sensitivity'
mm[mm$variable%like%'spec',]$variable <- 'is subclonal specificity'
mm[mm$variable%like%'clus_ccf',]$variable <- 'mean cluster CCF error'
mm[mm$variable%like%'clus_num',]$variable <- 'cluster number error'

g1 <- ggplot(mm[mm$variable%like%'cluster num',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-3, 3) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated cluster number') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g2 <- ggplot(mm[mm$variable%like%'cluster CCF',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated cluster CCF') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g3 <- ggplot(mm[mm$variable%like%'mean var',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated CCF') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g4 <- ggplot(mm[mm$variable%like%'mult',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('optimal - estimated multiplicity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g5 <- ggplot(mm[mm$variable%like%'sens',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(0,1) +
        ylab('sensitivity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_brewer(palette = colpal)
g6 <- ggplot(mm[mm$variable%like%'spec',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(0,1) +
        ylab('specificity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=7)) +
        scale_colour_brewer(palette = colpal)

grid.arrange(g1, g3, g5,
             g2, g4, g6, ncol=3)
```


```{r calc_all_metrics}
metrics <- NULL
for(p1 in seq(0.1,0.9,0.1)) {
    p2 <- 1 - p1
    minor <- min(p1, p2)
    major <- max(p1, p2)
    mix <- paste(p1, p2, sep='-')

    x <- sv_data[sv_data$mix%in%mix, ]
    metrics <- rbind(metrics, calc_metrics_3clus(x, mix, 'ccube_SV', type='sv'))

    x <- pyc_all[pyc_all$mix%in%mix, ]
    metrics <- rbind(metrics, calc_metrics_3clus(x, mix, 'pyclone', type='pyc'))

    x <- snv_data[snv_data$mix%in%mix, ]
    metrics <- rbind(metrics, calc_metrics_3clus(x, mix, 'ccube_SNV', type='snv'))

    x <- cna_data[cna_data$mix%in%mix,]
    mean_ccf_error <- mean(as.numeric(x$ccf) - x$true_ccf, na.rm=T)
    is_subclonal_truth <- x$true_ccf < 1
    is_subclonal_bb <- x$ccf < cutoff_cna
    is_subclonal_sensitivity <- sum(is_subclonal_truth & is_subclonal_bb) / sum(is_subclonal_truth)
    is_subclonal_specificity <- 1 - sum(!is_subclonal_truth & is_subclonal_bb) / nrow(x)
    metrics <- rbind(metrics, data.frame(mix, mean_mult_error=NA, mean_ccf_error, clus_num_error=NA,
                                        is_subclonal_sensitivity, is_subclonal_specificity,
                                        clus_ccf_error=NA, method='battenberg'))
}

for(i in c('4clus', '5clus')) {
    if(i == '4clus') {
        cn <- read.delim(paste0(bb_out, '4clus_subclones.txt'))[,1:13]
    } else {
        cn <- read.delim(paste0(bb_out, '5clus_subclones.txt'))[,1:13]
    }
    cn <- cn[!((cn$frac1 == 1 & cn$nMaj1_A == 1 & cn$nMin1_A == 1) |
               (cn$chr %in% c('X','Y') &  cn$frac1 == 1 & cn$nMaj1_A == 1)), ]
    cx <- GRanges(seqnames = cn$chr, ranges = IRanges(start = cn$startpos, end = cn$endpos))
    olaps <- findOverlaps(cx, p1x, maxgap=boundary, type='equal')
    tmp <- cn[queryHits(olaps),]
    tmp$sample <- p1cn[subjectHits(olaps),]$sample
    tmp$subjectHits <- subjectHits(olaps)

    tmp1 <- tmp[grep('only',tmp$sample),]
    fracs <- apply(tmp1, 1,
                function(x){
                    maj <- p1cn[x['subjectHits'],'major']
                    min <- p1cn[x['subjectHits'],'minor']
                    cmaj1 <- as.numeric(x['nMaj1_A'])
                    cmin1 <- as.numeric(x['nMin1_A'])
                    cmaj2 <- as.numeric(x['nMaj2_A'])
                    cmin2 <- as.numeric(x['nMin2_A'])
                    if(maj == cmaj1 & min == cmin1){
                        return(x['frac1_A'])
                    } else if(!is.na(cmaj2) & (maj == cmaj2 & min == cmin2)) {
                        return(x['frac2_A'])
                    } else {return(NA)}
                })
    tmp$true_ccf <- 1
    if(i == '4clus') {
        tmp[tmp$sample%in%'bM_only' & tmp$chr %% 2 == 1, 'true_ccf'] <- 0.2
        tmp[tmp$sample%in%'bM_only' & tmp$chr %% 2 == 0, 'true_ccf'] <- 0.6
        tmp[tmp$sample%in%'gM_only' & tmp$chr %% 2 == 1, 'true_ccf'] <- 0.4
    } else {
        tmp[tmp$sample%in%'bM_only' & tmp$chr %% 2 == 1, 'true_ccf'] <- 0.8
        tmp[tmp$sample%in%'bM_only' & tmp$chr %% 2 == 0, 'true_ccf'] <- 0.6
        tmp[tmp$sample%in%'gM_only' & tmp$chr %% 2 == 1, 'true_ccf'] <- 0.2
        tmp[tmp$sample%in%'gM_only' & tmp$chr %% 2 == 0, 'true_ccf'] <- 0.4
    }
    tmp$ccf <- tmp$frac1_A
    tmp[grep('only',tmp$sample), 'ccf'] <- fracs
    x <- tmp[!is.na(tmp$ccf),]

    mean_ccf_error <- mean(as.numeric(x$ccf) - x$true_ccf, na.rm=T)
    is_subclonal_truth <- x$true_ccf < 1
    is_subclonal_bb <- x$ccf < cutoff_cna
    is_subclonal_sensitivity <- sum(is_subclonal_truth & is_subclonal_bb) / sum(is_subclonal_truth)
    is_subclonal_specificity <- 1 - sum(!is_subclonal_truth & is_subclonal_bb) / nrow(x)
    metrics <- rbind(metrics, data.frame(mix=i, mean_mult_error=NA, mean_ccf_error, clus_num_error=NA,
                                        is_subclonal_sensitivity, is_subclonal_specificity,
                                        clus_ccf_error=NA, method='battenberg'))
}

x <- load_hiclus_svs('4clus', basedir)
metrics <- rbind(metrics, calc_metrics_hiclus(x, '4clus', nclus=4, type='svs', truth='ground'))

x <- load_hiclus_svs('5clus', basedir)
metrics <- rbind(metrics, calc_metrics_hiclus(x, '5clus', nclus=5, type='svs', truth='ground'))

x <- load_hiclus_snvs('4clus', basedir)
metrics <- rbind(metrics, calc_metrics_hiclus(x, '4clus', nclus=4, type='snvs', truth='ground'))

x <- load_hiclus_snvs('5clus', basedir)
metrics <- rbind(metrics, calc_metrics_hiclus(x, '5clus', nclus=5, type='snvs', truth='ground'))

x <- pyc_all[pyc_all$mix%in%'4clus',]
metrics <- rbind(metrics, calc_metrics_hiclus(x, '4clus', nclus=4, type='pyc', truth='ground'))

x <- pyc_all[pyc_all$mix%in%'5clus',]
metrics <- rbind(metrics, calc_metrics_hiclus(x, '5clus', nclus=5, type='pyc', truth='ground'))

all_metrics <- metrics
```

## Figure 3

Performance across SVclone (SVs and SNVs), battenberg and PyClone.

```{r Figure3, fig.height=8, fig.width=12}
mm <- melt(metrics)
mm$variable <- as.character(mm$variable)
mm[mm$variable%like%'mean_mult_error',]$variable <- 'mean multiplicity error'
mm[mm$variable%like%'mean_ccf',]$variable <- 'mean variant CCF error'
mm[mm$variable%like%'sens',]$variable <- 'is subclonal sensitivity'
mm[mm$variable%like%'spec',]$variable <- 'is subclonal specificity'
mm[mm$variable%like%'clus_ccf',]$variable <- 'mean cluster CCF error'
mm[mm$variable%like%'clus_num',]$variable <- 'cluster number error'

mixes <- sapply(mm$mix[grep('-',mm$mix)],
                function(x){mix <- strsplit(x, '-')[[1]];
                   paste(as.numeric(mix[1])*100, as.numeric(mix[2])*100, sep='-') })
mm$mix[grep('-',mm$mix)] <- as.character(mixes)
mm$mix <- factor(mm$mix, levels=c('10-90', '20-80', '30-70', '40-60', '50-50',
                                  '60-40', '70-30', '80-20', '90-10', '4clus', '5clus'))

# remove pyclone as technically it doesn't estimate mult
mm <- mm[!(mm$method%in%'pyclone' & mm$variable%in%'mean multiplicity error'),]

cols <- brewer.pal(4, 'RdBu')
names(cols) <- unique(mm$method)

g1 <- ggplot(mm[mm$variable%like%'cluster num',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-3, 3) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated cluster number') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=7)) +
        scale_colour_manual(values = cols)
g2 <- ggplot(mm[mm$variable%like%'cluster CCF',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated cluster CCF') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=7)) +
        scale_colour_manual(values = cols)
g3 <- ggplot(mm[mm$variable%like%'mean var',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated CCF') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=7)) +
        scale_colour_manual(values = cols)
g4 <- ggplot(mm[mm$variable%like%'mult',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('optimal - estimated multiplicity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=7)) +
        scale_colour_manual(values = cols)
g5 <- ggplot(mm[mm$variable%like%'sens',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(0,1) +
        ylab('sensitivity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=7)) +
        scale_colour_manual(values = cols)
g6 <- ggplot(mm[mm$variable%like%'spec',], aes(mix, value, group=method, colour=method)) +
        geom_line() + geom_point() + facet_grid(~variable) + theme_bw() + ylim(0,1) +
        ylab('specificity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=7)) +
        scale_colour_manual(values = cols)

grid.arrange(g1, g3, g5,
             g2, g4, g6, ncol=3)
```

## Supplementary Figure 8

Performance of post-assignment using a joint SV + SNV model, compared with SV and SNV results wihout post-assign. PyClone results added for reference.

```{r SuppFigure8, fig.height=8, fig.width=12}
metrics <- NULL
for(p1 in seq(0.1,0.9,0.1)) {
    p2 <- 1 - p1
    mix <- paste(p1, p2, sep='-')
    svs <- sv_data[sv_data$mix%in%mix, ]
    metrics <- rbind(metrics, calc_metrics_3clus(svs, mix, 'SVs', type = 'sv'))

    snvs <- snv_data[snv_data$mix%in%mix, ]
    metrics <- rbind(metrics, calc_metrics_3clus(snvs, mix, 'SNVs', type = 'snv'))

    snvs <- pyc_all[pyc_all$mix%in%mix, ]
    metrics <- rbind(metrics, calc_metrics_3clus(snvs, mix, 'PyClone', type = 'pyc'))
}

for(i in seq(0.1, 0.9, 0.1)) {
    mix <- paste(i, 1-i, sep='-')
    basename <- paste('001bM_p', i, '_001gM_p', 1-i, sep='')

    sv_filt_file <- paste0(resultFolder, '/', basename, '/', basename, '_filtered_svs.tsv')
    pp_file <- paste0(resultFolder, '/', basename, '/purity_ploidy.txt')

    ccubeResultSVRDataFile = paste0(resultFolder, '/', basename, '/ccube_out/ccube_sv_results.RData')
    load(ccubeResultSVRDataFile)

    mix_name <- paste(i * 100, (1-i) * 100, sep='-')
    ccubeResultSNVRDataFile = paste0(clonal_basedir, basename, '/ccube_out/snvs/',
                                     basename, '_ccube_snv_results.RData')
    load(ccubeResultSNVRDataFile)

    referenceRes <- snvRes$res
    mydata <- doubleBreakPtsRes$ssm

    postAssignRes <- RunPostAssignPipeline(snvRes = referenceRes, svRes = doubleBreakPtsRes$res, mydata = mydata, verbose = F)
    svs <- get_all_sv_meta_3clus(postAssignRes$ssm, i, resultFolder, wd_truth)

    ssm <- snvRes$ssm
    postAssignRes <- RunPostAssignPipeline(svRes = postAssignRes$res, mydata = ssm, verbose = F)

    # combine SVs and SNVs
    snvs <- get_all_snv_meta_3clus(ssm = postAssignRes$ssm, i, resultFolder, wd_truth)

    sv_ccf <- apply(svs[,c('ccube_ccf1', 'ccube_ccf2')], 1, mean)
    ccube_mult <- apply(svs[,c('ccube_mult1', 'ccube_mult2')], 1, mean)
    true_cn <- apply(svs[,c('true_cn1', 'true_cn2')], 1, mean)
    svd <- data.frame(ccube_ccf=sv_ccf, ccube_mult=ccube_mult,
                      true_cn=true_cn, true_ccf=svs$true_ccf,
                      ccube_ccf_mean=svs$ccube_ccf_mean)
    snd <- data.frame(ccube_ccf=snvs$ccube_ccf, ccube_mult=snvs$ccube_mult,
                      true_cn=snvs$true_cn, true_ccf=snvs$true_ccf,
                      ccube_ccf_mean=snvs$ccube_ccf_mean)
    both <- rbind(svd, snd)
    metrics <- rbind(metrics, calc_metrics_3clus(both, mix, 'SV + SNV post-assigned', type='snv'))
}

mm <- melt(metrics)
mm$mix <- sapply(mm$mix, function(x){tmp <- as.numeric(strsplit(x, '-')[[1]]); paste(tmp[1]*100, tmp[2]*100, sep='-')})
mm$variable <- as.character(mm$variable)
mm[mm$variable%like%'mean_mult_error',]$variable <- 'mean multiplicity error'
mm[mm$variable%like%'mean_ccf',]$variable <- 'mean variant CCF error'
mm[mm$variable%like%'sens',]$variable <- 'is subclonal sensitivity'
mm[mm$variable%like%'spec',]$variable <- 'is subclonal specificity'
mm[mm$variable%like%'clus_ccf',]$variable <- 'mean cluster CCF error'
mm[mm$variable%like%'clus_num',]$variable <- 'cluster number error'

# remove pyclone as technically it doesn't estimate mult
mm <- mm[!(mm$method%in%'PyClone' & mm$variable%in%'mean multiplicity error'),]

cols <- c("#E41A1C", "#377EB8", "#ff7f00", "#984EA3")
names(cols) <- unique(mm$method)

g1 <- ggplot(mm[mm$variable%like%'cluster num',], aes(mix, value, group=method, colour=method)) +
        geom_line(alpha=0.7) + geom_point(alpha=0.7) + facet_grid(~variable) + theme_bw() + ylim(-3, 3) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated cluster number') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_manual(values = cols)
g2 <- ggplot(mm[mm$variable%like%'cluster CCF',], aes(mix, value, group=method, colour=method)) +
        geom_line(alpha=0.7) + geom_point(alpha=0.7) + facet_grid(~variable) + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated cluster CCF') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_manual(values = cols)
g3 <- ggplot(mm[mm$variable%like%'mean var',], aes(mix, value, group=method, colour=method)) +
        geom_line(alpha=0.7) + geom_point(alpha=0.7) + facet_grid(~variable) + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('true - estimated CCF') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_manual(values = cols)
g4 <- ggplot(mm[mm$variable%like%'mult',], aes(mix, value, group=method, colour=method)) +
        geom_line(alpha=0.7) + geom_point(alpha=0.7) + facet_grid(~variable) + theme_bw() + ylim(-0.5, 0.5) +
        geom_hline(yintercept=0, colour='grey') + ylab('optimal - estimated multiplicity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_manual(values = cols)
g5 <- ggplot(mm[mm$variable%like%'sens',], aes(mix, value, group=method, colour=method)) +
        geom_line(alpha=0.7) + geom_point(alpha=0.7) + facet_grid(~variable) + theme_bw() + ylim(0,1) +
        ylab('sensitivity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
        scale_colour_manual(values = cols)
g6 <- ggplot(mm[mm$variable%like%'spec',], aes(mix, value, group=method, colour=method)) +
        geom_line(alpha=0.7) + geom_point(alpha=0.7) + facet_grid(~variable) + theme_bw() + ylim(0,1) +
        ylab('specificity') + xlab('') +
        theme(legend.position='bottom', axis.text.x = element_text(size=7)) +
        scale_colour_manual(values = cols)

grid.arrange(g1, g3, g5,
             g2, g4, g6, ncol=3)
```

## Supplementary Figure 9

PyClone precision traces.

```{r SuppFigure9, fig.height=6, fig.width=12}
cols <- brewer.pal(9, 'Set1')
precisions <- NULL
for(i in seq(0.1, 0.9, 0.1)) {
    x <- read.delim(gzfile(paste0(base_pyc_results, "/001bM_p", i, "/trace/precision.tsv.bz2")), header=F)
    precisions <- rbind(precisions, data.frame(iter=as.numeric(rownames(x)),
                                               trace=as.numeric(x$V1), mix=as.character(i)))
}

tsize=20
ggplot(precisions, aes(iter, trace, group=mix, colour=mix)) +
    geom_line() + theme_bw() + xlab('iteration') + ylab('PyClone precision trace') +
    theme(axis.text = element_text(size=tsize),
          axis.title = element_text(size=tsize),
          legend.text = element_text(size=tsize),
          legend.title = element_text(size=tsize))
```



```{r run_cn_perturbation_experiments}
outdir <-  '001_svclone_results/001_CNdiv/001bM_p0.7_001gM_p0.3/'

het_frac <- 0.5
frac_divergence <- 0.3
magnitude <- 1

p1 <- i <- 0.7; p2 <- 1 - p1; mix <- paste(p1, p2, sep='-')
basename <- paste('001bM_p', p1, '_001gM_p', p2, sep='')

# uncomment below to generate data
# SVclone will then need to be run on output
#
# svin <- read.delim(paste0(basedir, basename, '/', basename, '_filtered_svs.tsv'), sep='\t')
# het_cns <- table(svin$gtype1==svin$gtype2)[1]
# non_het_cns <- as.numeric(round(het_cns/het_frac) - het_cns)
# SEED <- 123; set.seed(SEED); svin <- svin[sample(nrow(svin)),] #shuffle
# svin <- rbind(svin[svin$gtype1==svin$gtype2,][1:non_het_cns,],
#               svin[svin$gtype1!=svin$gtype2,])
# 
# # generate reduced normal set
# system(paste('mkdir -p', outdir))
# system(paste('mkdir -p', paste0(outdir, basename)))
# write.table(svin, file=paste0(outdir, basename, '/', basename, '_filtered_svs.tsv'), quote=F, sep='\t', row.names=F)
# pp_file <- paste0(basedir, basename, '/purity_ploidy.txt')
# system(paste('cp', pp_file, paste0(outdir, basename, '/purity_ploidy.txt')))
# 
# generate_de_cn_experiment(svin, basedir, basename, outdir, magnitude, two_sided=FALSE)
# generate_de_cn_experiment(svin, basedir, basename, outdir, magnitude, two_sided=TRUE)
# generate_de_cn_experiment(svin, basedir, basename, outdir, magnitude * -1, two_sided=FALSE)
# generate_de_cn_experiment(svin, basedir, basename, outdir, magnitude * -1, two_sided=TRUE)
# generate_de_cn_experiment(svin, basedir, basename, outdir, frac_divergence, two_sided=FALSE, frac=TRUE)
# generate_de_cn_experiment(svin, basedir, basename, outdir, frac_divergence, two_sided=TRUE, frac=TRUE)
# 
# generate_se_cn_experiment(svin, basedir, basename, outdir, 0) #normals
# generate_se_cn_experiment(svin, basedir, basename, outdir, frac_divergence, frac=TRUE)
# generate_se_cn_experiment(svin, basedir, basename, outdir, magnitude * -1)
# generate_se_cn_experiment(svin, basedir, basename, outdir, magnitude)
```

```{r load_cn_perturbation_experiments}
svd_data <- NULL
metrics <- NULL

svs <- get_svde_data(basedir, outdir, basename)
svs$experiment <- 'dual-end'; svd_data <- rbind(svs, svd_data)
metrics <- calc_metrics_3clus(svs, mix, 'dual-end', type='sv')
metrics$side <- 'unmodified'
metrics$model <- 'dual'
metrics$type <- 'unmodified'

for(side in c('oneside', 'twoside')) {
    for(type in c('minusCN', 'plusCN', 'frac')) {
        outdir <- paste0('001_svclone_results/001_CNdiv/', paste(basename, side, type, sep='_'))
        svs <- get_svde_data(basedir, outdir, basename)
        
        type_label <- ifelse(type=='minusCN', 'CN - 1', 'CN + 1')
        type_label <- ifelse(type=='frac', 'frac ± 0.3', type_label)
        side_label <- ifelse(side=='oneside', '(one side)', '(both sides)')
        label <- paste(type_label, side_label)
        svs$experiment <- label
        
        metric <- calc_metrics_3clus(svs, mix, label, type='sv')
        metric$side <- side; metric$model <- 'dual'; metric$type <- type_label
        metrics <- rbind(metrics, metric)
        
        svd_data <- rbind(svs, svd_data)
    }
}

svs_data <- NULL
for(side in c('side1', 'side2')) {
    for(type in c('minusCN', 'plusCN', 'frac', '')) {
        dirname <- ifelse(type!='', paste(basename, side, type, sep='_'), paste(basename, side, sep='_'))
        outdir <- paste0('001_svclone_results/001_CNdiv/', dirname)
        svs <- get_svse_data(basedir, outdir, basename)
        
        type_label <- ifelse(type=='minusCN', 'CN - 1', 'CN + 1')
        type_label <- ifelse(type=='frac', 'frac ± 0.3', type_label)
        type_label <- ifelse(type=='', '', type_label)
        label <- ifelse(type=='', paste('single-end', side), paste(type_label, side))
        svs$experiment <- label
        
        metric <- calc_metrics_3clus(svs, mix, label, type='snv')
        metric$side <- side; metric$model <- 'single'; metric$type <- ifelse(type=='', 'unmodified', type_label)
        metrics <- rbind(metrics, metric)
        
        svs_data <- rbind(svs, svs_data)
    }
}
```


## Supplementary Figure 11

Comparison of clustering metrics in the unmodified dual-end model under various background copy-number perturbations (+1 and -1 from major allele copy-number, and +0.3 or -0.3 from the copy-number fraction (for subclonal copy-numbers). 

```{r SuppFigure11, fig.height=7, fig.width=11}
mm <- melt(metrics)
mm$variable <- as.character(mm$variable)
mm[mm$variable%like%'mean_mult_error',]$variable <- 'mean multiplicity error'
mm[mm$variable%like%'mean_ccf',]$variable <- 'mean variant CCF error'
mm[mm$variable%like%'sens',]$variable <- 'is subclonal sensitivity'
mm[mm$variable%like%'spec',]$variable <- 'is subclonal specificity'
mm[mm$variable%like%'spec',]$variable <- 'is subclonal specificity'
mm[mm$variable%like%'clus_ccf',]$variable <- 'mean cluster CCF error'
mm[mm$variable%like%'clus_num',]$variable <- 'cluster number error'
mm$variable <- factor(mm$variable, levels=c('cluster number error', 'mean variant CCF error',
                                            'is subclonal specificity', 'mean cluster CCF error',
                                            'mean multiplicity error', 'is subclonal sensitivity'))

mm$experiment <- 'CN CCF deviation'
mm[mm$type%like%'CN', 'experiment'] <- 'CN deviation'
mm[mm$type%in%'unmodified', 'experiment'] <- 'unmodified'

mm$side <- 'unmodified'
mm$side[mm$method%like%'one'] <- 'one side'
mm$side[mm$method%like%'two'] <- 'both sides'
mm$side[mm$method%like%'side1'] <- 'side 1'
mm$side[mm$method%like%'side2'] <- 'side 2'

dual <- mm[mm$model%like%'dual',]
dual$method <- factor(dual$method, levels=c('frac ± 0.3 (both sides)', 'frac ± 0.3 (one side)',
                                            'CN - 1 (both sides)', 'CN - 1 (one side)',
                                            'CN + 1 (both sides)', 'CN + 1 (one side)', 'dual-end'))

cols <- brewer.pal(3, 'Dark2')
names(cols) <- c('unmodified', 'CN CCF deviation', 'CN deviation')

ggplot(dual, aes(method, value, group=method, fill=experiment)) +
    geom_bar(stat='identity') + theme_bw() + 
    geom_hline(yintercept=0, colour='grey')+ xlab('') + ylab('') +
    theme(legend.position='bottom', axis.text.x = element_text(size=8), legend.title = element_blank()) +
    scale_fill_manual(values = cols) + facet_wrap(~variable, scales='free_x') +
    coord_flip() + geom_hline(data=dual[dual$method=='dual-end',], aes(yintercept=value), colour='grey')
```

## Supplementary Figure 12

Comparison of clustering metrics for the dual and single-end models under copy-number perturbation.

```{r SuppFigure12, fig.height=7, fig.width=11}
single <- mm[mm$model%like%'single',]
single <- rbind(single, mm[mm$method=='dual-end',])
single$experiment <- as.character(sapply(as.character(single$method), function(x){strsplit(x, " side")[[1]][1]}))

x <- mm[(mm$model=='dual' & mm$side=='one side') | mm$model=='single' | mm$type=='unmodified',]
x <- data.table(x)[,mean(value),by=c('type', 'variable', 'model')]
x$type <- factor(x$type, levels=c('frac ± 0.3', 'CN + 1', 'CN - 1', 'unmodified'))

ggplot(x, aes(type, V1,  fill=model)) +
    geom_bar(stat='identity', position='dodge') + theme_bw() + 
    geom_hline(yintercept=0, colour='grey')+ xlab('') + ylab('') +
    theme(legend.position='bottom', axis.text.x = element_text(size=8)) +
    scale_fill_brewer(palette = 'Dark2') + facet_wrap(~variable, scales='free_x') +
    coord_flip() + geom_hline(data=x[x$type%like%'unmodified',], aes(yintercept=V1), colour='grey')
```