---
title: "Figures for Famdenovo Output"
output:
  pdf_document: default
  html_document:
    df_print: paged
    toc: true
    toc_float: true
toc: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r Libs and Options, warning=FALSE, message=FALSE, echo=FALSE}
# clean
rm(list = ls(all=TRUE))

# options
options(stringsAsFactors = F)

# libraries
library(knitr)
library(pROC)
library(ROCR)
library(tidyverse)
library(ggrepel)
library(beeswarm)
library(Hmisc)

# reading tables
# Famdenvo.TP53.cohort.mutation.data.tsv
cohort_tp53_table <- read.table('Famdenvo.TP53.cohort.mutation.data.tsv', header = TRUE, sep = "\t", stringsAsFactors = FALSE)
# Famdenovo.TP53.cohort.cancer.data.tsv
cancer_tp53_table <- read.table('Famdenovo.TP53.cohort.cancer.data.tsv', header = TRUE, sep = "\t", stringsAsFactors = FALSE)
```


```{r silent reading, message=FALSE, echo=FALSE, warning=FALSE}
# run fans function
FamFunc <- function(dat) {
  fams <- unique(dat$fam.id)
  extract <-  data.frame(fam.id = fams)
  for(i in fams) {
    extract[extract$fam.id == i, "truth"] <- ifelse(sum(dat[dat$fam.id == i, "state"] == "denovo"), 1, 0)
    extract[extract$fam.id == i, "prob"] <- max(dat[dat$fam.id == i, "prob"])
  }
  return(extract)
}

# MDA + NCI + DFCI + CHOP
#comb4 <- rbind(mda.family, nci.family, dfci.family, chop.family)
# re-organize data table
comb_df <- filter(cohort_tp53_table, cohort_tp53_table$truth!="unknown")
comb_df2 <- select(comb_df, c(id="Nu.ID", fam.id="Nu.Fam", "prob", state="truth"))
comb_df3 <- distinct(comb_df2)
comb4 <- FamFunc(comb_df3)

# comparison criteria
# read in 2 criteria
roc.relaxed <- read.csv("result_relaxed.csv", header=T)
roc.stringent <- read.csv("result_stringent.csv", header=T)

# get age points to annotate on plot
age_points_stringent = roc.stringent[roc.stringent$age.cut>30 & roc.stringent$age.cut<45,]
age_points_relaxed = roc.relaxed[roc.relaxed$age.cut>40 & roc.relaxed$age.cut<70,]
age_points_relaxed$tpr[age_points_relaxed$age.cut==60]=age_points_relaxed$tpr[age_points_relaxed$age.cut==60]+0.02
```


\newpage

## 1. ROC
### 1.1 Structure of the data sets

          The data sets in section 1 have the following structure.

```{r, warning=FALSE, message=FALSE}
#  (mda.family, nci.family, chop.family, dfci.family)
# structure for the data set using mda as an example
knitr::kable(data.frame(
  Variables = colnames(comb4),
  Meaning = c("Family ID", "Truth", "De Novo Probability")
), align = 'c')

# these are the first 3 rows
knitr::kable(head(comb4, n=3), align = 'c')
```


### 1.2 ROC curve

```{r Figure 3 ROC, warning=FALSE, message=FALSE}
# get roc object
rocobj = roc(comb4$truth, comb4$prob)
rocobj$fpr <- 1-rocobj$specificities
rocobj$tpr <- rocobj$sensitivities

# get ci
ci_rocobj = ci.thresholds(rocobj, boot.n=1000)
rocobj_ci_specificity <- as.vector(ci_rocobj$specificity[,2])
ci_rocobj_fpr <- 1- rocobj_ci_specificity
ci_rocobj_upper_bound <- as.vector(ci_rocobj$sensitivity[,3])
ci_rocobj_lower_bound <- as.vector(ci_rocobj$sensitivity[,1])

# plot the ROC
par(mar=c(4,4,1,1))
plot(y=rocobj$sensitivities, x=1-rocobj$specificities,ylab="True Positive Rate",
     xlab="False Positive Rate", xlim = c(0.03,0.99), ylim = c(0.03,0.99), type='l')
polygon(c(ci_rocobj_fpr , rev(ci_rocobj_fpr )), c(ci_rocobj_upper_bound, rev(ci_rocobj_lower_bound)), col = "cadetblue1", border = NA)
lines(y=rocobj$sensitivities, x=1-rocobj$specificities, lwd=3)
points(c(1 - 0.9514, 1 - 0.9306, 1 - 0.9167, 1 - 0.8750), c(0.7143, 0.7619, 0.7857, 0.8810), col = "red", pch = 19, cex = 1.2) # cutoff = 0.3, 0.2, 0.1, 0.05
points(roc.relaxed$fpr, roc.relaxed$tpr, col="purple", pch=23)
points(roc.stringent$tpr~roc.stringent$fpr, col = "darkorange", pch = 8, cex = 1.2)
legend('bottomright', legend=c('Famdenovo', 'Criterion 1', 'Criterion 2','95% CI'),
       pch=c(15, 23,8,15), pt.cex=c(1,1,1,2), box.lwd = 0,box.col = "white",bg = "white",
       col=c('black',  'purple', 'darkorange', 'cadetblue1'), inset=0.01)
