################################################

#  File Name:Figure3.R

#  Author: Baoyan Bai

#################################################

####Integrate three different mutation files:

#1. original mutation calling file from Sigve

#   this files has been filtered by Ole-Christian'script and by Jillian's TCGA PON 

#2. raw vcf files from TSD. Many mutations in serial biopsies were in fact called, but rejcted for different reasons. 

#   I started to work this. Noticed some of mutations were still not called. But if looked at coverage, variant allele can be detected

#3. Pyclone result. Some mutations are missing due to CN=0 by ascat. I tried to recovery those mutations

#

# two mutations targeting on TBL1XR1 in P46 have wrong AF by IGV inspection. Those mutations were also manually added

#######Aim of 

#######1. calcuate which mutation is clonal, subclonal

#######2. Check whether's possible to detect the variants in control DNA and which group of mutations can be detected in blood DNAs

####data.frame empty.raw: raw mutation file.

####this file: from raw VCF files. add additional columns to indicate the mutation is shared by all biopsy for a given patient.

####data.frame empty pyclone result

###flcombined imported from 16_combined/fl42_combined_updated010719.txt the latest mutation file, original file from Sigve

library(plyr)

library(dplyr)

library(ggplot2)

setwd("G:/FL_resequncing/FL_exome_final/fl_latest/11_pyclone/BED_bam/counts_new/CPC_codingPlusUTR_pyclone/All_codingPlusUTR/All_codingPlusUTR_i1")

flcombined<- read.csv(file="../fl_latest/16_combined/fl44_reseq_final_060619.txt", sep="\t",header=T, stringsAsFactors = F)

mutation<- unique(flcombined%>%select(Patient, Gene, SYMBOL,VAR_ID, CDS_CHANGE,HGVSp_short,Region,

                                      VARIANT,VARIANT_CLASS))

mutation$chr<- as.character(lapply(strsplit(as.character(mutation$VAR_ID), "\\_"),"[",1))

mutation$star<- as.character(lapply(strsplit(as.character(mutation$VAR_ID), "\\_"),"[",2))

mutation$mutation_id<- paste(mutation$chr, mutation$star,sep="_")

###process JP6 seperately because JP6 does not have raw VCF file

JP6_anno<- mutation%>%filter(Patient=="P6B")

#########annotate pyclone result

pyclone.anno<- merge(empty, mutation, by=c("Patient", "mutation_id"))

###raw vcf df

###mutation calling

empty.raw$mutation_id<- paste(empty.raw$V1, empty.raw$V2,sep="_")

raw.vcf<- empty.raw%>%select(Patient,SAMPLE,mutation_id, VAR_ID, V7,V8, num_biopsy,sharedbyAll )

raw.vcf.anno<- merge(raw.vcf, mutation, by=c("Patient", "VAR_ID"))

##extract only coding

raw.vcf.anno.coding<- raw.vcf.anno%>%filter((Region=="Coding")| (VARIANT=="3'UTR") | (VARIANT=="5'UTR"))

#########

fl32<- subset(raw.vcf.anno.coding, Patient %in% unique(pyclone.anno$ Patient))

names(fl32)[names(fl32)=="mutation_id.x"]<- "mutation_id"

pyclone.df<- pyclone.anno%>%select(SAMPLE, mutation_id, cellular_prevalence,variant_allele_frequency, mean, class)

fl32.pyclone<- merge(fl32, pyclone.df, by=c("SAMPLE", "mutation_id"), all.x = T)

fl32.pyclone.removed<- fl32.pyclone%>%filter(is.na(cellular_prevalence))

missed_gene<- c("TAS1R3","ATAD3C","MMP23B","MORN1","HES2","TRAV8-6","SNRPA1","DLL4","VPS39","PAQR5","CEMIP","AMELY","ANHX","ATAD3C","B3GALT6","BRCA2","CBLB","CEMIP","DLL4","ESPN","FAM46A","HES2","IL5RA","IL9R","LCE1E","LEPROTL1","MAP1A","MMP23B","MORN1","PAQR5",

                "PPIAL4G","RAB7A","RORC","SNRPA1","TAS1R3","TNFRSF14","TRAV8-6","TRMT10A","TSPY1","USF3","USP9Y","VPREB1","VPS39","ZBTB20","ZNF268")

