---
editor_options:
  chunk_output_type: console
output:
  html_document:
    df_print: paged
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,fig.width = 7, fig.height = 7)
```

## Supplemental Data 1 for Boot et al., "Characterization of colibactin-associated mutational signature in an Asian oral squamous cell carcinoma and in other mucosal tumor types"

These are the scripts used to search for the presence of the novel mutational signature (SBS_A^n^T)
in the data set of 19,184 whole-exomes and 4,645 whole-genomes 
that was compiled for the PCAWG on Mutational Signatures (Alexandrov et al., 2018).

### Background 

SBS_A^n^T is characterized by
extremely strong enrichment for adenines 4bp and 3bp 5' of mutated thymines.

We will use this characteristic
to search for this mutational signature in the whole-exome and whole-genome data.

To test for enrichment we perform binomial tests, comparing: 
a) the proportion of T>N mutations with an adenine 3 or 4bp 5' against 
b) the proportion thymines in the human genome or exome that has an adenine 3 or 4bp 5'

The null hypothesis is that the proportion of T>N mutations with A at -3 (or -4) is 
equal to the proportion of thymines in the human genome/exome with A at -3 (or -4)

### Methods / code 

We retrieved the variant calls for 19,184 whole-exomes and 4,645 whole-genomes, and counted:

- T2Ncount:		Sum of the number of T>A, T>C and T>G single nucleotide substitutions
- ANNTcount:		Total number of T>A, T>C and T>G mutations where there is an adenine 3bp 5' of the mutated T
- ANNNTcount:		Total number of T>A, T>C and T>G mutations where there is an adenine 4bp 5' of the mutated T

We load these data into data frame *df*, which is keyed by the name of the sample and also contains the field *dataType* indicating whether the data are whole-genome or whole-exome.

```{r load data}
df<-read.csv("input/mutCounts_allExomes+allGenomes.txt",sep="\t",as.is=T)
```

```{r minimum T2N mutation count}
## keep only samples with >= 50 (exome) or >= 500 (genome) thymine mutations
df<-df[c(which(grepl("WGS",df$dataType) & df$T2Ncount>499),
         which(grepl("WES",df$dataType) & df$T2Ncount>49)), ]

## load the tri- tetra- and pentanucleotide abundances. After load we have
## variable *opp* which contains abundance of tri, tetra, and penta nucleotides
## for the human genome and exome.
load("input/trinucContent_hs37d5.RData")
```

```{r helper function}
source('src/mSigTools.v0.13.R')
rev.comp <- function(seq) { revc(toupper(seq)) }
```

The following function tests for statistically significant enrichment for the 
adenines at the -3 and -4 position relative to mutations from thymines. 

```{r function for statistical analysis}
enrichmentTest<-function(df){
  df$pval_ANNT2ANNN<-NA
  df$pval_ANNNT2ANNNN<-NA
  
  # Create columns for the reverse complements of tetra and penta nucleotides
  tmpOpp4<-opp$`4bp`
  for(ii in 1:nrow(tmpOpp4)){tmpOpp4$revComp[ii]<-rev.comp(tmpOpp4$NA.[ii])}
  tmpOpp5<-opp$`5bp`
  for(ii in 1:nrow(tmpOpp5)){tmpOpp5$revComp[ii]<-rev.comp(tmpOpp5$NA.[ii])} 
  
  for(i in 1:nrow(df)){
    dataType<-df$dataType[i]
    if(grepl("WGS",dataType)) dataType <- "genome"
    if(grepl("WES",dataType)) dataType <- "exome..agilent.V6."
    
    # test1: 
    #
    # - H0: proportion of T>N mutations with A at -3 is the same as the
    # proportion of all T's in the genome with A at -3. 
    # 
    # - HA: proportion of T>N mutations with A at -3 is > the above
    
    obs.mut<-df$T2Ncount[i]         ## total number of T>N mutations
    obs.mut.target<-df$ANNTcount[i] ## total number of T>N mutations with an A 3bp 5'
    
    abundance_T2N<-sum(tmpOpp4[substr(tmpOpp4$NA.,4,4) == "T" |
                                 substr(tmpOpp4$revComp,1,1) == "A",
                               dataType])
    
    abundance_ANNT2ANNN<-sum(tmpOpp4[(substr(tmpOpp4$NA.,4,4) == "T" &
                                        substr(tmpOpp4$NA.,1,1) == "A")
                                     |
                                       (substr(tmpOpp4$revComp,1,1) == "A" &
                                          substr(tmpOpp4$revComp,4,4) == "T"),
                                     dataType])
    
    stopifnot(mode(obs.mut.target) == "numeric")
    stopifnot(mode(obs.mut) == "numeric")
    testResult<-binom.test(obs.mut.target,
                           obs.mut,
                           abundance_ANNT2ANNN/abundance_T2N,
                           alternative = "greater")
    df$pval_ANNT2ANNN[i]<-testResult$p.value
    
    #  test2:
    #
    # - H0: proportion of T>N mutations with A at -4 is the same as the
    # proportion of all T's in the genome with A at -4.
    #
    # - HA: proportion of T>N mutations with A at -4 is > the above
    
    obs.mut<-df$T2Ncount[i]         ## total number of T>N mutations
    obs.mut.target<-df$ANNNTcount[i] ## total number of T>N mutations with an A 3bp 5'
    abundance_T2N<-sum(tmpOpp5[substr(tmpOpp5$NA.,5,5) == "T" |
                                 substr(tmpOpp5$revComp,1,1) == "A",
                               dataType])
    
    abundance_ANNNT2ANNNN<-sum(tmpOpp5[(substr(tmpOpp5$NA.,5,5) == "T" &
                                          substr(tmpOpp5$NA.,1,1) == "A" ) 
                                       |
                                         (substr(tmpOpp5$revComp,1,1) == "A" &
                                            substr(tmpOpp5$revComp,5,5) == "T")
                                       ,dataType])
    
    stopifnot(mode(obs.mut.target) == "numeric")
    stopifnot(mode(obs.mut) == "numeric")
    testResult<-binom.test(obs.mut.target,
                           obs.mut,
                           abundance_ANNNT2ANNNN/abundance_T2N,
                           alternative = "greater")
    df$pval_ANNNT2ANNNN[i]<-testResult$p.value
  }
  return(df)
}