abline(h = roc.relaxed$tpr[9], col = "gray", lty = 2) # age 45
abline(h = roc.relaxed$tpr[13], col = "gray", lty = 2) # age 65
text(y=age_points_stringent$tpr,x=age_points_stringent$fpr+0.035,labels=age_points_stringent$age.cut, col="darkorange", cex=0.9)
text(y=age_points_relaxed$tpr,x=age_points_relaxed$fpr-0.035,labels=age_points_relaxed$age.cut, col="purple", cex=0.9)
text(x=c(1 - 0.9514-0.03, 1 - 0.9306-0.04, 1 - 0.9167-0.04, 1 - 0.8750-0.04), y=c(0.7143, 0.7619, 0.7857+0.02, 0.8810+0.01),labels = c(0.3, 0.2, 0.1, 0.05) ,col = "red", cex=0.9) 
```


## 2. "de novo" probabilities in validation and discovery sets


```{r, message=FALSE, warning=FALSE, echo=FALSE}
# validation data
valid_data <- select(cohort_tp53_table, c('fam.id'='Nu.Fam', 'id'='Nu.ID', 'prob', 'institution',state='truth'))
valid_data <- filter(valid_data, valid_data$state != 'unknown')
# per cohort
mda.valid <- distinct(filter(valid_data, valid_data$institution=='MDA_pros'))
chop.valid <- distinct(filter(valid_data, valid_data$institution=='CHOP'))
nci.valid <- distinct(filter(valid_data, valid_data$institution=='NCI'))
dfci.valid <- distinct(filter(valid_data, valid_data$institution=='DFCI'))

# validation
extract.valid.func <- function(dat){
  fams <- unique(dat$fam.id)
  extract <- data.frame(fam.id=fams)
  for (i in fams) {
    extract[extract$fam.id == i, "prob"] <- max(dat[dat$fam.id == i, "prob"])
    extract[extract$fam.id == i, "truth"] <- ifelse(sum(dat[dat$fam.id == i, "state"] == "denovo"), T, F)
    }
  extract$pred <- extract$prob >= 0.2
  extract$ifcorr <- extract$pred == extract$truth
  return(extract)
}

mda.valid.extract <- extract.valid.func(mda.valid)
nci.valid.extract <- extract.valid.func(nci.valid)
dfci.valid.extract <- extract.valid.func(dfci.valid)
chop.valid.extract <- extract.valid.func(chop.valid)

# discovery data
disc_data <- select(cohort_tp53_table, c('fam.id'='Nu.Fam', 'id'='Nu.ID', 'prob', 'institution',state='truth'))
disc_data <- filter(disc_data, disc_data$state == 'unknown')
# per cohort
mda.disc <- distinct(filter(disc_data, disc_data$institution=='MDA_pros'))
chop.disc <- distinct(filter(disc_data, disc_data$institution=='CHOP'))
nci.disc <- distinct(filter(disc_data, disc_data$institution=='NCI'))
dfci.disc <- distinct(filter(disc_data, disc_data$institution=='DFCI'))

extract.disc.func <- function(dat){
  fams <- unique(dat$fam.id)
  fam.repr <- data.frame(fam.id=fams)
  for (i in fams) fam.repr[fam.repr$fam.id == i, "prob"] <- max(dat[dat$fam.id == i, "prob"])
  return(fam.repr)
}

mda.disc.extract <- extract.disc.func(mda.disc)
nci.disc.extract <- extract.disc.func(nci.disc)
dfci.disc.extract <- extract.disc.func(dfci.disc)
chop.disc.extract <- extract.disc.func(chop.disc)
```


### 2.1 Structure of the data sets
          The data sets in section 2 have the following structure.


```{r}
#  (mda.valid.extract, nci.valid.extract, chop.valid.extract, dfci.valid.extract)
# structure for the data set using mda.valid as an example
knitr::kable(data.frame(
  Variables = colnames(mda.valid.extract),
  Meaning = c("Family ID","De Novo Probability", "True State", "Prediction", "Is The Prediction Correct")
), align = 'c')

# these are the first 3 rows
knitr::kable(head(mda.valid.extract, n=3), align = 'c')