fl32.missed<- subset(fl32.pyclone.removed, SYMBOL %in%missed_gene)

fl32.ig<- setdiff(fl32.pyclone.removed, fl32.missed)

#######have to manually add clonal informations for missed genes

write.table(fl32.missed,file="clonal_gene/fl32.missed_genes.txt",sep="\t",row.names=F, quote=F )

########read manually revised genes

fl32.missed<- read.table(file="clonal_gene/fl32.missed_genes_manulrevised.txt",sep="\t", quote="",header=T,stringsAsFactors = F)

fl32.pyclone.kept<- fl32.pyclone%>%filter(!is.na(cellular_prevalence))

fl32.pyclone.new<- rbind(fl32.pyclone.kept, fl32.missed)

########defining mutation order: early dominant, dominant, subclonal

###

biopsy_empty_new<- fl32.pyclone.new[-FALSE,]

for (i in unique(fl32.pyclone.new$SAMPLE)){

  biopsy<- subset(fl32.pyclone.new, SAMPLE==i)

  biopsy_specific<- biopsy%>%filter(variant_allele_frequency>0)

  biopsy_specific$clonal<- ifelse(biopsy_specific$cellular_prevalence>=0.5 | biopsy_specific$mean>=0.5, "dominant", "subclonal")

  biopsy_empty_new<- rbind(biopsy_empty_new, biopsy_specific)

}

removed<- fl32.pyclone.new%>%filter(variant_allele_frequency==0)

##remove"TBL1XR1", discard_rescue_TBL1XR1

recover<- subset(removed, class=="discard_rescue_TBL1XR1")

recover$clonal<- "dominant"

fl32.missed$clonal<- fl32.missed$class

######

fl32.all<- rbind( biopsy_empty_new, recover,fl32.missed)

mutation_af<- ddply(fl32.all, .(mutation_id, Patient), summarise, min_CF=min(cellular_prevalence))

fl32.all.af<- merge(fl32.all, mutation_af,by=c("Patient","mutation_id"))

fl32.all.af$clonal[fl32.all.af$min_CF>0.5 & fl32.all.af$sharedbyAll=="Y" ]<-"early_dominant"

###corrected on April 5th 2019

fl32.all.af$clonal[fl32.all.af$class=="CPC_dominant_lost_LOH" & fl32.all.af$clonal=="dominant"]<-"early_dominant"

fl32.all.af$clonal[fl32.all.af$class=="CPC_dominant"]<- "early_dominant"

fl32.final<- fl32.all.af[,-c(16,17,18)]

##########process JP6 ,see previous

JP6_final<- JP6_df[,-c(4,5,7,9,11,13,14,15,16,17,18,19)]

JP6_final<- ddply(JP6_final, .(mutation_id), mutate, min_CF=min(cellular_prevalence))

fl32.all.final<-rbind(fl32.final,JP6_final)

##########

write.table(fl32.all.final, file="clonal_gene/fl32.all.coding_040319.txt",sep="\t", row.names=F,quote=F)

write.table(fl32.all.final, file="clonal_gene/fl32.all.coding_040519.txt",sep="\t", row.names=F,quote=F)

####plotting

fl32<- read.table(file="clonal_gene/fl32.all.coding_040519.txt",sep="\t",header=T,quote="",stringsAsFactors = F)

fl32.dominant<- fl32%>%filter(clonal !="subclonal")

fl32.dominant.ns<- fl32.dominant%>%filter(Region=="Coding" & VARIANT!="Silent")%>%filter(class!="discard")

fl32.dominant.ns.df<- unique(fl32.dominant.ns%>%select(Patient, SYMBOL,clonal))

fl32.dominant.ns.df.count<- ddply(fl32.dominant.ns.df, .(SYMBOL),mutate, num_pat=length(unique(Patient)))

clinic<- read.xlsx(file="../../../../../../2_clinic_infor/FL_Baoyan_clinical_file_FINAL_020719 - survival.xlsx",sheetIndex = 1)

clinic$group<- ifelse(is.na(clinic$Date_transf), "nFL", "tFL")

clinic.reduced<- clinic%>%select(Patient_ID, group, POD24)

names(clinic.reduced)<- c("Patient", "group", "POD24")

