library(dplyr)

library(plyr)

library(clonevol)

library(fishplot)

library(reshape)

####coverage: from reseq BAM files 

####step1: Using mpileup to extract coverage for each site. (Prepared 27/10/18), default parameter

####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 (29/10/18)

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

###       using Fishplot to visulaize

###P21, 2 biopies, nFL

#T1:LN_neck_unknown

#T2:LN_left_inguinal

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

pyclone_out<- read.table(file="ouput_pyclone/output_200X_i2/P21_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-0021-T01`/2)*100,

                  T2_vaf=(P_case$`GCF0150-0021-T02`/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

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("P21_1","P21_2")

samples1<-c("P21_1\nPretreat1","P21_2\nPretreat2")

names(vafs)[2:3] = 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/P21")

##

##first: use all clusters, no consensus models

vafs$P21_2[vafs$cluster=="2"]<-vafs$P21_2[vafs$cluster=="2"]-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,ignore.clusters = c(5,6,8:10),

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

                          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 5:7 consensus model based on manual review

vafs_used<- subset(vafs, !cluster %in% c(5,6,8:10))

vafs_used$cluster[vafs_used$cluster=="7"]<-5

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/P21/P21_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,

                   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/P21/',

                   out.format = 'pdf',

                   overwrite.output = TRUE,

                   width = 10,

                   height = 4,

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

                   panel.widths = c(1.5,2.5,1.5,2.5,2))

###removing cell.frac annotation

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,

                   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/P21/',

                   out.format = 'pdf',

                   overwrite.output = TRUE,

                   width = 10,

                   height = 4,

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

                   panel.widths = c(1.5,2.5,1.5,2.5,2))

##generating fish plot

f<- generateFishplotInputs(results = res)

fishes=createFishPlotObjects(f)

samples2<-c("P21_1\nPrimary\nLN_neck_unknown","P21_2\nPretreat\nLN_left_inguinal")

pdf('./clonevol/P21/P21_fish_200x_pyclone_anno_loc.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(samples2)), vlab=samples2, pad.left=0.5)

}

dev.off()

pdf("./clonevol/P21/P21_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/P21/P21_variants.pairwise.plot')

pdf('./clonevol/P21/P21_flow.pdf')

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

                  sample.names = c('Pretreat1', 'Pretreat2'))

dev.off()

####checking coverage f

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

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

#P21_1<- read.table(file="input_pyclone_271018_newpara/200X/GCF0150-0021-N01_200X/GCF0150-0021-T01.txt",sep="\t",header=T)

#P21_2<- read.table(file="input_pyclone_271018_newpara/200X/GCF0150-0021-N01_200X/GCF0150-0021-T02.txt",sep="\t",header=T)

#Pcase<- rbind(P21_1, P21_2)

########

#min (dp) =min(c1inP4_1$var_counts+c1inP4_1$ref_counts)=273

#max (dp) =max(c1inP4_1$var_counts+c1inP4_1$ref_counts)=648

#median (dp) =median(c1inP4_1$var_counts+c1inP4_1$ref_counts)=380

#mean (dp) =mean(c1inP4_1$var_counts+c1inP4_1$ref_counts)=417

library(ggplot2)

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

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

  geom_histogram(position="dodge")+

  facet_grid(~sample)

dev.off()