# (mda.disc.extract, nci.disc.extract, chop.disc.extract, dfci.disc.extract)
# structure for the data set using mda disc as an example
knitr::kable(data.frame(
  Variables = colnames(mda.disc.extract),
  Meaning = c("Family ID", "De Novo Probability")
), align = 'c')

# these are the first 3 rows
knitr::kable(head(mda.disc.extract, n=3), align = 'c')
```

```{r Figure 4, message=FALSE, warning=FALSE}
# ------ De Novo Probability ------ #
# validation set
valid <- rbind(data.frame(name="MDA (26/82)", prob=mda.valid.extract$prob, truth=mda.valid.extract$truth, ifcorr=mda.valid.extract$ifcorr), 
               data.frame(name="NCI (10/66)", prob=nci.valid.extract$prob, truth=nci.valid.extract$truth, ifcorr=nci.valid.extract$ifcorr), 
               data.frame(name="DFCI (4/30)", prob=dfci.valid.extract$prob, truth=dfci.valid.extract$truth, ifcorr=dfci.valid.extract$ifcorr), 
               data.frame(name="CHOP (2/8)", prob=chop.valid.extract$prob, truth=chop.valid.extract$truth, ifcorr=chop.valid.extract$ifcorr))

# plot the distribution on the validation set
ggplot(valid, aes(factor(name, levels=c("MDA (26/82)", "NCI (10/66)", "DFCI (4/30)", "CHOP (2/8)")), prob), cex.lab = 1) +
  geom_violin(color="gray") + geom_jitter(height=0, width=0.2, aes(color=truth), size = 3) + 
  geom_hline(yintercept=0.2, color="gray", linetype="dashed") +
  labs(title="TP53 de novo probability in the validation set", x="", y="TP53 de novo probability") +
  theme(panel.background=element_blank(),
        plot.title=element_text(size = 20, hjust=0.5),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 10)) + 
  guides(color=FALSE)

# discovery
disc <- rbind(data.frame(name="MDA (23/58)", prob=mda.disc.extract$prob), 
              data.frame(name="NCI (7/12)", prob=nci.disc.extract$prob), 
              data.frame(name="DFCI (7/61)", prob=dfci.disc.extract$prob), 
              data.frame(name="CHOP (3/7)", prob=chop.disc.extract$prob))
disc$pred <- disc$prob >= 0.2

# plot the distribution on the discovery set
ggplot(disc, aes(factor(name, levels=c("MDA (23/58)", "NCI (7/12)", "DFCI (7/61)", "CHOP (3/7)")), prob)) +
  geom_violin(color="gray") + 
  geom_point(position=position_jitter(width = 0.2),shape = 21, color=ifelse(disc$pred == TRUE, 'gray50', 'grey70'), size = 3, fill = ifelse(disc$pred == TRUE, 'gray50', 'white'), stroke = 1)+
  geom_hline(yintercept=0.2, color="gray", linetype="dashed") +
  labs(title="TP53 de novo probability in the discovery set", x="", y="TP53 de novo probability") +
  theme(panel.background=element_blank(),
        plot.title=element_text(size = 20, hjust=0.5),  
        axis.text = element_text(size = 10), 
        axis.title = element_text(size = 10))+ 
  guides(color=FALSE)
```


## 3. Parental and diagnosis age analysis
### 3.1 Structure of the data set in this section

```{r, message=FALSE, warning=FALSE, echo=FALSE}
# make same table with github data
cancer.cohort <- distinct(inner_join(x=cohort_tp53_table, y=cancer_tp53_table, by=c('Nu.ID','Nu.Fam')))
cancer.cohort <- filter(cancer.cohort, cancer.cohort$proband=="proband")# & !is.na(cancer.github$cancer.first) )

# make needed columns
cancer.cohort$state <- cancer.cohort$truth
cancer.cohort$state[cancer.cohort$truth=="unknown"] <- sapply(as.vector(as.numeric(cancer.cohort$prob[cancer.cohort$truth=="unknown"])),
                                                              function(x){
                                                               if(x>=0.2){return('denovo')} 
                                                                else{return('familial')}})