fl32.dominant.clinic<- merge(fl32.dominant.ns.df.count,clinic.reduced,by="Patient", all.x=T )

#########

discard_genes<- c("TTN","MUC16","MUC12")

fl32.dominant.clinic.clean<- subset(fl32.dominant.clinic, !SYMBOL %in% discard_genes)

#######taking genes mutated in at least 4 cases

dominant_4cases<- fl32.dominant.clinic.clean%>%filter(num_pat>=4)

dominant_4cases$SYMBOL<- reorder(dominant_4cases$SYMBOL, dominant_4cases$num_pat)

dominant_4cases$SYMBOL_ordered <- with(dominant_4cases, reorder(SYMBOL, SYMBOL, function(x) -length(x)))

####sepearte nFL and tFL

dominant_4cases.nFL<- dominant_4cases%>%filter(group=="nFL")

tmp.nFL<-ddply(dominant_4cases.nFL, .(Patient,SYMBOL), mutate, status=length(Patient))

tmp.nFL$clonal_correct<-tmp.nFL$clonal

tmp.nFL$clonal_correct[tmp.nFL$status=="2"]<-"Both"

nFL_plot<- unique(tmp.nFL[, -c(3)])

nFL_plot$clonal_correct[nFL_plot$clonal_correct=="dominant"]<- "Late"

nFL_plot$clonal_correct[nFL_plot$clonal_correct=="early_dominant"]<- "CPC"

nFL_plot$clonal_correct<-factor(nFL_plot$clonal_correct, levels=c("Both","Late","CPC"))

nFL<-ggplot(nFL_plot, aes(x=SYMBOL,fill=clonal_correct,drop=FALSE))+geom_bar(width=0.5) +

  scale_fill_manual(values=c("dark blue","pink","purple" ))+theme_bw()+coord_flip()+

  theme(legend.position = "none",

        axis.text.y=element_blank(),

        axis.ticks.y = element_blank(),

        axis.text.x = element_text(size = 18,face="bold"),

        legend.text=element_text(size=8,face="bold"),

        legend.title = element_text(size=8, face="bold"),

        plot.margin = unit(c(0.2,0,0.2,0.2), "cm"),

        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),

        panel.border= element_rect(color="blue",size=2))+

  scale_x_discrete(drop=FALSE)+

  scale_y_continuous(breaks = c(0,5,10,15))+xlab("")+ylab("")

############separate tFL group

dominant_4cases.tFL<- dominant_4cases%>%filter(group=="tFL")

tmp.tFL<-ddply(dominant_4cases.tFL, .(Patient,SYMBOL), mutate, status=length(Patient))

tmp.tFL$clonal_correct<-tmp.tFL$clonal

tmp.tFL$clonal_correct[tmp.tFL$status=="2"]<-"Both"

tFL_plot<- unique(tmp.tFL[, -c(3)])

tFL_plot$clonal_correct[tFL_plot$clonal_correct=="dominant"]<- "Late"

tFL_plot$clonal_correct[tFL_plot$clonal_correct=="early_dominant"]<- "CPC"

tFL_plot$clonal_correct<-factor(tFL_plot$clonal_correct, levels=c("Both","Late","CPC"))

tFL<-ggplot(tFL_plot, aes(x=SYMBOL,fill=clonal_correct))+geom_bar(width=0.5) +

  scale_fill_manual(values=c("dark blue","pink","purple" ))+coord_flip()+theme_bw()+

  theme(legend.position = "none",

        axis.text.y=element_blank(),

        axis.ticks.y = element_blank(),

        axis.text.x = element_text(size = 18,face="bold"),

        legend.text=element_text(size=8,face="bold"),

        legend.title = element_text(size=8, face="bold"),

        plot.margin = unit(c(0.2,0.2,0.2,-0.1), "cm"),

        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),

        panel.border=element_rect(color="red", size=2))+xlab("")+ylab("")+

  scale_x_discrete(drop=FALSE)+scale_y_continuous(breaks = c(0,5,10,15))

pdf(file="CPC_plots/24genes_freq_groups_20191203.pdf",width=4, height=6)

plot_grid(nFL, 

          tFL, align="h"

          , nrow = 1,rel_widths=c(1,1))

