library(dplyr)

library(plyr)

library(clonevol)

library(fishplot)

library(reshape)

####step1: Using mpileup to extract coverage for each site. (Prepared 13/08/18)

####step2: Prepare input for pyclone, using segments from ASCAT

####step3: Run Pyclone twice. The first run is to identify founding cluster, 

####       the second run adding --tumour_content based on the cellular prevalence for each biopsy (21/08/18)

###step4: using clonevol to build evolution tree and model (in this script) (22/08/18)

###       using Fishplot to visulaize

###P75, 4 biopsies, nFL

#T1:LN_right_neck

#T2:LN_right_inguinal

#T3:LN_left_supraclav

#T4:EN_right_subcutaneous

setwd("G:/FL_exome_final/fl_latest/11_pyclone/BED_bam/counts_new/pyclone_output_301018/")

pyclone_out<- read.table(file="ouput_pyclone/output_200X_i2/P75_200X_i2/tables/loci.tsv",sep="\t",header=T)

mutations<-cast(pyclone_out[,1:4], mutation_id~sample_id, value="cellular_prevalence")

id<-unique(pyclone_out[, c(1,3)] )

P_case<- merge(id, mutations,by="mutation_id")

vafs = data.frame(cluster=P_case$cluster_id,

                  T1_vaf=(P_case$`GCF0150-0075-T01`/2*100),

                  T2_vaf=(P_case$`GCF0150-0075-T02`/2*100),

                  T3_vaf=(P_case$`GCF0150-0075-T03`/2*100),

                  T4_vaf=(P_case$`GCF0150-0075-T04`/2*100),

                  stringsAsFactors=F)

vafs$cluster<- vafs$cluster+1

###needs manully check density plot to identify which cluster is founding cluster

##clonevol requires founding cluster=1

##in this case founding cluster is #1, 

max_id<- max(vafs$cluster)+1

vafs$cluster[vafs$cluster==1]<-max_id

vafs$cluster[vafs$cluster==3]<-1

vafs$cluster[vafs$cluster==max_id]<-3

samples<-c("P75_1","P75_2", "P75_3", "P75_4")

samples1<-c("P75_1\nPretreat\nLN_right_neck","P75_2\nRelapsed_1\nLN_right_inguinal", 

            "P75_3\nRelapsed_2\nLN_left_supraclav", "P75_4\nRelapsed_3\nEN_right_subcutaneous")

names(vafs)[2:5] = samples 

##step 2: run infer.clonal.models, run twice: 1. include all cluster. 

###########2. manual review density plot and exlcude cluster with small number of mutations

dir.create("./clonevol/P75_final")

##

##first: use all clusters, no consensus models

#vafs$P44_4[vafs$cluster==7]<-vafs$P44_4[vafs$cluster==7]-1

res = infer.clonal.models(variants=vafs, cluster.col.name="cluster", vaf.col.names=samples,

                          subclonal.test="bootstrap", subclonal.test.model="non-parametric",

                          founding.cluster=1,

                          cluster.center="mean", num.boots=1000,

                          min.cluster.vaf=0.01, sum.p=0.01, alpha=0.01)

res = infer.clonal.models(variants=vafs, cluster.col.name="cluster", vaf.col.names=samples,

                          subclonal.test="bootstrap", subclonal.test.model="non-parametric",

                          founding.cluster=1,  

                          cluster.center="mean", num.boots=1000,ignore.clusters = c(2,8:11),

                          min.cluster.vaf=0.01,

                          sum.p=0.01, alpha=0.01)

#Finding consensus models across samples...

#Found  0 consensus model(s)

#Found 0 consensus model(s)

##second: ignore cluster ... consensus model based on manual review

#removing the clusters

vafs_used<- subset(vafs, !cluster %in% c(2,8:11))

#####rename clusters, cluster identities must be contiguous integer starting at 1

vafs_used$cluster[vafs_used$cluster==7]<-2

res = infer.clonal.models(variants=vafs_used, cluster.col.name="cluster", vaf.col.names=samples,

                          subclonal.test="bootstrap", subclonal.test.model="non-parametric",

                          founding.cluster=1,

                          cluster.center="mean", num.boots=1000,

                          min.cluster.vaf=0.01, sum.p=0.01, alpha=0.01)

res<-convert.consensus.tree.clone.to.branch(res, branch.scale = 'sqrt')

pdf("./clonevol/P75_final/P75_trees.pdf", useDingbats = FALSE)

plot.all.trees.clone.as.branch(res, branch.width = 0.5, node.size = 1, node.label.size = 0.5)

dev.off()