```

Find samples with mutations from thymine that have statistically significant
enrichment for adenines at the -3 and -4 positions and with proportions of
adenines above determined cutoffs.

```{r perform analysis}
## add the counts for sample 62074759
df<-rbind(df,'62074759T'=list("WGS",34905,27375,22695,18549))
df<-enrichmentTest(df)

## only look at tumors that fit all criteria
result<-df[ p.adjust(df$pval_ANNT2ANNN,method="BH") < 0.05 &
              p.adjust(df$pval_ANNNT2ANNNN,method="BH") < 0.05 ,]

## plot scatter of enrichment for adenines 3 and 4bp 5' of mutated thymines
plot(df$ANNTcount/df$T2Ncount,df$ANNNTcount/df$T2Ncount,
     pch=16,cex=0.5,
     main="Identification of samples with enrichment of adenines\n(3 and 4bp 5' of mutated thymines)",
     xlab="Proportion thymines with adenine 3bp 5'",
     ylab="Proportion thymines with adenine 4bp 5'")
points(result$ANNTcount/result$T2Ncount,result$ANNNTcount/result$T2Ncount,
       pch=16,cex=0.5,col="magenta")
points(result$ANNTcount[nrow(result)]/result$T2Ncount[nrow(result)],
       result$ANNNTcount[nrow(result)]/result$T2Ncount[nrow(result)],
       pch=8,cex=0.8,col="red",lwd=1.5)
## expected ratio of A's 5' of thymines (abundance in the genome)
points(c(325868953,3890707)/c(1074123658,14758195),
       c(318941964,3724456)/c(1074123658,14758195),
       pch=c(3,4),cex=0.8,col="lightgreen",lwd=3)
legend("topleft",
       legend=c("not significant","significant","62074759",
                "expected_exomes","expected_genomes"),
       col=c("black","magenta","red","lightgreen","lightgreen"),
       pch=c(16,16,8,4,3),
       pt.lwd=c(1,1,1.5,2,2),
       pt.cex=1.1,cex=0.7,bty="n",ncol=1)

selected <- result[result$ANNTcount/result$T2Ncount > 0.4 | 
                   result$ANNNTcount/result$T2Ncount > 0.4,]

## format to print output table to html
selected$pval_ANNT2ANNN<-format(selected$pval_ANNT2ANNN,digits=3)
selected$pval_ANNNT2ANNNN<-format(selected$pval_ANNNT2ANNNN,digits=3)
library(knitr)
knitr::kable(selected,row.names = T,
             align = "c",
             caption = paste("Whole-exome and whole-genome sequenced samples",
                             "with likely AnT exposure based on enrichment of",
                             "adenines 5' of mutated thymines"))

write.table(selected,file="candidatesW_AnT.txt",sep="\t",quote=F,row.names = T)

```





