Datasets
Models
Downloads: 1
Diff of
/Figures/clonevol/clonevol_P75_021118.R
[000000]
..
[4ff317]
Switch to side-by-side view
--- a +++ b/Figures/clonevol/clonevol_P75_021118.R @@ -0,0 +1,312 @@ +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")