plot.clonal.models(res,

                   # box plot parameters

                   box.plot = TRUE,

                   fancy.boxplot = TRUE,

                   fancy.variant.boxplot.highlight = 'is.driver',

                   fancy.variant.boxplot.highlight.shape = 21,

                   fancy.variant.boxplot.highlight.fill.color = 'red',

                   fancy.variant.boxplot.highlight.color = 'black',

                   fancy.variant.boxplot.highlight.note.col.name = 'gene',

                   fancy.variant.boxplot.highlight.note.color = 'blue',

                   fancy.variant.boxplot.highlight.note.size = 2,

                   fancy.variant.boxplot.jitter.alpha = 1,

                   fancy.variant.boxplot.jitter.center.color = 'grey50',

                   fancy.variant.boxplot.base_size = 12,

                   fancy.variant.boxplot.plot.margin = 1,

                   fancy.variant.boxplot.vaf.suffix = '.VAF',

                   # 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 = TRUE,

                   mtcab.event.sep.char = ',',

                   mtcab.branch.text.size = 1,

                   mtcab.branch.width = 0.75,

                   mtcab.node.size = 3,

                   mtcab.node.label.size = 1,

                   mtcab.node.text.size = 1.5,

                   # 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 = './clonevol/P75_final/',

                   out.format = 'pdf',

                   overwrite.output = TRUE,

                   width = 12,

                   height = 7,

                   # vector of width scales for each panel from left to right

                   panel.widths = c(3,4,2,4,2))

plot.clonal.models(res,

                   # box plot parameters

                   box.plot = TRUE,

                   fancy.boxplot = TRUE,

                   fancy.variant.boxplot.highlight = 'is.driver',

                   fancy.variant.boxplot.highlight.shape = 21,

                   fancy.variant.boxplot.highlight.fill.color = 'red',

                   fancy.variant.boxplot.highlight.color = 'black',

                   fancy.variant.boxplot.highlight.note.col.name = 'gene',

                   fancy.variant.boxplot.highlight.note.color = 'blue',

                   fancy.variant.boxplot.highlight.note.size = 2,

                   fancy.variant.boxplot.jitter.alpha = 1,

                   fancy.variant.boxplot.jitter.center.color = 'grey50',

                   fancy.variant.boxplot.base_size = 12,

                   fancy.variant.boxplot.plot.margin = 1,

                   fancy.variant.boxplot.vaf.suffix = '.VAF',

                   # 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 = TRUE,

                   mtcab.event.sep.char = ',',

                   mtcab.branch.text.size = 1,

                   mtcab.branch.width = 0.75,

                   mtcab.node.size = 3,

                   mtcab.node.label.size = 1,

                   mtcab.node.text.size = 1.5,

                   # 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 = TRUE,

                   # output figure parameters

                   out.dir = './clonevol/P75_final/',

                   out.format = 'pdf',

                   overwrite.output = TRUE,

                   width = 12,

                   height = 7,

                   # vector of width scales for each panel from left to right

                   panel.widths = c(3,4,2,4,2))

##generating fish plot

f<- generateFishplotInputs(results = res)

fishes=createFishPlotObjects(f)

pdf('./clonevol/P75_final/P75_fish_200x_pyclone_anno.pdf', width=14, height=7)

for (i in 1:length(fishes)){

  fish = layoutClones(fishes[[i]])

  fish = setCol(fish,f$clonevol.clone.colors)

  fishPlot(fish,shape="spline", title.btm="P1", cex.title=0.7,cex.vlab = 1.4,

           vlines=seq(1, length(samples1)), vlab=samples1, pad.left=0.5)

}

dev.off()

##other plots

pdf("./clonevol/P75_final/P75_box.pdf", width=3, height=3,useDingbats = FALSE, title='')

pp<-plot.variant.clusters(vafs_used,

                          cluster.col.name = 'cluster',

                          show.cluster.size = FALSE,

                          cluster.size.text.color = 'blue',

                          vaf.col.names = samples,

                          vaf.limits = 70,

                          sample.title.size = 20,

                          violin = FALSE,

                          box = FALSE,

                          jitter = TRUE,

                          jitter.shape = 1,

                          jitter.size = 3,

                          jitter.alpha = 1,

                          jitter.center.method = 'median',

                          jitter.center.size = 1,

                          jitter.center.color = 'darkgray',

                          jitter.center.display.value = 'none',

                          highlight = 'is.driver',

                          highlight.shape = 21,

                          highlight.color = 'blue',

                          highlight.fill.color = 'green',

                          highlight.note.col.name = 'gene',

                          highlight.note.size = 2,

                          order.by.total.vaf = FALSE)

dev.off()

plot.pairwise(vafs_used, col.names = samples,

              out.prefix = './clonevol/P75_final/P75_variants.pairwise.plot')

pdf('./clonevol/P75_final/P75_flow.pdf')

plot.cluster.flow(vafs_used, vaf.col.names = samples,

                  sample.names = samples1)

dev.off()

Pcase<-do.call("rbind", lapply( list.files("input_pyclone_271018_newpara/200X/GCF0150-0075-N01_200X/",full=TRUE),

                                read.table, header=TRUE, sep="\t"))

library(ggplot2)

pdf("clonevol/P75_final/coverage.pdf")

ggplot(Pcase, aes(x=(var_counts+ref_counts)))+

  geom_histogram(position="dodge")+

  facet_grid(~sample)

dev.off()

save.image("G:/FL_resequncing/FL_exome_final/fl_latest/11_pyclone/BED_bam/counts_new/pyclone_output_301018/P75.RData")