Main

Inferring the clonal structures by the Clonevol (https://github.com/hdng/clonevol).

Prepare data

library(clonevol)
library(ggrepel)
## Loading required package: ggplot2
library(tidyverse)
## Warning in system("timedatectl", intern = TRUE): running command 'timedatectl'
## had status 1
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1
## ✓ purrr   0.3.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
mutdata.vaf1 = readRDS("data/mutdata.vaf1.rds")
new.label = readRDS(file = "data/new.label.rds")



mutdata = mutdata.vaf1 %>%
  left_join(new.label)
## Joining, by = c("mutid", "sciClone", "kmeans")
mutdata[, 2:34] = mutdata[, 2:34]*100

mutdata$cluster = mutdata$Clone2

#Visualizing the variant clusters
set.colors = c("#C6C6C6", "#6FDCBF","#5AA36E","#E99692","#B4D985","#EA7D23","#E53CFB","#4B85B7","#8439FA","#BD8736","#B3B371","#A7C7DE","#EE97FC","#57C222","#BFABD0","#44589B", "#794C18", RColorBrewer::brewer.pal(n = 10, name = "Paired") )

clone.colors = set.colors[ 1:length(unique(mutdata$cluster))]
cluster.col.name  = "cluster"

head(mutdata)

Labeling the driver mutations.

#Stomach drivers

GC01 = readRDS("data/GC01.rds")

GC01 = GC01 %>% mutate(Chromosome =  str_c("chr", Chromosome) ) %>%
  mutate(mutid = str_c(Chromosome, Start_Position, Reference_Allele, sep = ":")) 


drivers = read.delim("data/IntOGen-Drivers-Cancer_Genes.tsv", header = T) %>%
  filter(CANCER_TYPE %in% c("BRCA","COREAD","LUAD", "LUSC") ) %>%
  pull(SYMBOL) %>% unique()

mutinfo = GC01 %>%
  select(Chromosome, Start_Position, End_Position, Hugo_Symbol, Variant_Classification, Protein_Change, Reference_Allele, Tumor_Seq_Allele2 ) %>%
  unique.data.frame() %>%
  mutate(#Chromosome = str_c("chr", Chromosome),
    Protein_Change = ifelse(Variant_Classification == "Splice_Site", "splice", Protein_Change),
    Protein_Change = ifelse(Variant_Classification == "3'UTR", "3UTR", Protein_Change),
    Protein_Change = ifelse(Variant_Classification == "5'UTR", "5UTR", Protein_Change),
    gene_site = str_c(Hugo_Symbol, "_" ,Protein_Change),
    is.driver = ifelse( (Hugo_Symbol %in% drivers) & ( !Variant_Classification %in% c("Silent", "3'Flank", "IGR", "Intron", "RNA") ) , TRUE, FALSE)  ) %>%
  
  mutate( mutid  = str_c(Chromosome , Start_Position , Reference_Allele , sep = ":") )


mutdata = right_join(mutinfo, mutdata ) %>%
  arrange(cluster) %>%
  as.data.frame() 
## Joining, by = "mutid"
#add ref and var columns for merge data.

sampleNames = c( paste0("Breast_", 1:5), paste0("Lung_", 1:5),
                 paste0("Coad_", 1:5), paste0("OveryLM_", 1:5), 
                 paste0("OveryRM_", 1:6), paste0("UterusM_", c(1:7)) )

var.data = round(mutdata[, sampleNames])
ref.data = 100 - var.data

colnames(var.data) = paste0(colnames(var.data), ".var")
colnames(ref.data) = paste0(colnames(ref.data), ".ref")

mutdata = cbind(mutdata, var.data, ref.data)

The clonevol infers the clonal history

Visualzied the mutation clusters

  1. Plot the mutations in each cluster and each samples.
message("<2> Analysis using clonevol")
## <2> Analysis using clonevol
output = "Met2" #Output folder names

#saveRDS(mutinfo, file = sprintf("%s/%s.mutinfo.rds", output))

sampleNames = c( paste0("Coad_", 1:5), paste0("OveryLM_", 1:5), 
                 paste0("OveryRM_", 1:6),paste0("UterusM_", c(1, 3)) )


vaf.col.names = sampleNames

show.ccf = F

pp <- plot.variant.clusters(mutdata,  
                            show.cluster.size = F,
                            show.cluster.label= F,
                            cluster.col.name = cluster.col.name, 
                            vaf.col.names = vaf.col.names,
                            violin = FALSE,
                            box = FALSE,
                            jitter = TRUE,
                            jitter.shape = 1,
                            variant.class.col.name =cluster.col.name,
                            #vaf.limits = 
                            jitter.color = clone.colors,
                            jitter.size = 1.2,
                            jitter.alpha = 1,
                            jitter.width = 0.2,
                            jitter.center.method = "median",
                            jitter.center.size = 1,
                            jitter.center.color = "darkgray",
                            jitter.center.display.value = "none",
                            display.plot = F,
                            horizontal = T,
                            order.by.total.vaf = F, 
                            highlight = 'is.driver',
                            highlight.shape = 21,
                            highlight.color = 'blue',
                            highlight.fill.color = 'green',
                            highlight.size = 2.5,
                            
                            highlight.note.col.name = NULL,
                            highlight.note.size = 2,
                            highlight.note.color = "blue",
                            highlight.note.angle = 0,
                            founding.cluster = 1, 
                            ccf = show.ccf
)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
#if(show.ccf){max.limits = 180 }else{max.limits = 80}
max.limits = 110

pp[[1]] = pp[[1]] +
  theme(axis.title.y = element_blank()) +
  scale_y_continuous(position = "right", limits  = c(0, max.limits)) +
  theme(axis.text.x  = element_text(angle = 270, size = 8) )
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
for(ii in 2:length(pp)){
  pp[[ii ]] = pp[[ii]] +
    scale_y_continuous(position = "right", limits  = c(0, max.limits)) +
    theme(axis.text.x  = element_text(angle = 270, size = 8 ),
          axis.text.y = element_blank()
    )
}
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
#annotate the driver genes.

if(any(mutdata$is.driver) ){
  
  labels = pp[[1]]$data
  #select colors for annotation.
  clone.colors.sel =unique( clone.colors[ labels$cluster[labels$is.driver]])
  
  pp1 = labels %>% filter(is.driver) %>%
    mutate(cluster_1 = factor(cluster)) %>%
    ggplot(aes(y = cluster, x = 1, label = gene_site)) + theme_classic() +
    geom_point(aes( color = cluster_1), size = 3)+
    geom_text_repel(
      nudge_x      = 0.15,
      direction    = "y",
      hjust        = 0,
      segment.size = 0.2,
      size = 3,
    ) +
    ylim(0,  length(clone.colors) +1 ) +
    xlim(1, 0.8) +
    scale_color_manual( values= clone.colors.sel, guide="none") +
    theme(
      axis.line  = element_blank(),
      axis.ticks = element_blank(),
      axis.text  = element_blank(),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
      plot.title   = element_text(hjust = 0.5)
    )
  
  pp[[ length(pp) +1 ]] = pp1
  
}

pdf( sprintf( "%s/%s.box.pdf", output, output), width = 2*length(pp), height = 4)
ggpubr::ggarrange(plotlist = pp, ncol = length(pp), align  = "h" )
dev.off()
## svg 
##   2
ggpubr::ggarrange(plotlist = pp, ncol = length(pp), align  = "h" )

  1. Plotting mean/median of clusters across samples (cluster flow)
p1 = plot.cluster.flow(mutdata,
                  cluster.col.name = cluster.col.name,
                  vaf.col.names = vaf.col.names,
                  sample.names =vaf.col.names,
                  colors = clone.colors,
                  y.title = "Variant Allele Frequency %"
) +
  theme(axis.text.x = element_text(angle = 90) )
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
pdf( sprintf( "%s/%s.flow.pdf", output, output) , width=8, height=5 )
p1
dev.off()
## svg 
##   2
p1

Clonal Structure inference and visualization

This script first infers the clonal orders in each sample by the pigeonhole principle. Then get the consensus trees across samples. Upon the built clonal trees, the scripts plot the clonal evolutionary history in each sample, the whole clonal structures for each sample, and the inferred clonal trees in the tumor.

source("data/auxiliary_functions.R")

# Inferring the main structures of mutation clusters
sample.groups = mapply(function(x) x[1], strsplit(vaf.col.names, "_")  )
names(sample.groups) = vaf.col.names

sel = c(1:18) 
y = plot.tree( output  = output , mutdata = mutdata, 
               ccf.col.names = vaf.col.names[sel],
               sample.groups = sample.groups[sel],
               cancer.initiation.model = "monoclonal",
               founding.cluster = 1 , ignore.clusters = 4 , cluster.col.name = "cluster",
               subclonal.test.model = "non-parametric", 
               #scale.monoclonal.cell.frac = FALSE,
               sum.p = 0.01, alpha = 0.05, weighted = FALSE,
               plt.pairwise = F,
               highlight.note.col.name = NULL,
               highlight = "is.driver", highlight.CCF = TRUE
)
## infer.clonal.models
## Calculate VAF as CCF/2
## Sample 1: Coad_1 <-- Coad_1
## Sample 2: Coad_2 <-- Coad_2
## Sample 3: Coad_3 <-- Coad_3
## Sample 4: Coad_4 <-- Coad_4
## Sample 5: Coad_5 <-- Coad_5
## Sample 6: OveryLM_1 <-- OveryLM_1
## Sample 7: OveryLM_2 <-- OveryLM_2
## Sample 8: OveryLM_3 <-- OveryLM_3
## Sample 9: OveryLM_4 <-- OveryLM_4
## Sample 10: OveryLM_5 <-- OveryLM_5
## Sample 11: OveryRM_1 <-- OveryRM_1
## Sample 12: OveryRM_2 <-- OveryRM_2
## Sample 13: OveryRM_3 <-- OveryRM_3
## Sample 14: OveryRM_4 <-- OveryRM_4
## Sample 15: OveryRM_5 <-- OveryRM_5
## Sample 16: OveryRM_6 <-- OveryRM_6
## Sample 17: UterusM_1 <-- UterusM_1
## Sample 18: UterusM_3 <-- UterusM_3
## Using monoclonal model
## Note: all VAFs were divided by 100 to convert from percentage to proportion.
## Generating non-parametric boostrap samples...
## Coad_1 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 2,3,4,5,6,7,8 
## User ignored clusters:  4 
## Coad_1 : 2 clonal architecture model(s) found
## 
## Coad_2 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 3,4,5,6,7,8 
## User ignored clusters:  4 
## Coad_2 : 6 clonal architecture model(s) found
## 
## Coad_3 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 3,4,5,6,7,8 
## User ignored clusters:  4 
## Coad_3 : 6 clonal architecture model(s) found
## 
## Coad_4 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 2,3,4,5,6,7,8 
## User ignored clusters:  4 
## Coad_4 : 2 clonal architecture model(s) found
## 
## Coad_5 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 3,4,5,6,7,8,10 
## User ignored clusters:  4 
## Coad_5 : 2 clonal architecture model(s) found
## 
## OveryLM_1 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,6,7,9,10 
## User ignored clusters:  4 
## OveryLM_1 : 12 clonal architecture model(s) found
## 
## OveryLM_2 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,6,7,9,10 
## User ignored clusters:  4 
## OveryLM_2 : 10 clonal architecture model(s) found
## 
## OveryLM_3 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,6,7,9,10 
## User ignored clusters:  4 
## OveryLM_3 : 13 clonal architecture model(s) found
## 
## OveryLM_4 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,6,7,9,10 
## User ignored clusters:  4 
## OveryLM_4 : 12 clonal architecture model(s) found
## 
## OveryLM_5 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,6,7,9,10 
## User ignored clusters:  4 
## OveryLM_5 : 13 clonal architecture model(s) found
## 
## OveryRM_1 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,7,8,9,10 
## User ignored clusters:  4 
## OveryRM_1 : 13 clonal architecture model(s) found
## 
## OveryRM_2 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,7,8,9,10 
## User ignored clusters:  4 
## OveryRM_2 : 13 clonal architecture model(s) found
## 
## OveryRM_3 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,7,8,9,10 
## User ignored clusters:  4 
## OveryRM_3 : 12 clonal architecture model(s) found
## 
## OveryRM_4 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,7,8,9,10 
## User ignored clusters:  4 
## OveryRM_4 : 10 clonal architecture model(s) found
## 
## OveryRM_5 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,7,8,9,10 
## User ignored clusters:  4 
## OveryRM_5 : 12 clonal architecture model(s) found
## 
## OveryRM_6 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 4,5,7,8,9,10 
## User ignored clusters:  4 
## OveryRM_6 : 12 clonal architecture model(s) found
## 
## UterusM_1 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 2,3,5,6,7,8,9,10 
## User ignored clusters:  4 
## UterusM_1 : 1 clonal architecture model(s) found
## 
## UterusM_3 : Enumerating clonal architectures...
## Determining if cluster VAF is significantly positive...
## Exluding clusters whose VAF < min.cluster.vaf=0.01
## Non-positive VAF clusters: 2,3,5,6,7,8,9,10 
## User ignored clusters:  4 
## UterusM_3 : 1 clonal architecture model(s) found
## 
## Finding consensus models across samples...
## Found  6 consensus model(s)
## Generating consensus clonal evolution trees across samples...
## Found 6 consensus model(s)
## Scoring models...
## Pruning consensus clonal evolution trees....
## Seeding aware pruning is:  off
## Number of unique pruned consensus trees: 6
## convert.consensus.tree.clone.to.branch
## plot.clonal.models
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:dplyr':
## 
##     as_data_frame, groups, union
## The following objects are masked from 'package:purrr':
## 
##     compose, simplify
## The following object is masked from 'package:tidyr':
## 
##     crossing
## The following object is masked from 'package:tibble':
## 
##     as_data_frame
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union
## Plotting pruned consensus trees...
## Output plots are in: Met2
#PLot the clonal trees.

pdf(file = sprintf("%s/%s.trees.pdf", output, output), width = 6, height = 6)

plot.all.trees.clone.as.branch(y, branch.width = 0.5,
                               node.size = 2, node.label.size = 0.5,
                               tree.rotation    =180,
                               angle = 20
)

dev.off()
## svg 
##   2
plot.all.trees.clone.as.branch(y, branch.width = 0.5,
                               node.size = 2, node.label.size = 0.5,
                               tree.rotation    =180,
                               angle = 20
)

plot.clonal.models(y,
                   # box plot parameters
                   box.plot = FALSE,
                   fancy.boxplot = FALSE,
                   
                   # bell plot parameters
                   clone.shape = "bell",
                   bell.event = TRUE,
                   bell.event.label.color = "blue",
                   bell.event.label.angle = 60,
                   clone.time.step.scale = 1,
                   bell.curve.step = 2,
                   
                   # node-based consensus tree parameters
                   merged.tree.plot = FALSE,
                   tree.node.label.split.character = NULL,
                   tree.node.shape = "circle",
                   tree.node.size = 30,
                   tree.node.text.size = 0.5,
                   merged.tree.node.size.scale = 1.25,
                   merged.tree.node.text.size.scale = 2.5,
                   merged.tree.cell.frac.ci = FALSE,
                   
                   # branch-based consensus tree parameters
                   merged.tree.clone.as.branch = FALSE,
                   
                   # cellular population parameters
                   cell.plot = TRUE,
                   num.cells = 100,
                   cell.border.size = 0.25,
                   cell.border.color = "black",
                   clone.grouping = "horizontal",
                   
                   #meta-parameters
                   scale.monoclonal.cell.frac = TRUE,
                   show.score = FALSE,
                   cell.frac.ci = TRUE,
                   disable.cell.frac = FALSE,
                   
                   # output figure parameters
                   out.dir = sprintf("%s/%s", output, "merge"),
                   out.format = "pdf",
                   overwrite.output = TRUE,
                   width = 9,
                   #height = ,
                   # vector of width scales for each panel from left to right
                   panel.widths = c(4,2))
## Plotting pruned consensus trees...
## Output plots are in: Met2/merge
save(mutdata, clone.colors, cluster.col.name , output, sampleNames, y, sample.groups, file = sprintf("%s/clonevol.data.rda", output) )

Merge multi-samples into a meta sample.

Merge multiple samples into a single sample. This can be used to merge multiregional sequencing (MRS) into one sample that could represent whole tumors.

Note: CCF should be in the range between 0-100

#merge coad
sel = 1:5
y.merge = merge.samples(y, samples = vaf.col.names[sel], 
                        new.sample = "Coad", new.sample.group = "Coad", 
                        ref.cols = str_c(vaf.col.names[sel], ".ref"), 
                        var.cols = str_c(vaf.col.names[sel], ".var") )
## Note: all VAFs were divided by 100 to convert from percentage to proportion.
## Generating non-parametric boostrap samples...
## Estimating CCF of clones for merged sample...
## Merging clonal evolution trees across samples...
## WARN: Model 1 : Clone(s) 1 (existed in one of Coad_1+Coad_2+Coad_3+Coad_4+Coad_5 ) now lost in Coad 
## ** merged sample Coad reannotated
## WARN: Model 2 : Clone(s) 1 (existed in one of Coad_1+Coad_2+Coad_3+Coad_4+Coad_5 ) now lost in Coad 
## ** merged sample Coad reannotated
## WARN: Model 3 : Clone(s) 1 (existed in one of Coad_1+Coad_2+Coad_3+Coad_4+Coad_5 ) now lost in Coad 
## ** merged sample Coad reannotated
## WARN: Model 4 : Clone(s) 1 (existed in one of Coad_1+Coad_2+Coad_3+Coad_4+Coad_5 ) now lost in Coad 
## ** merged sample Coad reannotated
## WARN: Model 5 : Clone(s) 1 (existed in one of Coad_1+Coad_2+Coad_3+Coad_4+Coad_5 ) now lost in Coad 
## ** merged sample Coad reannotated
## WARN: Model 6 : Clone(s) 1 (existed in one of Coad_1+Coad_2+Coad_3+Coad_4+Coad_5 ) now lost in Coad 
## ** merged sample Coad reannotated
## Pruning consensus clonal evolution trees....
## Seeding aware pruning is:  off
## Number of unique pruned consensus trees: 6
#merge overyLM
sel = 6:10
y.merge = merge.samples(y.merge, samples = vaf.col.names[sel], 
                        new.sample = "OveryLM", new.sample.group = "OveryLM", 
                        ref.cols = str_c(vaf.col.names[sel], ".ref"), 
                        var.cols = str_c(vaf.col.names[sel], ".var") )
## Note: all VAFs were divided by 100 to convert from percentage to proportion.
## Generating non-parametric boostrap samples...
## Estimating CCF of clones for merged sample...
## Merging clonal evolution trees across samples...
## WARN: Model 1 : Clone(s) 2 (existed in one of OveryLM_1+OveryLM_2+OveryLM_3+OveryLM_4+OveryLM_5 ) now lost in OveryLM 
## ** merged sample OveryLM reannotated
## WARN: Model 2 : Clone(s) 2 (existed in one of OveryLM_1+OveryLM_2+OveryLM_3+OveryLM_4+OveryLM_5 ) now lost in OveryLM 
## ** merged sample OveryLM reannotated
## WARN: Model 3 : Clone(s) 2 (existed in one of OveryLM_1+OveryLM_2+OveryLM_3+OveryLM_4+OveryLM_5 ) now lost in OveryLM 
## ** merged sample OveryLM reannotated
## WARN: Model 5 : Clone(s) 2 (existed in one of OveryLM_1+OveryLM_2+OveryLM_3+OveryLM_4+OveryLM_5 ) now lost in OveryLM 
## ** merged sample OveryLM reannotated
## Pruning consensus clonal evolution trees....
## Seeding aware pruning is:  off
## Number of unique pruned consensus trees: 3
#merge overyRM
sel = c(11:13, 14, 15, 16)
y.merge = merge.samples(y.merge, samples = vaf.col.names[sel], 
                        new.sample = "OveryRM", new.sample.group = "OveryRM", 
                        ref.cols = str_c(vaf.col.names[sel], ".ref"), 
                        var.cols = str_c(vaf.col.names[sel], ".var") )
## Note: all VAFs were divided by 100 to convert from percentage to proportion.
## Generating non-parametric boostrap samples...
## Estimating CCF of clones for merged sample...
## Merging clonal evolution trees across samples...
## Pruning consensus clonal evolution trees....
## Seeding aware pruning is:  off
## Number of unique pruned consensus trees: 2
#merge UterusM
sel = c(17, 18)
y.merge = merge.samples(y.merge, samples = vaf.col.names[sel], 
                        new.sample = "UterusM", new.sample.group = "UterusM", 
                        ref.cols = str_c(vaf.col.names[sel], ".ref"), 
                        var.cols = str_c(vaf.col.names[sel], ".var") )
## Note: all VAFs were divided by 100 to convert from percentage to proportion.
## Generating non-parametric boostrap samples...
## Estimating CCF of clones for merged sample...
## Merging clonal evolution trees across samples...
## Pruning consensus clonal evolution trees....
## Seeding aware pruning is:  off
## Number of unique pruned consensus trees: 2
plot.clonal.models(y.merge,
                   # box plot parameters
                   box.plot = FALSE,
                   fancy.boxplot = FALSE,
                   
                   # bell plot parameters
                   clone.shape = "bell",
                   bell.event = TRUE,
                   bell.event.label.color = "blue",
                   bell.event.label.angle = 60,
                   clone.time.step.scale = 1,
                   bell.curve.step = 2,
                   
                   # node-based consensus tree parameters
                   merged.tree.plot = TRUE,
                   tree.node.label.split.character = NULL,
                   tree.node.shape = "circle",
                   tree.node.size = 30,
                   tree.node.text.size = 0.5,
                   merged.tree.node.size.scale = 1.25,
                   merged.tree.node.text.size.scale = 2.5,
                   merged.tree.cell.frac.ci = FALSE,
                   
                   # branch-based consensus tree parameters
                   merged.tree.clone.as.branch = FALSE,
                   
                   # cellular population parameters
                   cell.plot = TRUE,
                   num.cells = 100,
                   cell.border.size = 0.25,
                   cell.border.color = "black",
                   clone.grouping = "horizontal",
                   
                   #meta-parameters
                   scale.monoclonal.cell.frac = TRUE,
                   show.score = FALSE,
                   cell.frac.ci = TRUE,
                   disable.cell.frac = FALSE,
                   
                   # output figure parameters
                   out.dir = sprintf("%s/%s", output, "merge1"),
                   out.format = "pdf",
                   overwrite.output = TRUE,
                   width = 9,
                   #height = ,
                   # vector of width scales for each panel from left to right
                   panel.widths = c(4,4,4))
## Plotting pruned consensus trees...
## Output plots are in: Met2/merge1
#merge MRS data into one.
save(y.merge , file = sprintf("%s/%s.y.merge.rda", output, output))