dev.off()

#######try to plot them together, waterfall plot, frequency plot, clinic data

all<- rbind(nFL_plot, tFL_plot)

patient_order_nFL<- unique(fl32.dominant.clinic.clean$Patient[fl32.dominant.clinic.clean$group=="nFL"])

patient_order_tFL<- unique(fl32.dominant.clinic.clean$Patient[fl32.dominant.clinic.clean$group=="tFL"])

all$Patient<- factor(all$Patient, levels=c(patient_order_nFL,patient_order_tFL))

all$clonal_correct<- factor(all$clonal_correct,levels=c("Both", "CPC", "Late", "NA"))

all.tmp<- all%>%select(Patient, SYMBOL, clonal_correct)

all.tmp.cast<- cast(all.tmp, Patient~SYMBOL)

all.melt.again<- melt(all.tmp.cast, id.vars="Patient")

all.melt.again$value<- factor(all.melt.again$value, levels=c("Both", "CPC", "Late"))

all.plot<- ggplot(all.melt.again, aes(x=Patient, y=SYMBOL, fill=value))+

  geom_tile(color="gray", size=0.1)+ scale_fill_manual(values=c("dark blue","purple","pink"),drop=FALSE)+scale_x_discrete(drop=FALSE)+

  theme(panel.background = element_blank(),

        axis.text.x=element_blank(),

        axis.ticks.x=element_blank(),

        axis.text.y =element_text(size=15,face="bold"),

        legend.position = "NULL",

        plot.margin = unit(c(0.5,-0.5,-0.5,0.5), "cm"))+ xlab("")+ylab("")

clinic.32<-unique(fl32.dominant.clinic.clean%>%select(Patient, group, POD24))

names(clinic.32)<- c("Patient", "Group", "POD24")

clinic.32$POD24<- factor(clinic.32$POD24, levels=c(0,1))

clinic.32$Patient<- factor(clinic.32$Patient, levels=c(patient_order_nFL,patient_order_tFL))

clinic.32.melt<- melt(clinic.32, id.vars = "Patient")

clinic_plot<- ggplot(clinic.32.melt, aes(x=Patient, y=variable, fill=value))+

  geom_tile()+

  scale_fill_manual(values=c("blue", "red", "grey","black"))+

  theme(panel.background = element_blank(),

        axis.text.x = element_blank(),

        axis.text.y =element_text(size=16,face="bold"),

        axis.ticks.x=element_blank(),

        legend.position = "none",

        plot.margin = unit(c(-0.1,0.5,-0.5,0.5), "cm"))+

  xlab("")+ylab("")

library(egg)

####legend

tFL_plot$clonal_correct<- factor(tFL_plot$clonal_correct, levels=c("Both", "CPC", "Late"))

pdf(file="CPC_plots/Figure4B_20191112_legend.pdf",width=12,height=13)

ggplot(tFL_plot, aes(x=SYMBOL,fill=clonal_correct))+geom_bar(width=0.5) +

  scale_fill_manual(values=c("dark blue","purple","pink"))+coord_flip()+theme_bw()+

  theme(

    axis.text.y=element_blank(),

    axis.ticks.y = element_blank(),

    axis.text.x = element_text(size = 18,face="bold"),

    legend.text=element_text(size=14,face="bold"),

    legend.title = element_text(size=16, face="bold"),

    plot.margin = unit(c(0.2,0.2,0.2,-0.2), "cm"),

    panel.grid.major = element_blank(), panel.grid.minor = element_blank())+xlab("")+ylab("")+labs(fill = "Timing")+

  scale_x_discrete(drop=FALSE)+panel_border(color="red", size=2)+scale_y_continuous(breaks = c(0,5,10,15))

dev.off()

save.image("G:/FL_resequncing/FL_exome_final/fl_latest/11_pyclone/BED_bam/counts_new/CPC_codingPlusUTR_pyclone/All_codingPlusUTR/All_codingPlusUTR_i1/Figure4B_20191112.RData")

###updated on 2019.11.15

pdf(file="CPC_plots/24genes_groups_20191203.pdf",width=6, height=5.8)

plot_grid(all.plot, 

          clinic_plot, align="v"

          , nrow = 2,rel_heights =c(14,2))

dev.off()