cancer.cohort$new_state <- sapply(as.vector(cancer.cohort$state), function(x){
  if(x=='denovo'){return('de novo')}else{return(x)}
})
# make groups
cancer.cohort$groups <- paste(cancer.cohort$gender, cancer.cohort$state, sep = " ")
cancer.cohort$maternal_age <- cancer.cohort$maternal.age
cancer.cohort$paternal_age <- cancer.cohort$paternal.age
cancer.cohort$LFS.Spectrum <- as.character(cancer.cohort$cancer.first)
cancer.cohort$LFS.Spectrum[ !is.na(cancer.cohort$cancer.first) ] <-sapply(as.vector(cancer.cohort$cancer.first[!is.na(cancer.cohort$cancer.first) ]), function(x){
  if(x==1){return('sts')}else if(x==2){return('ost')}else if(x==3){return('brain')}
  else if(x==4){return('breast')} else if(x==5){return('acc')} else if(x==6){return('leukemia')}
  else if(x==7){return('lung')} else if(x==8){return('choroid')} else if(x==20){'other.lfs'}
  else if (x==50) {return('non.lfs')}else{return(NA)}
})
probands3 <- cancer.cohort
```

        
        The data set used in this section considers only the probands of each family 
        and has the following structure.
        
```{r, message=FALSE, warning=FALSE}
# first 3 rows
knitr::kable(head(probands3, n=3), align = 'c')
```
### 3.2 Several plots regarding age of diagnosis and parental age at birth of the proband

```{r Figure 5, message=FALSE, warning=FALSE}
#probands 2 mother
beeswarm(as.numeric(probands3$maternal_age) ~ probands3$new_state,
         pwpch = ifelse(probands3$state == "denovo", 16, 1),
         method="swarm", ylim = c(10, 60),
         pwcol = ifelse(probands3$state == "denovo", 1, 2),
         xlab = "", ylab="") 
title(ylab="Maternal age (years)")
legend("topright", legend = c("de novo", "familial"), col = c(1, 2))

# probands 2 father
beeswarm(as.numeric(probands3$paternal_age) ~ probands3$new_state,
         pwpch = ifelse(probands3$state == "denovo", 16, 1),
         method="swarm", ylim = c(10, 60),
         pwcol = ifelse(probands3$state == "denovo", 1, 2),
         xlab = "", ylab="") 
title(ylab="Paternal age (years)")
legend("topright", legend = c("de novo", "familial"), col = c(1, 2))

# all probands
par(mar=c(4,4,1,1))
beeswarm(as.numeric(probands3$diag.age.first) ~ probands3$LFS.Spectrum,#probands3$cancer.first,
         pwpch = ifelse(probands3$state == "denovo", 16, 1), 
         pwcol = ifelse(probands3$gender == 0, 3, 4), method="swarm",
         ylim = c(0, 80), xlab="", ylab="", las=2)
title(ylab="Age at diagnosis (years)")
legend("topright", legend = c("de novo, females", "de novo, males", "familial, females", "familial, males"), inset = c(0.25, 0) )
```

              
## 4. Supplementary: Distribution of family size in LFS and BRCA cohorts


```{r Supp Figure 1, warning=FALSE, message=FALSE, echo=FALSE}
# CGN-BRCA data
#BRCA_CGN <- read.table(file = "BRCA_CGN_Families_Size.txt", header = T, sep = '\t')

# how to plot starting from the Github tables
LFS_Families <- distinct(select(cohort_tp53_table, c('Nu.Fam', 'Fam.Size', 'truth')))
LFS_Families <- distinct(filter(LFS_Families, LFS_Families$truth != 'unknown'))
LFS_Families  <- distinct(select(LFS_Families, c('Nu.Fam', 'Fam.Size')))
LFS_Family_Size <- data.frame(table(LFS_Families$Fam.Size))
colnames(LFS_Family_Size) <- c('lfs_pedigree_size', 'lfs_families')
```

### 5.1 Supplementary: Structure of the data sets in this section

          These data sets contain the size of each family belonging to the LFS cohorts used for our
          Famdenovo.TP53 project and the BRCA families from the CGN cohort used for Famdenovo.BRCA. 
          Lets take a look at the structure of these data sets.

```{r, message=FALSE, warning=FALSE}
# LFS data set structure
knitr::kable(data.frame(Variable = colnames(LFS_Family_Size),
                        Meaning = c("Number of Members Available in Pedigree", "Number of Pedigrees with This Amount of Available Members")), align = 'c')
# first 3 rows
knitr::kable(head(LFS_Family_Size, n =3),align = 'c' )
```


### 5.2 Supplementary: Plots of Family Size Distributions

```{r, warning=FALSE, message=FALSE}
# --------- Distribution of Family Size --------- #
# LFS Data plot
plot(x=as.vector(LFS_Family_Size$lfs_pedigree_size), y=as.vector(LFS_Family_Size$lfs_families), type = 'h', col = "black",
      main="LFS Data Family Size Distribution", ylab = 'No. of Families',
     xlim = c(0,150), ylim = c(0,10), xlab = "Family Size")
  minor.tick(nx = 5, ny = 0)
lines(x=as.vector(LFS_Family_Size$lfs_pedigree_size), y=as.vector(LFS_Family_Size$lfs_families),
      type = 'p', pch = 19, cex = 2, col = "Blue")
```





