---
title: "Genome-wide discovery of somatic coding and non-coding mutations in pediatric endemic and sporadic Burkitt lymphoma"
author: "Bruno M. Grande^☯^, Daniela S. Gerhard^☯^, Nicholas B. Griner, Aixiang Jiang, Eric Y. Zhao, Luka Culibrk, Nicole Knoetze, Jeremy S. Abramson, Thomas Alexander, Hilary Allen, Jeffrey Bethony, Jay Bowen, Corey Casper, John Kim Choi, Tanja M. Davidsen, Maureen Dyer, Julie M. Gastier-Foster, Patee Gesuwan, Timothy C. Greiner, Thomas G. Gross, Nancy Lee Harris, Yiwen He, John D. Irvin, Elaine S. Jaffe, Steven J. Jones, Jean Paul Martin, Marie-Reine Martin, Sam M. Mbulaiteye, Charles G. Mullighan, Andrew J. Mungall, Constance Namirembe, Karen Novik, Ariela Noy, Abraham Omoding, Jackson Orem, John T. Sandlund, Roland Schmitz, Wyndham H. Wilson, Yussanne Ma, Marco A. Marra, Ryan D. Morin^☯^, Louis M. Staudt^☯^"
csl: "blood.csl"
bibliography: "bibliography.bib"
output:
  bookdown::word_document2:
    reference_docx: "reference.docx"
    toc: no
  bookdown::gitbook:
    dev: "svg"
    css: "manuscript.css"
    number_sections: no
    split_bib: yes
    split_by: "none"
---

☯ These authors have contributed equally.

```{r setup, echo=FALSE}
blgsp_criterion <- rprojroot::has_file("blgsp.Rproj")
project <- blgsp_criterion$find_file

knitr::opts_chunk$set(
  warning = FALSE, message = FALSE, echo = FALSE, comment = "", dev = "png",
  dpi = 300, fig.width = 7, fig.height = 5, fig.path='figures/')
```

```{r load}
source(project("src/R/load.R"))
```


# Abstract

Burkitt lymphoma (BL) is an aggressive B-cell non-Hodgkin lymphoma that may be
curable for sporadic pediatric cases in high-resource countries, but it remains
a deadly disease for endemic cases in equatorial Africa and older sporadic
cases. 
A number of genomic studies have focused on BL, but many questions remained out
of reach due to limited samples sizes and/or lack of multi-dimensional
molecular data, such as the elusive role of the Epstein–Barr virus (EBV).
We have assembled an unprecedented large and diverse BL cohort consisting of 
`r num$ebl` endemic and `r num$sbl` sporadic pediatric cases, and performed
whole genome, RNA, and microRNA sequencing on the tumor samples. 
We also re-analyzed genomic data from `r num$icgc` sporadic pediatric BL cases
included in the ICGC MMML-Seq project.
Through an integrative analysis of this dataset, we describe the
genome-wide landscape of somatic coding and non-coding mutations in BL.
We uncover genes as being significantly mutated in BL for the first time,
namely _SIN3A_, _HIST1H1E_, _CHD8_, _USP7_, and _RFX7_. 
We also identify sources of molecular heterogeneity, including differences in 
AID-mediated somatic hypermutation (SHM), germinal center cell-of-origin, and
genome-wide and gene-wise mutation rates.
Notably, we report `r smgs_assoc[variable == "ebv_status", sum(qval < 0.1)]`
genes or gene modules that are significantly differentially mutated based on 
EBV infection status, none of which are differentially mutated according to
clinical variant status or EBV genome type.  
The data presented here not only advances our understanding of BL biology 
through the identification of novel recurrently mutated genes, but it
demonstrates that EBV infection status is a more meaningful classification
criterion, resulting in more genetically distinct subtypes.


# Introduction

Burkitt lymphoma (BL) is an aggressive B-cell non-Hodgkin lymphoma. 
In high-resource countries, it is considered rare and curable in children, but
it remains a deadly disease as the most common pediatric cancer in equatorial
Africa [@Magrath2012-xr]. 
Survival is poor among African children with BL due to late stage at 
presentation and limitations in the ability to support intensive 
chemotherapeutic regimens necessary to induce durable remissions 
[@Buckle2016-hd]. 
Genetically, it features a hallmark chromosomal rearrangement that juxtaposes 
the _MYC_ proto-oncogene in proximity to an immunoglobulin heavy (_IGH_)or 
light chain (_IGK_/_IGL_) enhancer [@Dalla-Favera1982-lp; @Taub1982-fa].
The uneven geographic distribution of BL in Africa has prompted investigation
into etiological agents such as viruses and parasites 
[@Burkitt1962-ts; @Magrath2012-xr].
Cases arising in malaria-endemic regions, where the disease is most common, 
constitute the clinical variant known as endemic BL (eBL). 
In addition to an association with malaria, virtually every eBL case is 
infected by the Epstein–Barr virus (EBV), a ubiquitous virus that infects 90%
of the global adult population [@Cohen2000-kc].
Given EBV's potent ability of transforming B cells _in vitro_ [@Miller1982-jx],
researchers have long speculated on its oncogenic role in BL formation, but 
many questions remain unanswered.
Notably, little is known about what distinguishes BL-associated EBV from 
strains that infect healthy individuals who do not develop lymphoma.
Some have suggested a hit-and-run mechanism whereby EBV would infect B cells,
induce durable epigenetic changes, and then switch to a latent state
[@Gromminger2012-xa].
However, this hypothesis must be reconciled with the fact that EBV persists in
BL cells despite spontaneous episome loss in dividing cells [@Nanbo2007-fh]. 
One can only presume that EBV still provides some selective advantage to 
infected B cells, which outcompete cells that lose the virus.
EBV must achieve this effect with little to no protein expression—with the
exception of _EBNA1_—to avoid detection by the immune system.
Addressing these questions must start by improving our understanding of the 
biological differences associated with or induced by EBV in BL.

The second clinical variant, sporadic BL (sBL), has a much lower incidence and
is restricted to regions outside of malaria-endemic areas, _e.g._ North America
and Europe [@Sant2010-mu; @Morton2006-zh]. 
EBV infection rates are markedly lower among sBL cases at around 10-20% 
[@Schmitz2014-tz].
Given the extreme rate of cellular proliferation in BL, tumors tend to 
dramatically respond to intensive albeit toxic chemotherapy [@Costa2013-mv].
In the presence of proper supportive care, cure rates in children with sBL are
often above 90% [@Costa2013-mv].
However, treating adult and elderly patients remains a significant challenge, 
with overall survival rates being as low as 30% to 60% [@Costa2013-mv].
This major difference in response to therapy is often attributed to the reduced
ability of older patients of coping with intensive treatment regimens 
[@Jacobson2014-am].
Poor survival rates also apply to the third clinical variant, 
immunodeficiency-related BL, which is mostly caused by HIV infection 
[@Orem2009-dx; @Kojima2014-sg; @Puga_L2015-zi]. 
Thus, the need for an improved molecular understanding of this disease is true
even for countries with the resources to administer standard therapy and 
manage its toxicity. 
A number of genomic studies have advanced our understanding of BL, but many
questions remained out of reach due to limited samples sizes and/or lack of
multi-dimensional molecular data [@Schmitz2012-df; @Richter2012-ay;
@Love2012-qk; @Abate2015-iu; @Kaymaz2017-xm].
The Burkitt Lymphoma Genome Sequencing Project (BLGSP) is well-positioned to 
conduct an integrative molecular characterization of a large comprehensive BL
cohort including an unprecedented representation of endemic cases. 
In this first landmark BLGSP study, we focus on the genome-wide discovery of 
somatic coding and non-coding mutations in `r num$total` pediatric BL cases. 
We hypothesized on the existence of biological differences between BL subtypes
and suspected that EBV infection status would result in more significant
contrast than clinical variant status.
Through this work, we have identified sources of molecular heterogeneity in 
BL, including differences in AID-mediated somatic hypermutation (SHM),
germinal center cell-of-origin, and genome-wide and gene-wise mutation rates.
Many of these differences were either greater or exclusively significant 
according to EBV infection status rather than clinical variant status.
The data presented here demonstrates that EBV infection status is a more
meaningful classication criterion than the current clinical variant system,
resulting in more genetically distinct subtypes.
<!---
We also identify novel coding and non-coding mutation patterns in both 
established and previously unreported BL genes.
--->


# Results

## Clinical and molecular characteristics of the BLGSP pediatric discovery cohort

<!---
TODO
- Update sequencing statistics. 
- Carefully check the EBV type, esp. for the case with the deletion.
--->

The BLGSP pediatric discovery cohort consists of `r num$ebl` endemic BL (eBL)
cases from Uganda and `r num$sbl` sporadic BL (sBL) cases from the United 
States, amounting to `r num$blgsp` patients in total.
Endemic and sporadic EBV infection rates are consistent with previous reports
(`r num$ebl_ebv`% and `r num$sbl_ebv`%, respectively)
Conversely, our cohort includes `r num$ebl_ebvneg` EBV-negative eBL cases and
`r num$sbl_ebvpos` EBV-positive sporadic cases, which will assist in 
attributing biological variation to either clinical variant status or EBV 
infection status.
Tumor samples underwent deep whole genome sequencing (WGS; mean depth 82X;
range 55-96), ribo-depleted RNA sequencing (RNA-seq; mean 199M reads; range
132-255M), and microRNA sequencing (miRNA-seq; mean 13M reads; range 2.7-32M).
Constitutive DNA was sequenced (mean depth 41X; range 30-51) for all patients
to distinguish between somatic and germline mutations.
Being the putative cells-of-origin for BL, centroblasts and centrocytes from
six pediatric tonsil donors underwent the same RNA-seq and miRNA-seq protocols
to act as normal comparators for gene expression analyses.
We supplemented certain genomic analyses with re-analyzed WGS data from
`r num$icgc` pediatric sBL cases included in the International Cancer Genome 
Consortium (ICGC) Molecular Mechanisms in Malignant Lymphoma by Sequencing 
(MMML-Seq) project [@Richter2012-ay; @Hezaveh2016-pw].
Clinical and molecular characteristics are summarized in Table
\@ref(tab:table-univariate-stats).


(ref:caption-table-univariate-stats) – Univariate table listing the clinical and molecular characteristics for the pediatric cohorts of the BLGSP discovery and ICGC MMML-Seq projects.

```{r table-univariate-stats, results="asis"}
metadata$patient %>%
  left_join(metadata$genome[tissue_status == "Tumor", .(patient, ff_or_ffpe)],
            by = "patient") %>%
  left_join(myc_summary, by = "patient") %>%
  filter(!grepl("^TS", patient)) %>%
  mutate(
    clinical_variant = fct_relevel(clinical_variant, "Sporadic", "Endemic"),
    ebv_type = recode(ebv_type, None = "Negative"),
    age_group = case_when(
      age_at_initial_pathologic_diagnosis <= 5 ~ "0 - 5",
      age_at_initial_pathologic_diagnosis <= 10 ~ "6 - 10",
      age_at_initial_pathologic_diagnosis <= 15 ~ "11 - 15",
      age_at_initial_pathologic_diagnosis <= 20 ~ "16 - 20"),
    age_group = fct_relevel(age_group, "0 - 5", "6 - 10", "11 - 15"),
    sv = fct_infreq(sv),
    project = recode(project,
                     BLGSP = "BLGSP Discovery",
                     ICGC = "ICGC MMML-Seq")) %>%
  univariateTable(project ~ clinical_variant + ebv_status + ebv_type +
                    age_group + ff_or_ffpe + sv, .,
                  freq.format = "count(x) (percent(x)%%)", digits = 0,
                  clinical_variant = "Clinical variant", ebv_type = "EBV type",
                  ebv_status = "EBV status", age_group = "Age group (years)",
                  ff_or_ffpe = "Tumor biopsy", sv = "MYC translocation") %>%
  summary(show.pvalues = FALSE) %>%
  kable(caption = "(ref:caption-table-univariate-stats)")
```


## AID-mediated patterns of non-coding mutations

<!---
TODO
- Update wavelet results once the PAX5 enhancer peak is back.
- Provide percentage of non-coding mutations in clusters.
--->

Unlike previous genomic studies focused on BL, our dataset is sufficiently
large to characterize the genome-wide landscape of non-coding mutations.
We leveraged the Rainstorm method [@Arthur2017-kj] to discern regions enriched
in silent mutations.
This method located `r num$wavelets` clusters of non-coding mutations with an
median size of `r num$wavelets_size` bp (range `r num$wavelets_size_min`-
`r num$wavelets_size_max` bp).
In total, `r num$wavelets_olap` clusters (`r num$wavelets_olap_percent`%)
overlap at least one gene, and the median distance from the nearest 
transcription start site (TSS) is `r num$wavelets_dist_tss` bp.
If we only consider clusters that are mutated in `r wavelet_threshold`
or more cases, `r num$wavelets_near_tss` (`r num$wavelets_near_tss_percent`%)
are within `r wavelet_dist` bp of a TSS (Supplementary Figure
\@ref(fig:histogram-distfromtss)).
The most recurrently mutated clusters are altered in up to
`r num$wavelets_num_mutated_max` cases 
(`r num$wavelets_num_mutated_max_percent`%) and overlap established targets of
physiologic (_i.e._ _IGH_, _IGK_, _IGL_, and _MYC_) and aberrant 
(_e.g._ _BCL7A_, _BCL6_, _BACH2_, _TCL1A_, and _BTG2_ [@Jiang2012-ah]) somatic
hypermutation (SHM; Figure \@ref(fig:heatmap-ismutated-topwavelets)).
Here, we consider _MYC_ a target of physiologic SHM given its proximity to an
immunoglobulin (Ig) locus after the canonical IG-_MYC_ translocation 
[@Raffeld1995-ra; @Rabbitts1983-qg].
Altogether, a majority of these clusters can be attributed to AID-mediated SHM
given their targets and proximity to the nearest TSS [@Pasqualucci2001-ap].
Physiologic SHM is observed in nearly all cases, especially at the _IGH_ and
_MYC_ loci.
On the other hand, signs of aberrant SHM appear to be more prevalent in endemic 
or EBV-positive cases (Supplementary Figure
\@ref(fig:boxplot-covariate-ashmcount)).
This is consistent with the observation that AID transcript levels are higher 
in endemic or EBV-positive tumors (Figure 
\@ref(fig:boxplot-covariates-aidexpr)) according to the BLGSP RNA-seq data. 
The fold change in AID expression is slightly higher according to EBV infection 
status (`r num$aid_fc_ebv`) than clinical variant status (`r num$aid_fc_cv`),
hinting that EBV may be modulating AID expression.
We confirmed that AID expression is indeed higher in EBV-positive tumors
regardless of geographic origin (Supplementary Figure
\@ref(fig:boxplot-ebvtype-aidexpr-cv)).


(ref:caption-heatmap-ismutated-topwavelets) – Heat map showing the mutation rate of the most recurrently altered clusters of non-coding variants (altered in at least 15 cases). Clusters are labeled according to the nearest TSS (rows). Clusters that are associated with the same TSS are combined for plotting. Light gray denotes no mutations, and darker shades of gray indicate the presence of up to two mutations or more than two mutations, respectively.

```{r heatmap-ismutated-topwavelets, fig.height=5, fig.cap="(ref:caption-heatmap-ismutated-topwavelets)"}
wavelets_mat %>%
  {.[. > 2] <- 2; .}  %>%
  .[rowSums(. > 0) >= 15, ] %>%
  .[order(!rownames(.) %in% c("MYC", "IGH", "IGK", "IGL"), 
          -rowSums(. > 0)),] %>%
  .[, order(annotations$genome[colnames(.), "ebv_type"],
            rnorm(ncol(.)))] %>%
  pheatmap(annotation_col = annotations$genome, annotation_colors = colours,
           scale = "none", show_rownames = TRUE, cluster_rows = FALSE,
           gaps_row = 4, gaps_col = sum(annotations$genome$ebv_type == "None"),
           color = c(`0` = "grey80", `1` = "grey50", `2` = "grey30"),
           cluster_cols = FALSE, legend = FALSE)
```


(ref:caption-boxplot-covariates-aidexpr) – Box plots showing the variance-stabilized expression of AID in germinal center cells and tumor samples (classified according to either clinical variant status or EBV infection status).

```{r boxplot-covariates-aidexpr, fig.cap="(ref:caption-boxplot-covariates-aidexpr)"}
aid_expr %>%
  ggplot(aes(value, expr, fill = value)) +
  geom_boxplot() +
  facet_grid(~ variable, scales = "free_x") +
  scale_fill_manual(values = flatten_chr(colours), breaks = NULL) +
  scale_x_discrete(labels = c(None = "Negative",
                             Centroblasts = "Centro-\nblasts",
                             Centrocytes = "Centro-\ncytes")) +
  labs(x = NULL, y = "AID expression (variance-stabilized)")
```


While several clusters of non-coding mutations overlap known targets of
aberrant SHM, many affect genes that have not yet been linked to BL.
One cluster overlaps a validated enhancer that is known to regulate the
expression of _PAX5_, which plays an important role in B-cell differentiation.
Mutations in this enhancer were first described in chronic lymphocytic leukemia
(CLL) [@Puente2015-gx], although we observe a higher mutation rate in BL
(`r num$pax5_percent`% versus 11%).
Unlike in CLL, the mutation status of this enhancer is not significantly 
associated with _PAX5_ expression (Supplementary Figure
\@ref(fig:boxplot-pax5waveletstatus-expr-nearbygenes-wilcox)).
Although, the trend toward higher _PAX5_ expression when the enhancer is 
mutated leads us to suspect AID-mediated SHM of this presumed active enhancer.
<!---
As in CLL, the mutation status of this enhancer is significantly associated
with _PAX5_ expression (Supplementary Figure
\@ref(fig:boxplot-pax5waveletstatus-expr-nearbygenes-wilcox)).
Tumors with a mutated enhancer have higher _PAX5_ expression, leading us to 
suspect that this AID-mediated SHM of this presumably active enhancer.
--->
Another target of non-coding mutations is the miR-142 locus.
Variants affecting the mature sequence of miR-142 have been reported in DLBCL
and FL in up to 20% of cases [@Kwanhian2012-qb; @Hezaveh2016-pw], but not in BL
despite being investigated in one of those studies.
We observe `r num$mir142_3p_muts` cases (`r num$mir142_3p_muts_percent`%)
with mutations in the mature sequence of miR-142-3p.
Most mutations in this cluster are actually located one kilobase upstream of 
miR-142.
In total, these mutations affect `r num$mir142_percent`% of our cases.
Figure \@ref(fig:heatmap-ismutated-topwavelets) shows that they
are exclusive to EBV-positive cases (P-value `r num$mir142_ebvpos_pval`,
Fisher's exact test).
It is possible that these variants disrupt a regulatory element that modulates
the expression of the miR-142 or other nearby genes.
While we did not observe a difference in expression for miR-142-5p
(P-value `r num$mir142_5p_pval`) or miR-142-3p (P-value `r num$mir142_3p_pval`)
based on the mutation status of this cluster, a number of nearby genes were
differentially expressed (Supplementary Figure
\@ref(fig:boxplot-mir142waveletstatus-expr-nearbygenes-wilcox)).


(ref:caption-custom-noncoding-pax5) – Non-coding mutations in the _PAX5_ enhancer. The cluster in question is indicated by a black outline on the right-hand side. _PAX5_ is located on the left-hand side. Note the overlap between the cluster of non-coding mutations and the peak in histone 3 lysine 27 (H3K27) acetylation.

```{r custom-noncoding-pax5, fig.cap="(ref:caption-custom-noncoding-pax5)"}
include_graphics("images/noncoding_pax5.png")
```


(ref:caption-custom-noncoding-mir142) – Non-coding mutations in proximity of _miR-142_.

```{r custom-noncoding-mir142, fig.cap="(ref:caption-custom-noncoding-mir142)"}
include_graphics("images/noncoding_mir142.png")
```


Lastly, a cluster of non-coding mutations seem to target the promoter region of
the _PVT1_ non-coding locus, which encodes the _PVT1_ long non-coding RNA
(lncRNA) and multiple miRNAs including miR-1204.
_PVT1_ is 60 kilobases downstream from the _MYC_ locus and often harbours 
IG-_MYC_ translocations due to this proximity.
It is also a target of Myc-mediated upregulation [@Carramusa2007-oc], which may
account for its high expression in BL tumors.
The median expression values in transcripts per million (TPM) for _MYC_ and 
_PVT1_ are in the 99th percentile.
Paradoxically, this locus has also been identified as a p53-inducible target
gene [@Barsotti2012-nl].
The cluster in question is mutated in `r num$pvt1_percent`% of cases, and its
mutation status is significantly associated with higher _PVT1_ expression
(P-value `r num$pvt1_pval`).
Mutations in this cluster are also significantly associated with EBV-positive 
tumors (P-value `r num$pvt1_ebvpos_pval`, Fisher's exact test; Figure
\@ref(fig:heatmap-ismutated-topwavelets)).
Given p53's role in inducing the expression of this non-coding locus, we sought
to identify any enrichment of mutations in p53 binding sites genome-wide,
defined by chromatin immunoprecipitation paired-end tag sequencing (ChIP-PET)
peaks [@Wei2006-qw].
One such peak overlapped the _PVT1_ promoter region and cluster of non-coding
mutations, and was deemed significantly mutated by these non-coding variants
(Benjamini-Hochberg (BH)–adjusted P-value 0.0508, OncodriveFML).
While the driver status for these mutations needs to be functionally validated,
unpublished work has shown that _PVT1_ acts as a tumor suppressor gene in BL by
lowering _MYC_ expression. 
These mutations highlight a potential role for AID in promoting BL formation
through non-coding mutagenesis of biologically relevant targets.
<!---
despite published reports supporting an oncogenic 
role in other cancer types [@Guan2007-sk; @Tseng2014-ol].
--->


(ref:caption-custom-noncoding-pvt1) – Non-coding mutations in the _PVT1_ promoter region. Note the p53 chromatin immunoprecipitation (ChIP) peak that overlaps the mutations in this cluster.

```{r custom-noncoding-pvt1, fig.cap="(ref:caption-custom-noncoding-pvt1)"}
include_graphics("images/noncoding_pvt1.png")
```


## Differential AID expression is not related to cell-of-origin

Figure \@ref(fig:boxplot-covariates-aidexpr) shows that the difference in AID
expression between EBV-positive and -negative tumors approaches the difference
between germinal center centroblasts and centrocytes.
This led us to hypothesize that these differences in AID expression may be
caused by differences in the tumors' germinal center cell-of-origin.
To test this hypothesis, a supervised approach was devised to compare the
expression profile of each tumor with that of centroblasts and centrocytes
using a weighted voting algorithm [@Golub1999-qz; @Victora2012-pa].
The algorithm used 118 genes that are differentially expressed
between dark zone and light zone cells.
<!--- `r nrow(coo$genes)` --->
This resulted in a centroblast score for each tumor scaled between –1 and 1,
where positive and negative scores indicate a greater similarity to 
centroblasts and centrocytes, respectively.
The majority of tumors (`r num$coo_cb_percent`%) are more similar to cells
from the germinal center dark zone in terms of gene expression.
This is consistent with previous reports that BL tumors derive from
centroblasts [@Dave2006-pq; @Caron2009-nf].
However, we note the existence of a minority of BL tumors 
(`r num$coo_cc_percent`%) that are more similar to centrocytes.
Heterogeneity in the cell-of-origin has been previously described; although in
that study, most BL tumors were more similar to centrocytes [@Victora2012-pa].
The difference in centroblast score between EBV-positive and -negative cases is
small but significant (P-value `r num$coo_pval_ebv`, Wilcoxon test).
The same is true between clinical variants (P-value `r num$coo_pval_cv`).
However, the difference in centroblast score does not have the same magnitude
as the difference between centroblasts and centrocytes.
Therefore, it seems unlikely that cell-of-origin alone can account for the
significantly greater AID expression in EBV-positive tumors.
Of note, recent work has shown that EBV is capable of inducing AID expression 
in lymphoblastic cell lines (LCLs) [@Kalchschmidt2016-je].
EBV achieves this effect through durable epigenetic changes of the AID gene 
locus, which is mediated by the protein encoded by the latent gene _EBNA3C_.
Our findings seem to extend this observation to BL tumors _in vivo_ and
potentially implicate EBV in promoting AID-mediated mutagenesis in 
EBV-positive tumors.
<!---
One limitation of this approach is that it does not distinguish between B cells
undergoing differentiation in the germinal center for the first time and
B cells re-transiting through the germinal center for additional affinity
maturation.
--->


(ref:caption-boxplot-covariates-centroblastscore) – Box plots showing the centroblast score for germinal center cells and tumor samples (classified according to either clinical variant status or EBV infection status). Positive and negative scores correspond to greater similarity to centroblasts and centrocytes, respectively.

```{r boxplot-covariates-centroblastscore, fig.cap="(ref:caption-boxplot-covariates-centroblastscore)"}
coo$ps_df %>%
  mutate(
    signif = case_when(
      variable == "Clinical variant" & num$coo_pval_cv<0.05 ~ "P-value < 0.05",
      variable == "EBV type" & num$coo_pval_ebv<0.05 ~ "P-value < 0.05",
      variable == "Germinal center" & num$coo_pval_gc<0.05 ~ "P-value < 0.05",
      TRUE ~ "Not significant"),
    variable = fct_relevel(variable, "Germinal center")) %>% 
  ggplot(aes(value, ps)) +
  geom_boxplot(aes(alpha = signif, fill = value), outlier.alpha = 0) +
  geom_quasirandom(aes(alpha = signif), size = 1, color = "grey20") +
  geom_hline(yintercept = 0, linetype = 2, color = "grey40") +
  facet_grid(~ variable, scales = "free_x") +
  scale_x_discrete(labels = c(None = "Negative",
                             Centroblasts = "Centro-\nblasts",
                             Centrocytes = "Centro-\ncytes")) +
  scale_fill_manual(values = flatten_chr(colours), breaks = NULL) +
  scale_alpha_manual(limits = c("P-value < 0.05", "Not significant"), 
    values = c(`P-value < 0.05` = 1, `Not significant` = 0.3),
    guide = guide_legend(title.position = "top")) +
  place_legend("bottom") +
  coord_cartesian(ylim = c(-1, 1)) +
  labs(x = NULL, y = "Centroblast score", alpha = "Wilcoxon test")
```


## Structural deregulation of _MYC_ expression

<!---
TODO 
- I want to take a close look at non-cacnonical translocations near the IG 
  loci as well as MYC, BCL2 and BCL6. 
- I also want to create a summary Circos figure displaying the translocations
  between the MYC and IG loci where links are coloured by EBV infection status.
--->

Simple IG-_MYC_ translocations were identified in all but two cases, with a
majority involving _IGH_ and the remainder associated with _IGL_ and _IGK_
(Table \@ref(tab:table-univariate-stats)).
In one case (BL123), we detected two t(3;8)(q27;q24) rearrangements near
_BCL6_ and _MYC_. 
Similar events have been previously described in BL and DLBCL
[@Sandlund2006-tm; @Wang2007-lf; @Bertrand2007-nx; @Huret2009-lo].
Interestingly, the breakpoints for each translocation flank focal gains at
each locus.
Based on the orientation of the translocations, these regions may form part
of a double minute consisting of the 3' portion of _MYC_ (after codon 79) and 
a portion of the _BCL6_ super-enhancer.
The effect of truncating the first 79 codons of _MYC_ is not clear though.
A deletion mutant screen showed that deleting codons 3–53 resulted in a 
partially functioning protein, while deleting codons 7–91 led to inactivation
of _MYC_ [@Stone1987-pt].
_MYC_ expression is still significantly higher than the normal germinal center
cells (fold change `r num$myc_fc_bl123`).
We also note the uneven coverage across the exons in the RNA-seq data, with a
sharp increase after codon 79 (Supplementary Figure 
\@ref(fig:custom-myc-bl123-codon79)).
The remaining case (BL277) harbors a complex three-way rearrangement 
(t(8;14;17)(q24;q32;q24)) involving _MYC_, _IGH_, and a locus proximal to 
_GNA13_ on chromosome 17.
The fold change in _MYC_ expression over that in centroblasts and centrocytes 
is also high (fold change `r num$myc_fc_bl277`).
We also inspected the expression of genes within one megabase of the 
breakpoint on chromosome 17 (Supplementary Figure
\@ref(fig:boxplot-chr17gene-zscore-bl277)).
_GNA13_ is the only gene in this region with an established role in BL, a
_bona fide_ tumor suppressor gene.
While its expression is paradoxically high in this tumor, a somatic splice site
mutation that causes the skipping of exon 2 has a variant allele fraction
of 86%, indicating that this region has underwent loss-of-heterozygosity.
The RNA-seq data confirms that most _GNA13_ transcripts lack exon 2 and are 
thus predicted to be non-functional (data not shown).

Despite _MYC_ being deregulated in all tumors by structural alterations, one 
EBV-positive endemic case harbored both an _IGH_-_MYC_ translocation and a 
chromosome 11q aberration with the characteristic proximal gain and distal loss.
It is not possible to determine the chronology of these structural variations
with the available data, which would determine whether the _MYC_ translocation
superseded the 11q aberration.
The 11q abnormality has been previously associated with _MYC_-negative cases
[@Salaverria2014-up] and prompted the creation of a new entity, Burkitt-like
lymphoma with 11q aberration, in the 2016 revision of the World Health
Organization (WHO) classification of lymphoid neoplasms [@Swerdlow2017-rn].
Recent work has demonstrated the existence of sporadic BL tumors harboring both
events [@Grygalewicz2017-cq; @Rymkiewicz2018-zf], and we extend this
observation to an EBV-positive endemic case.


## IGH gene usage for productive allele

<!---
TODO
- Check if there are datasets available for normal B-cell repertoires so we can
  assign P-values to the enrichment in certain genes.
--->

Given the important role for the B-cell receptor (BCR) in BL, we sought to
characterize the immunoglobulin (Ig) heavy (_IGH_) and light (_IGK_/_IGL_)
chains of the productive (_i.e._ expressed) allele.
The MiXCR tool was applied to the RNA-seq data from our BLGSP cohort in 
sensitive mode in order to identify Ig clonotypes 
[@Bolotin2017-hm; @Bolotin2015-az].
Dominant clones were defined as those with a clonal fraction above 30% with at
least 30 supporting reads to minimize spurious results (Supplementary Figure
\@ref(fig:stripplot-sample-clonefraction-isdominant)).
Following this approach, we identified a dominant heavy chain clone in
`r num$ig_heavy` cases (`r num$ig_heavy_percent`%) and a dominant light chain
clone in `r num$ig_light` cases (`r num$ig_light_percent`%). 
There are `r num$ig_heavy_mult` cases with more than one dominant heavy chain
clone. 
On the other hand, there are `r num$ig_light_mult_diff` cases that co-express 
a dominant _IGK_ clone and a dominant _IGL_ clone.
It is not clear whether this heterogenous expression pattern is due to tumor 
cells expressing both _IGK_ and _IGL_ light chains.
Ig gene usage in BL is non-uniform (Figure \@ref(fig:barplot-igvgenes-freq)), 
providing evidence for BCR steoreotypy.
There is no differential gene usage according to clinical variant or EBV
infection status (BH-adjusted P-value < 0.1, Fisher's exact test). 
We replicate previous work showing an enrichment for _IGHV4-34_ and _IGHV3-30_
among _IGHV_ genes [@Lombardo2017-dq; @Amato2016-jf; @Baptista2014-sh].
Furthermore, there is a striking enrichment for _IGKV3-20_ among _IGKV_ genes.
_IGKV3-20_ stereotypy has been described in other B-cell lymphoproliferations,
commonly in combination with _IGHV4-59_, and has been linked to rheumatoid 
factor activity [@Kostareli2012-pf].
While we see that combination in `r clone_pairs_df["IGHV4-59,IGKV3-20",]` 
cases, we observe _IGKV3-20_ pairing with _IGHV4-34_ more commonly in
`r clone_pairs_df["IGHV4-34,IGKV3-20",]` cases (Supplementary Table 
\@ref(tab:table-igheavylightchain-combinations)). 
If _IGKV3-20_ stereotypy is a result of its auto-reactivity, this would suggest
that antigen-dependent BCR signaling may complement antigen-independent (tonic)
signaling in promoting BL cell growth and survival [@Schmitz2014-tz].


(ref:caption-barplot-igvgenes-freq) – Frequency of immunoglobulin V gene usage among dominant heavy (_IGH_) and light (_IGK_/_IGL_) chain clones. Only displaying genes with a minimum frequency of four. No Ig genes were significantly associated with either clinical variant status, EBV infection status or EBV genome type.

```{r barplot-igvgenes-freq, fig.cap="(ref:caption-barplot-igvgenes-freq)", fig.height=7}
clones[is_dominant & class == "IG"] %>%
  .[, .(N = 1), by = .(sample, clinical_variant, ebv_status, 
                       ebv_type, topVHit)] %>%
  .[clinical_variant %in% c("Endemic", "Sporadic")] %>%
  .[, .SD[.N >= 4], by = topVHit] %>% 
  dcast.data.table(
      sample + clinical_variant + ebv_status + ebv_type ~ topVHit, 
      value.var = "N", fill = 0) %>% 
  melt.data.table(
      id.vars = c("sample", "clinical_variant", "ebv_status", "ebv_type"),
      variable.name = "topVHit", value.name = "N") %>% 
  melt.data.table(
      id.vars = c("sample", "topVHit", "N"),
      measure.vars = c("clinical_variant", "ebv_status", "ebv_type"),
      variable.name = "classification", value.name = "subtype") %>%
  .[!(classification == "ebv_type" & subtype == "None")] %>% {
      pvals <- .[, .(pval = fisher.test(as.matrix(xtabs(
                                ~ subtype + N, .SD)))$p.value), 
                 by = .(classification, topVHit)]
      qvals <- pvals[, qval := p.adjust(pval, "BH")]
      merge(., qvals, by = c("classification", "topVHit"))
  } %>% 
  .[, signif := qval < 0.1] %>%
  .[, subtype := fct_relevel(subtype, "Sporadic")] %>%
  .[, classification := fct_recode(classification,
                                   `Clinical variant` = "clinical_variant",
                                   `EBV status` = "ebv_status",
                                   `EBV type` = "ebv_type")] %>%
  .[N > 0] %>% 
  .[, chain := substr(topVHit, 1, 3)] %>% 
  .[, topVHit := fct_infreq(topVHit)] %>% 
  ggplot(aes(topVHit, fill = subtype, color = signif)) +
  geom_bar() +
  facet_grid(classification ~ chain, scales = "free", space = "free_x") +
  scale_fill_manual(
    values = flatten_chr(colours), 
    guide = guide_legend(order = 1, title.position = "top")) +
  scale_color_manual(
    values = c(`TRUE` = "grey0", `FALSE` = "grey60"),
    breaks = c("TRUE", "FALSE"),
    limits = c("TRUE", "FALSE"),
    labels = c("BH-adj. P-value < 0.1", "Not significant"),
    guide = guide_legend(nrow = 2, order = 2, title.position = "top",
      override.aes = list(fill = colours$ebv_type["Type 1"]))) +
  place_legend("bottom") +
  rotate_x_text() +
  labs(x = NULL, y = "Frequency", fill = "Clinical variant", 
       color = "Fisher's exact test")
```


## Mutation burden varies according to EBV infection status

<!---
TODO
- Attempt to explain mutation burden using AID expression, mutation signatures
  and telomere content.
- Double-check what Kaymaz did when comparing type 1 and 2.
--->

We compared the mutation burden between different BL subtypes to identify 
potential differences with respect to mutational processes. 
We only considered the BLGSP WGS data to avoid potential confounding 
effects due to differences in sequencing depth with the ICGC data.
Genome-wide, the mutation load is significantly higher in endemic or
EBV-positive cases compared to sporadic or EBV-negative cases, respectively
(Figure \@ref(fig:boxplot-subtype-mutcount-covariate-mutationtype-wilcoxon)).
The same is true for non-synonymous mutations in all protein-coding genes.
AID-mediated SHM might be the underlying cause for this significantly higher
mutation burden in EBV-positive tumors.
Interestingly, when we only consider non-synonymous mutations in significantly
mutated genes as described below, the direction of the difference is reversed
and remains significant only according to EBV infection status. 
This is despite the higher mutation burden genome-wide in EBV-positive cases.
We were unable to replicate a previously reported significant difference in the 
number of mutations in BL genes according to EBV genome type [@Kaymaz2017-xm]. 
This result supports an oncogenic role for EBV regardless of genome type and 
reinforces the notion that EBV infection status may be a more biologically 
meaningful method for classifying BL than geographic origin.


(ref:caption-boxplot-subtype-mutcount-covariate-mutationtype-wilcoxon) – Frequency of mutations genome-wide, non-synonymous mutations genome-wide, and putative driver mutations in significantly mutated genes in various disease subtypes.

```{r boxplot-subtype-mutcount-covariate-mutationtype-wilcoxon, fig.cap="(ref:caption-boxplot-subtype-mutcount-covariate-mutationtype-wilcoxon)", fig.width=8, fig.height=7}
maf %>%
  .[, .(num_all = .N, num_nonsyn = sum(Protein_Effect == "Nonsynonymous"),
        num_nonsyn_smg = sum(Protein_Effect == "Nonsynonymous" &
                               Hugo_Symbol %in% genes$smgs)),
    by = .(sample, project, clinical_variant, ebv_status, ebv_type, sex)] %>%
  .[clinical_variant %in% c("Endemic", "Sporadic") & project == "BLGSP"] %>%
  .[num_all < min(boxplot.stats(num_all, coef = 5)$out)] %>%
  melt(id.vars = c("sample", "num_all", "num_nonsyn", "num_nonsyn_smg"),
    measure.vars = c("clinical_variant", "ebv_status", "ebv_type", "sex"),
    variable.name = "classification", value.name = "subtype") %>%
  melt.data.table(id.vars = c("sample", "classification", "subtype"),
                  measure.vars = c("num_all", "num_nonsyn", "num_nonsyn_smg"),
                  variable.name = "mutation_type", value.name = "count") %>%
  .[!(classification == "ebv_type" & subtype == "None")] %>% {
    pvals <- .[, .(pval = wilcox.test(count ~ subtype, .SD)$p.value),
               by = .(mutation_type, classification)]
    pvals <- pvals[, qval := p.adjust(pval, "BH")]
    merge(., pvals, by = c("mutation_type", "classification"))
  } %>%
  .[, signif := qval < 0.1] %>%
  .[, subtype := fct_relevel(subtype, "Sporadic")] %>%
  .[, mutation_type := fct_recode(mutation_type,
    `All mutations` = "num_all",
    `Non-synonymous\nmutations` = "num_nonsyn",
    `Non-synonymous\nmutations in BL genes` = "num_nonsyn_smg")] %>%
  .[, classification := fct_recode(classification,
    `Clinical variant` = "clinical_variant",
    `EBV status` = "ebv_status",
    `EBV type` = "ebv_type",
    `Sex` = "sex")] %>%
  ggplot(aes(subtype, count, fill = subtype, alpha = signif, 
             colour = signif)) +
  geom_boxplot() +
  facet_grid(mutation_type ~ classification, scales = "free", 
             space = "free_x") +
  scale_fill_manual(values = flatten_chr(colours), breaks = NULL) +
  scale_alpha_manual(values = c(`TRUE` = 1, `FALSE` = 0.3),
                     breaks = NULL) +
  scale_color_manual(values = c(`TRUE` = "grey0", `FALSE` = "grey60"),
                     breaks = c("TRUE", "FALSE"),
                     labels = c("BH-adj. P-value < 0.1", "Not significant"),
                     guide = guide_legend(nrow = 1, title.position = "top")) +
  place_legend("bottom") +
  labs(x = NULL, y = "Mutation frequency", color = "Wilcoxon test")
```

<!---
Write this section once you have updated mutation signature results.


**Figures**

\@ref(fig:heatmap-exposurefraction-sample-signature)

\@ref(fig:boxplot-signatureexposure-clinicalvariant-wilcoxontest)


**Supplementary Figures**

\@ref(fig:corrplot-autocorrelation-wtsimutsigs)

\@ref(fig:barplot-mutationtype-proportion-signatures)


(ref:caption-heatmap-exposurefraction-sample-signature) – Exposure fractions across the cohort.

```{r heatmap-exposurefraction-sample-signature, fig.cap="(ref:caption-heatmap-exposurefraction-sample-signature)"}
wtsi$exposures_norm %>%
  pheatmap::pheatmap(
    cluster_rows = FALSE, annotation_col = annotations$genome,
    annotation_colors = colours, show_colnames = FALSE,
    color = colorRampPalette(brewer.pal(10, "Spectral")[9:2])(100))
```


(ref:caption-boxplot-signatureexposure-clinicalvariant-wilcoxontest) – A figure showing the differential exposure of mutation signatures in BL subtypes (clinical variant status or EBV infection status).

```{r boxplot-signatureexposure-clinicalvariant-wilcoxontest, fig.cap="(ref:caption-boxplot-signatureexposure-clinicalvariant-wilcoxontest)", fig.width=8, fig.height=7}
wtsi$exposures %>%
  .[n_signatures == selected_sig] %>%
  inner_join(metadata$genome, by = c("sample", "library_id")) %>%
  filter(clinical_variant %in% c("Endemic", "Sporadic")) %>%
  select(sample, signature, exposure, clinical_variant,
         ebv_status, ebv_type, sex) %>%
  gather("classification", "subtype", clinical_variant,
         ebv_status, ebv_type, sex) %>%
  filter(!(classification == "ebv_type" & subtype == "None")) %>% {
    pvals <-
      group_by(., signature, classification) %>%
      summarise(pval = wilcox.test(
        exposure ~ subtype, data = data.frame(exposure, subtype))$p.value) %>%
      mutate(qval = p.adjust(pval, "BH"),
             signif = qval < 0.1)
    inner_join(., pvals, by = c("signature", "classification"))} %>%
  mutate(subtype := fct_relevel(subtype, "Sporadic"),
         classification := fct_recode(classification,
           `Clinical variant` = "clinical_variant",
           `EBV status` = "ebv_status",
           `EBV type` = "ebv_type",
           `Sex` = "sex")) %>%
  ggplot(aes(subtype, exposure, fill = subtype,
             alpha = signif, colour = signif)) +
  geom_boxplot() +
  facet_grid(signature ~ classification, scales = "free", space = "free_x") +
  scale_fill_manual(values = flatten_chr(colours), breaks = NULL) +
  scale_alpha_manual(values = c(`TRUE` = 1, `FALSE` = 0.3),
                     breaks = NULL) +
  scale_color_manual(values = c(`TRUE` = "grey0", `FALSE` = "grey50"),
                     breaks = c("TRUE", "FALSE"),
                     labels = c("BH-adj. P-value < 0.1", "Not significant"),
                     guide = guide_legend(nrow = 1)) +
  # scale_y_log10() +
  place_legend("bottom") +
  labs(x = NULL, y = "Signature exposure", title = "Sample covariate",
       color = NULL)
```
--->

## Significantly mutated genes associate with EBV infection status

<!---
TODO
- Look at all possible pairs for mutual exclusivity.
- Also, perform differential mutation rate analysis while considering other
  mutation types.
- Mention other genes that are not mutated in our cohort
--->

We identified significantly mutated genes (SMGs) by employing an ensemble
approach that leveraged four methods and required support from at least
two [@Lawrence2013-ph; @Tamborero2013-xi; @Mularoni2016-ze;
@Gonzalez-Perez2012-jj] (Figure \@ref(fig:custom-oncoprint)).
Established BL genes figured on our list of SMGs, such as _DDX3X_, _ID3_,
_TCF3_, _ARID1A_, _SMARCA4_, _TP53_, _FOXO1_, and _CCND3_.
We were also able to replicate recently uncovered BL genes such as _TFAP4_ and
_KMT2D_.
On the other hand, we note the absence of somatic mutations in _CCNF_,
which has been described in endemic BL [@Abate2015-iu].
The previously described hotspot mutation (R451C) exists in two endemic cases 
but as germline events, and it also appears as a rare allele in the African 
1000 Genomes Project population in dbSNP (rs150447372).
In addition, our dataset has enabled the identification of genes that have
not yet been described as significantly mutated in BL, namely _SIN3A_,
_HIST1H1E_, _USP7_, _CHD8_, and _RFX7_.
Some of these novel genes carry out functions that are relevant to BL biology.
Notably, _SIN3A_ is a known antagonist of _MYC_ activity [@Nascimento2011-xq].
Thus, it is a likely tumor suppressor gene, especially given its diffuse 
pattern of truncating mutations that are predicted to deactivate the protein.
_USP7_ encodes a deubiquitinase that counteracts Mdm2-mediated ubiquitination
and degradation of p53 [@Li2002-rk].
Mutations in _USP7_ are nearly all truncating, suggesting it is another tumor 
suppressor gene involved in BL. 
Only a single case houses mutations in both _TP53_ and _USP7_, indicating that 
mutations in _USP7_ potentially represent an alternative mechanism for 
abrogating apoptosis.
Interestingly, the protein encoded by the EBV _EBNA1_ gene is known to disrupt 
the interaction between USP7 and p53 [@Holowaty2004-us].
Considering that mutations in _TP53_ and _USP7_ are significantly enriched in 
EBV-negative cases (Figure \@ref(fig:heatmap-diffmutrate-samplecovariate-gene)),
our data supports a role for EBV in inhibiting p53-mediated apoptosis.
Finally, _HIST1H1E_ encodes a histone component, and _CHD8_ is involved in
chromatin remodeling, contributing to the recurrent theme of epigenetic
disruption in B-cell lymphomagenesis.


(ref:caption-custom-oncoprint) – Oncoprint plot showing putative driver mutations in significantly mutated genes (rows) across the BL cohort. Mutations are labeled according to their predicted consequence on the protein.

```{r custom-oncoprint, fig.cap="(ref:caption-custom-oncoprint)"}
include_graphics("images/oncoprint.png")
```

<!---
DDX3X and FOXO11


(ref:caption-custom-gnai2) – Clustering of mutations in GNAI2-encoded protein structure.

```{r custom-gnai2, fig.cap="(ref:caption-custom-gnai2)"}
include_graphics("images/gnai2.png")
```
--->

In an attempt to infer biological differences between BL subtypes, many 
previous studies have compared the mutation rates of various BL genes
according to classification systems such as clinical variant status, EBV 
infection status, and more recently, EBV genome type.
However, these studies were ill-equipped to address this question because of 
small cohort sizes and/or issues relating to RNA-seq somatic variant calls.
On the other hand, our dataset is amenable for these types of analyses given 
its size and high-quality somatic variant calls.
Indeed, we were able to identify a number of SMGs are significantly 
differentially mutated when comparing either endemic and sporadic cases or 
EBV-positive and -negative cases (BH-adjusted P-value < 0.1, Fisher's exact 
test).
However, the only differences that survive multiple test correction are based 
on EBV infection status. 
In all cases, the genes are more commonly mutated in EBV-negative tumours:
_SMARCA4_, _ID3_, _CCND3_, and _TP53_ (with and without _USP7_).
This is liekly the driving force between the significant difference in the 
number of putative driver mutations according to EBV infection status presented
earlier (Figure 
\@ref(fig:boxplot-subtype-mutcount-covariate-mutationtype-wilcoxon)).
We did not identify any difference in gene-wise mutation rates between EBV type
1 and EBV type 2, which departs from what has been previously reported
[@Kaymaz2017-xm]. 


(ref:caption-heatmap-diffmutrate-samplecovariate-gene) – Differential mutation rates among significantly mutated genes in BL. Only showing genes with a significant differences according to at least one classification system. The labels indicate the subtype in which mutations are more common.

```{r heatmap-diffmutrate-samplecovariate-gene, fig.cap="(ref:caption-heatmap-diffmutrate-samplecovariate-gene)", fig.height=4.5, fig.width=7}
smgs_assoc[, .SD[any(pval < 0.05)], by = .(gene)] %>% 
  mutate(font = case_when(signif == "BH-adj. P-value < 0.1" ~ 2, TRUE ~ 1)) %>%
  ggplot(aes(variable, gene, fill = signif, label = higher_in)) +
  geom_tile(color = "grey99", size = 0.5) +
  geom_text(aes(alpha = signif, fontface = font)) +
  scale_x_discrete(
    limits = c("ebv_status", "clinical_variant", "ebv_type", "sex"),
    labels = c("EBV\nstatus", "Clinical\nvariant", "EBV\ntype", "Sex")) +
  scale_fill_manual(guide = guide_legend(title.position = "top"),
    values = c(`BH-adj. P-value < 0.1` = "grey60", `P-value < 0.05` = "grey82",
               "Not significant" = "grey95")) +
  scale_alpha_manual(values = c(`BH-adj. P-value < 0.1` = 1,
                                `P-value < 0.05` = 0.8, "Not significant" = 0),
                     breaks = NULL) +
  place_legend("bottom") +
  labs(x = NULL, y = "Gene/Module", fill = "Fisher's exact test")
```


# Summary and Future Directions

The BLGSP is an ongoing international collaborative project aimed at providing 
a comprehensive molecular portrait of BL across all subtypes.
Here, we identify novel candidate coding and non-coding driver mutations, 
nominating compelling targets for therapy. 
Perhaps the most promising candidate is _USP7_ given the putative role of 
mutations in this gene in abrogating apoptosis. 
We hypothesize that Mdm2 inhibitors may trigger synthetic lethality in 
_USP7_-mutated, _TP53_-wildtype tumors. 
In fact, others have shown that several BL cell lines are sensitive to Mdm2 
antagonists such as nutlin-3 [@Renouf2009-dk].
While only `r num$usp7_muts` cases (`r num$usp7_muts_percent`%) are affected by
_USP7_ mutations, the enrichment of mutations in _TP53_ and _USP7_ among 
EBV-negative cases raises the possibility that all EBV-positive, 
_TP53_-wildtype cases could be treated using the same strategy.
In this work, we also highlight several molecular differences between BL 
subtypes and show that classifying according to EBV infection status has a 
firmer basis in biology.
This is supported by differences in AID-mediated aberrant SHM, AID expression, 
germinal center cell-of-origin, and genome-wide and gene-wise mutation rates. 
These findings put into question the current classification system that relies
on geographic origin by showing that there are more similarities between cases
that are EBV-positive, regardless of where they arise.

Going forward, this project will continue to sequence additional BL cases, with
a focus on adult and HIV-positive patients.
This will allow us to deepen our understanding on the differences between
pediatric and adult BL, if there are any, as well as attempt to delineate the
role for HIV infection in promoting BL formation.
Given the increasing fraction of tumor samples being derived from
formalin-fixed, parafin-embedded (FFPE) tissue, the project also sequenced the 
paired fresh frozen (FF) and FFPE tissues from 35 tumors. 
If I consider the FF data as ground truth, this data will allow me to better 
characterize FFPE artifacts in somatic mutation calls as well as in gene 
expression data. 
To address this, I plan on leveraging machine learning methods to identify
features that discriminate between real variants (_i.e._ found in FF sample) 
and FFPE artifacts.
I will also explore different methods for correcting the difference between the
expression profiles from FF and FFPE tissues for the same underlying sample.
Currently, my timeline for defending my PhD thesis is 1.5 years. 


<!---
# Discussion

To be done.

The objective is to define molecular features that drive lymphomagenesis and 
distinguish disease subtypes, which can be translated to new therapeutic 
strategies deployable worldwide. 


# Acknowledgments

  * ICGC MMML-Seq Project
  * CIHR Travel Award
  * CCSRI Travel Award
--->


##### Page break

# Supplementary Material

## Supplementary Figures

(ref:caption-histogram-distfromtss) – Distance between cluster of non-coding mutations and nearest transcription start site (TSS). Only considering clusters within `r dist_tss_limit` bp of a TSS, representing `r num$wavelets_in_plot_percent`% of all clusters.

```{r histogram-distfromtss, fig.cap="(ref:caption-histogram-distfromtss)"}
wavelets %>%
  filter(abs(dist_tss) < dist_tss_limit,
         numberPatients >= 10) %>%
  ggplot(aes(dist_tss)) +
  geom_rect(aes(xmin = -dist_tss_limit, xmax = 0,
                ymin = 0, ymax = 80), fill = "#D2EEF7") +
  geom_rect(aes(xmin = 0, xmax = dist_tss_limit,
                ymin = 0, ymax = 80), fill = "#C8F7C6") +
  geom_histogram(binwidth = 2000, boundary = 0,
                 closed = "left", fill = "grey50") +
  geom_vline(xintercept = 0, size = 0.5, color = "black", linetype = 2) +
  coord_cartesian(xlim = c(-dist_tss_limit, dist_tss_limit)) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  labs(x = "Distance from nearest TSS", y = "Frequency")
```

##### Page break


(ref:caption-boxplot-covariate-ashmcount) – Number of mutated clusters of non-coding mutations per sample according to clinical variant or EBV infection status. Physiologic targets of SHM consist of the immunoglobulin genes and _MYC_.

```{r boxplot-covariate-ashmcount, fig.height=4, fig.width=6, fig.cap="(ref:caption-boxplot-covariate-ashmcount)"}
wavelets_num_mutated %>%
  filter(shm_type == "Aberrant SHM") %>% 
  ggplot(aes(value, num_mutated, fill = value, 
             alpha = signif, color = signif)) +
  geom_boxplot() +
  # facet_grid(shm_type ~ variable, scales = "free") +
  scale_fill_manual(values = flatten_chr(colours), breaks = NULL) +
  scale_alpha_manual(values = c(`TRUE` = 1, `FALSE` = 0.3),
                     breaks = NULL) +
  scale_color_manual(values = c(`TRUE` = "grey0", `FALSE` = "grey60"),
                     breaks = c("TRUE", "FALSE"),
                     limits = c("TRUE", "FALSE"),
                     labels = c("BH-adj. P-value < 0.1", "Not significant"),
                     guide = guide_legend(nrow = 1, title.position = "top"),
                     name = "Wilcoxon test") +
  place_legend("bottom") +
  scale_x_discrete(labels = c(None = "Negative")) +
  labs(x = NULL, y = "Number of aberrantly mutated clusters")
```

##### Page break


(ref:caption-boxplot-ebvtype-aidexpr-cv) – Box plots showing higher AID expression in EBV-positive cases regardless of geographic origin (endemic or sporadic).

```{r boxplot-ebvtype-aidexpr-cv, fig.cap="(ref:caption-boxplot-ebvtype-aidexpr-cv)"}
aid_expr %>%
  filter(tissue_status == "Tumor") %>%
  spread(variable, value) %>%
  ggplot(aes(`EBV type`, expr, fill = `EBV type`)) +
  geom_boxplot() +
  facet_grid(~ `Clinical variant` ~ .) +
  scale_fill_manual(values = flatten_chr(colours), breaks = NULL) +
  scale_x_discrete(labels = c(None = "Negative")) +
  labs(y = "AID expression")
```

##### Page break


(ref:caption-boxplot-pax5waveletstatus-expr-nearbygenes-wilcox) – Box plots showing expression of genes near the cluster of non-coding mutations that overlaps a _PAX5_ enhancer as a function of its mutation status. Box plots highlighted in red are significantly differentially expressed (Wilcoxon test).

```{r boxplot-pax5waveletstatus-expr-nearbygenes-wilcox, fig.cap="(ref:caption-boxplot-pax5waveletstatus-expr-nearbygenes-wilcox)"}
plot_nearby_gene_expr(
  dds$raw$vst, gene_annotations_gr,
  wavelets_gr[wavelets_gr$gene_symbol == PAX5],
  wavelets_mat_bin[PAX5,])
```

##### Page break


(ref:caption-boxplot-mir142waveletstatus-expr-nearbygenes-wilcox) – Box plots showing expression of genes near the cluster of non-coding mutations that overlaps miR-142 as a function of its mutation status. Box plots highlighted in red are significantly differentially expressed (Wilcoxon test).

<!---
miR-142 is not being shown because I don't have miRNA data for all cases yet.
--->

```{r boxplot-mir142waveletstatus-expr-nearbygenes-wilcox, fig.cap="(ref:caption-boxplot-mir142waveletstatus-expr-nearbygenes-wilcox)", fig.width=12}
plot_nearby_gene_expr(
  dds$raw$vst, gene_annotations_gr,
  wavelets_gr[wavelets_gr$gene_symbol == "MIR142"],
  wavelets_mat_bin["MIR142",], ntop = 10)
```

##### Page break


(ref:caption-boxplot-pvt1waveletstatus-expr-nearbygenes-wilcox) – Box plots showing expression of genes near the cluster of non-coding mutations that overlaps the _PVT1_ promoter region as a function of its mutation status. Box plots highlighted in red are significantly differentially expressed (Wilcoxon test).

```{r boxplot-pvt1waveletstatus-expr-nearbygenes-wilcox, fig.cap="(ref:caption-boxplot-pvt1waveletstatus-expr-nearbygenes-wilcox)"}
plot_nearby_gene_expr(
  dds$raw$vst, gene_annotations_gr,
  wavelets_gr[wavelets_gr$gene_symbol == "PVT1"],
  wavelets_mat_bin["PVT1",])
```

##### Page break


<!---
(ref:caption-heatmap-mrnaexpr-centroblastscoregenes) – Heatmap showing expression of genes used to calculate centroblast score

```{r heatmap-mrnaexpr-centroblastscoregenes, fig.cap="(ref:caption-heatmap-mrnaexpr-centroblastscoregenes)", eval=FALSE}
coo$order <- coo$ps$sample[order(coo$ps$ps)]
coo$order <- c(grep("^CC", metadata$mrna$sample, value = TRUE), coo$order)
coo$order <- c(coo$order, grep("^CB", metadata$mrna$sample, value = TRUE))
coo$order <- coo$order[coo$order %in% colnames(dds$raw$vst)]

dds$raw$vst[coo$genes$gene, coo$order] %>%
  pheatmap(annotation_col = annotations$mrna, annotation_colors = colours,
           cluster_cols = FALSE)
```


(ref:caption-boxplot-ebvtype-centroblastscore-cv) – Centroblast score is higher in EBV-positive cases in both endemic and sporadic cases.

```{r boxplot-ebvtype-centroblastscore-cv, fig.cap="(ref:caption-boxplot-ebvtype-centroblastscore-cv)"}
coo$ps_df %>%
  filter(project == "BLGSP") %>% 
  spread(variable, value) %>%
  ggplot(aes(`EBV type`, ps, fill = `EBV type`)) +
  geom_boxplot() +
  geom_hline(yintercept = 0, linetype = 2, color = "grey") +
  facet_grid(`Clinical variant` ~ .) +
  scale_fill_manual(values = flatten_chr(colours), breaks = NULL) +
  scale_x_discrete(labels = c(None = "Negative")) +
  coord_cartesian(ylim = c(-1, 1)) +
  labs(y = "Centroblast score")
```


##### Page break
--->


(ref:caption-custom-myc-bl123-codon79) – – Integrative Genomics Viewer (IGV) screenshot of tumor RNA-seq (top panel) and WGS (bottom panel) data for BL123 showing the translocation in _MYC_ at codon 79. The center line indicates the location of the breakpoint. The reads in the bottom track are sorted by the chromosome of the mate to emphasize the breakpoint. Note the increase in coverage in the RNA-seq data following the breakpoint.

```{r custom-myc-bl123-codon79, fig.cap="(ref:caption-custom-myc-bl123-codon79)"}
include_graphics("images/myc_bl123_codon79.png")
```

##### Page break


(ref:caption-boxplot-chr17gene-zscore-bl277) – Box plots showing the expression (Z score) of genes within 1 Mb of the chr17 breakpoint in BL277 involving with _IGH_ and _MYC_. The highlighted points in red display the gene expression in BL277.

```{r boxplot-chr17gene-zscore-bl277, fig.cap="(ref:caption-boxplot-chr17gene-zscore-bl277)"}
dds$raw$vst[genes$chr17_sv,] %>%
  as.data.frame() %>%
  rownames_to_column("gene") %>%
  as.data.table() %>%
  melt.data.table(id.vars = "gene", variable.name = "sample",
                  value.name = "expr") %>%
  .[, .SD[var(expr) / mean(expr) >= 0.01], by = gene] %>%
  .[, BL277 := sample == "BL277"] %>%
  .[, zscore := scale(.SD$expr), by = gene] %>%
  .[, gene := fct_relevel(gene, levels(fct_reorder(gene[BL277],
                                                   zscore[BL277])))] %>% 
  {
    ggplot(., aes(gene, zscore)) +
    geom_boxplot() +
    geom_point(data = .[(BL277)], aes(colour = "BL277")) +
    scale_color_manual(values = "red", name = "Sample") +
    place_legend(c(1, 1)) +
    rotate_x_text() +
    labs(x = "Genes near chr17 breakpoint in BL277",
         y = "Expression (Z score)")
  }
```

##### Page break


(ref:caption-stripplot-sample-clonefraction-isdominant) – Strip plot showing the clonal fraction of immunoglobulin heavy and light chains. The 30% minimum threshold in indicated by a dashed line. Translucent points above the dashed line failed to meet the required 30 supporting reads.

```{r stripplot-sample-clonefraction-isdominant, fig.cap="(ref:caption-stripplot-sample-clonefraction-isdominant)"}
clones[cloneCount > 2 & class == "IG"] %>% 
  .[, sample := fct_rev(fct_reorder(sample, cloneFractionPerChain, max))] %>%
  ggplot(aes(sample, cloneFractionPerChain, color = clinical_variant, 
             alpha = is_dominant)) + 
  geom_point() + 
  geom_hline(yintercept = 0.3, linetype = 2, color = "grey") +
  facet_grid(~ chain) +
  scale_color_manual(
    values = colours$clinical_variant,
    guide = guide_legend(nrow = 2, title.position = "top")) +
  scale_alpha_manual(
    name = "Dominant clone?",
    labels = c("Not dominant", "Dominant"), values = c(0.3, 1),
    guide = guide_legend(nrow = 2, title.position = "top")) +
  scale_x_discrete(labels = NULL) +
  place_legend("bottom") +
  labs(x = "Sample", y = "Clonal fraction", color = "Clinical variant")
```

##### Page break


(ref:caption-barplot-igallgenes-freq) – Ig gene usage distribution for all V, D, J and C genes.

```{r barplot-igallgenes-freq, fig.cap="(ref:caption-barplot-igallgenes-freq)", fig.height=8}
c("IGH", "IGK", "IGL") %>%
  map(function(x) {
    clones[(is_dominant)][chain == x] %>%
      melt.data.table(
        id.vars = c("biospecimen_id", "cloneId"), variable.name = "gene",
        measure.vars = c("topVHit", "topDHit", "topJHit", "topCHit")) %>%
      .[value != ""] %>%
      .[, .SD[.N >= 5], by = .(value)] %>%
      .[, gene := sub("top(.)Hit", "\\1", gene)] %>%
      .[, gene := suppressWarnings(fct_relevel(gene, "V", "D", "J"))] %>%
      ggplot(aes(fct_infreq(value))) +
      geom_bar() +
      facet_grid(~ gene, scales = "free_x", space = "free_x") +
      rotate_x_text() +
      labs(x = NULL, y = "Frequency", title = paste(x, "genes"))
  }) %>%
  plot_grid(plotlist = ., ncol = 1)
```

##### Page break


(ref:caption-scatterplot-ighmexpr-ighg1expr-cv-lmfit) – Anti-correlated _IGHM_ and _IGHG1_ expression.

```{r scatterplot-ighmexpr-ighg1expr-cv-lmfit, fig.cap="(ref:caption-scatterplot-ighmexpr-ighg1expr-cv-lmfit)"}
dds$raw$counts_corr[c("IGHM", "IGHG1"),] %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("sample") %>%
  inner_join(metadata$mrna, by = "sample") %>%
  ggplot(aes(IGHM, IGHG1)) +
  geom_point(aes(colour = clinical_variant)) +
  geom_smooth(method = "lm") +
  scale_colour_manual(values = colours$clinical_variant) +
  place_legend(c(1, 1))
```

##### Page break


(ref:caption-stripplot-mycigpartner-ighg1expr-clinicalvariant) – No association between IGG expression and which immunnoglobulin is rearranged with MYC.

```{r stripplot-mycigpartner-ighg1expr-clinicalvariant, fig.cap="(ref:caption-stripplot-mycigpartner-ighg1expr-clinicalvariant)"}
ighg_expr <- dds$raw$counts_corr["IGHG1",] %>% enframe("patient", "ighg_expr")

myc_summary[grepl("^IG", sv)] %>%
  .[, ig := substr(sv, 1, 3)] %>%
  inner_join(ighg_expr, by = "patient") %>%
  inner_join(metadata$mrna, by = "patient") %>%
  ggplot(aes(ig, ighg_expr, color = clinical_variant)) +
  geom_quasirandom() +
  scale_color_manual(values = colours$clinical_variant) +
  place_legend(c(1, 1)) +
  labs(x = "Immunoglobulin rearranged with MYC", y = "IGHG1 expression",
       color = "Clinical variant")
```

##### Page break

<!---
(ref:caption-scatterplot-reconstructionerror-stability-hasdups-maxaidcor) – Selecting an optimal number of signatures.

```{r scatterplot-reconstructionerror-stability-hasdups-maxaidcor, fig.cap="(ref:caption-scatterplot-reconstructionerror-stability-hasdups-maxaidcor)"}
wtsi$metrics %>%
  ggplot(aes(reconstruction_error, stability, color = has_dups,
             size = max_aid_cor, label = n_signatures)) +
  geom_point() +
  geom_text(nudge_y = 0.015, size = 5) +
  scale_color_manual(labels = c("Does not contain duplicates", "Contains duplicates"),
                     values = c("#d95f02", "#1b9e77")) +
  place_legend(c(1, 0)) +
  labs(x = "Reconstruction error", y = "Stability",
       color = "Sanger signatures", size = "Max correlation with AID expression")
```

##### Page break


(ref:caption-barplot-mutationtype-proportion-signatures) – Mutation signatures breakdown.

```{r barplot-mutationtype-proportion-signatures, fig.cap="(ref:caption-barplot-mutationtype-proportion-signatures)"}
contexts <- c(rep(unique(wtsi$signatures$subtype)[1:16], 3),
              rep(unique(wtsi$signatures$subtype)[17:32], 3))

wtsi$signatures %>%
  .[n_signatures == selected_sig] %>%
  .[, mutation_type := fct_inorder(mutation_type)] %>%
  ggplot() +
  geom_hline(yintercept = c(0, 0.05, 0.1), color = "grey", size = 0.3) +
  geom_col(aes(mutation_type, proportion, fill = type), width = 0.5) +
  geom_text(data = wtsi$top_sanger_sigs[n_signatures == selected_sig],
            mapping = aes(x = 0, y = 0.1, label = paste0("Similar to Sanger ", top_sig)),
            hjust = 0, vjust = 1, nudge_x = 1, nudge_y = -0.005) +
  facet_grid(signature ~ .) +
  scale_x_discrete(labels = contexts) +
  scale_y_continuous(breaks = c(0, 0.05, 0.1), labels = c("0%", "5%", "10%")) +
  scale_fill_manual(values = base_change_colours(),
                    guide = guide_legend(nrow = 1)) +
  rotate_x_text() +
  theme(legend.position = "top", legend.justification = 0.5,
        axis.text.x = element_text(size = 5)) +
  labs(x = "Mutation type", y = "Proportion", fill = "Base change")
```

##### Page break


(ref:caption-corrplot-autocorrelation-wtsimutsigs) – There is minimal correlation amongst the WTSI signatures.

```{r corrplot-autocorrelation-wtsimutsigs, fig.cap="(ref:caption-corrplot-autocorrelation-wtsimutsigs)"}
wtsi$signatures %>%
  .[n_signatures == selected_sig] %>%
  .[, mutations := paste0(type, ":", subtype)] %>%
  dcast(mutations ~ signature, value.var = "proportion") %>%
  df_to_matrix() %>%
  set_colnames(sub("V", "WTSI Signature ", colnames(.))) %>%
  cor(.) %>%
  corrplot(diag = FALSE)
```

##### Page break


(ref:caption-ridgeplot-tvaf-sample-project-fforffpe) – Variant allele fraction (VAF) distribution for SNVs and indels across all tumors.

```{r ridgeplot-tvaf-sample-project-fforffpe, fig.cap="(ref:caption-ridgeplot-tvaf-sample-project-fforffpe)", fig.height=10}
maf %>%
  ggplot(aes(t_vaf, fct_reorder(sample, t_vaf), fill = project)) +
  ggridges::geom_density_ridges(rel_min_height = 0.03) +
  facet_grid(ff_or_ffpe ~ ., scales = "free_y", space = "free_y") +
  scale_x_continuous(limits = c(0, 1)) +
  scale_fill_manual(values = colours$project)
```

##### Page break
--->

## Supplementary Tables

<!---
(ref:caption-table-coo-genes) – Genes used in weighted voting algorithm for COO prediction.

```{r table-coo-genes, results="asis"}
kable(coo$genes, caption = "(ref:caption-table-coo-genes)")
```

##### Page break
--->

(ref:caption-table-igheavylightchain-combinations) – Immunoglobulin heavy and light chain pairs. Only showing those with more than one occurrence.

```{r table-igheavylightchain-combinations}
kable(clone_pairs[N > 1], 
      col.names = c("Heavy IG chain", "Ligh IG chain", "Frequency"),
      caption = "(ref:caption-table-igheavylightchain-combinations)")
```

##### Page break

# References

```{r}
knitr::knit_exit()
```


# Results Outline


### EBV integration sites are rare

  * Only one detectable integration event affecting GAS6 gene
  * No effect on gene expression relative to other cases


### Recurrent copy number variations are minimal

  * Results are subject to change as we improve the methodology
  * Recurrent gains
      * Chromosomes 1q, 7, 12, 13
  * Very few recurrent deletions
      * 11q deletion in MYC-positive case, an event previously attributed to
        MYC-negative BL (Salaverria et al., 2014)
      * A few recurrent LOH events including 1p, 6p, 17p
  * Gene-centric copy number analysis
      * Overlay with significantly mutated gene list
      * miR-17~92
  * Figures
      * We can show the landscape of copy number gains and losses. We could
        split the plot based on clinical variant status, EBV infection status
        or age of diagnosis, depending on whether any events are found to be
        subtype-specific.
      * Additional figures will be included if the gene-centric approach yields
        any interesting results (e.g. focal amplification of miR-17-92).


### Intra-tumor heterogeneity

  * There is no compelling evidence of subclonal heterogeneity from the WGS data
  * To be confirmed using the deeper validation sequencing data


### Phylogenetic analysis of BL-associated EBV genomes

  * Perform a phylogenetic analysis of the EBV genomes from EBV-positive cases
  * Include as many EBV genomes from healthy individuals as possible to
    compare pathogenic EBV versus non-pathogenic EBV
  * Figures
    * A phylogenetic tree annotated with various sample covariates. For
      instance, in this tree, the sample labels are colored based on the
      EBV type (black, type 1; red, type 2).


### Established EBV latency states are not clearly followed in BL

  * No firm conclusions from analyzing EBV gene expression
  * Little heterogeneity in EBV gene expression among the cases
  * Established EBV latency states are not closely followed
    * EBNA1 is expression, as expected
    * LMP1, LMP2A, EBNA3A and EBNA3B are expressed in most samples
  * Would need to evaluate EBV protein levels to be conclusive


(ref:caption-heatmap-expr-patient-ebvlatencygenes)
Heatmap of the expression of latency-associated genes.

```{r heatmap-expr-patient-ebvlatencygenes, fig.cap="(ref:caption-heatmap-expr-patient-ebvlatencygenes)"}
ebv_pos <-
  metadata$mrna[ebv_status == "Positive", sample] %>%
  intersect(colnames(mirna$ebv))

latency_annot <-
  ebv_latency %>%
  mutate_at(vars(-gene), function(x) ifelse(x == 0, "Absent", "Present")) %>%
  as.data.frame() %>%
  column_to_rownames("gene")

latency_colours <-
  names(latency_annot) %>%
  set_names() %>%
  map(~ c(Absent = "#bfbfbf", Present = "#666666"))

ebv_latency %>%
  filter(gene %in% rownames(dds$raw$counts_corr)) %$%
  dds$raw$counts_corr[gene, ebv_pos] %>%
  {log10(1 + .)} %>%
  pheatmap(annotation_col = annotations$mrna,
           annotation_colors = c(colours, latency_colours),
           annotation_row = rev(latency_annot), scale = "none",
           cluster_rows = FALSE, show_rownames = TRUE,
           clustering_distance_cols = "euclidean")
```


(ref:caption-heatmap-expr-patient-ebvgenes)
Heatmap showing the expression of all EBV genes.

```{r heatmap-expr-patient-ebvgenes, fig.cap="(ref:caption-heatmap-expr-patient-ebvgenes)"}
log10(dds$raw$counts_corr[genes$ebv, ebv_pos] + 1) %>%
  pheatmap(scale = "none", annotation_col = annotations$mrna, annotation_colors = colours,
           clustering_distance_cols = "euclidean")
```


### Molecular Burkitt lymphoma

Confirm that BLGSP BL cases are indeed BL by running both mBL classifiers.

* Clustering on Hummel _et al._ genes.

```{r heatmap_expr_samples_mblhummel}
dds$raw$vst[rownames(dds$raw$vst) %in% genes$mbl_hummel,] %>%
  pheatmap(annotation_col = annotations$mrna, annotation_colors = colours,
           border_color = NA)
```

* Clustering on Dave _et al._ genes.

```{r heatmap_expr_samples_mbldave}
dds$raw$vst[rownames(dds$raw$vst) %in% genes$mbl_dave,] %>%
  pheatmap(annotation_col = annotations$mrna, annotation_colors = colours,
           border_color = NA)
```


### miRNA expression analysis



* Identify mRNAs with expression anti-correlated with differentially expressed miRNAs
* A subset of EBV miRNAs are detected in all EBV-positive cases

**Figures**

* First, the results from the miRNA:mRNA anti-correlation analysis can be visualized in a number of ways. This can be discussed once the first results are obtained.

```{r ma_plot_mirna}
mirna$de %>%
  mutate(factor = sub("EBV–", "EBV-", factor),
         factor = fct_relevel(factor, "Tumor vs Normal")) %>%
  ggplot(aes(log2(mean), lfc, color = signif)) +
  geom_point() +
  scale_color_manual(values = c(`TRUE` = "red", `FALSE` = "black"), breaks = NULL) +
  facet_grid(~ factor) +
  coord_cartesian(ylim = c(-10, 10)) +
  labs(x = "miRNA mean expression (log2 RPM)", y = "Log2 Fold Change")
```

```{r heatmap_mirna_expr_de}
mirna$rpms_filt[unique(mirna$de_ac$mirna), ] %>%
  {log10(. + 1)} %>%
  pheatmap(annotation_col = annotations$mirna, show_rownames = TRUE,
           clustering_distance_cols = "euclidean", annotation_color = colours)
```

```{r heatmap_mirna_expr_de_ac}
mirna$rpms_filt[unique(mirna$de_ac[(signif), mirna]), ] %>%
  {log10(. + 1)} %>%
  pheatmap(annotation_col = annotations$mirna, show_rownames = TRUE,
           clustering_distance_cols = "euclidean", annotation_color = colours)
```

* Second, the heatmap of EBV miRNA expression can be shown. The heatmap can also just show EBV miRNAs that are anti-correlated with interesting human mRNA targets.


# Other Results

AID activity is correlated with mutations in DDX3X and BCL7A.

```{r boxplot_mutationstatus_aidactivity_gene}
colSums(wavelets_mat) %>%
  enframe("sample", "aid_activity") %>%
  inner_join(smgs_status, ., by = "sample") %>%
  group_by(gene) %>%
  mutate(pval = wilcox.test(aid_activity[status], aid_activity[!status])$p.value) %>%
  ungroup() %>%
  mutate(gene = fct_reorder(gene, pval, median),
         signif = pval < (0.05 / n_distinct(gene))) %>%
  filter(dense_rank(pval) <= 5) %>%
  ggplot(aes(status, aid_activity, fill = signif)) +
  geom_boxplot() +
  facet_grid(~ gene) +
  scale_x_discrete(labels = c("No", "Yes")) +
  scale_fill_manual(values = c(`TRUE` = "dodgerblue", `FALSE` = "grey"), drop = FALSE,
                    labels = c("Significant", "Not significant"),
                    limits = c("TRUE", "FALSE")) +
  place_legend("top") +
  theme(panel.spacing.x = unit(1.5, "lines")) +
  labs(x = "Is mutated?", y = "AID activity", fill = "Wilcoxon test")
```


Mutations in SMARCA4 is correlated with AID expression.

```{r boxplot_mutationstatus_aidexpr_gene}
dds$raw$counts_corr["AICDA",] %>%
  enframe("sample", "aid_expr") %>%
  inner_join(smgs_status, ., by = "sample") %>%
  group_by(gene) %>%
  mutate(pval = wilcox.test(aid_expr[status], aid_expr[!status])$p.value) %>%
  ungroup() %>%
  mutate(gene = fct_reorder(gene, pval, median),
         signif = pval < (0.05 / n_distinct(gene))) %>%
  filter(dense_rank(pval) <= 5) %>%
  ggplot(aes(status, aid_expr, fill = signif)) +
  geom_boxplot() +
  facet_grid(~ gene) +
  scale_x_discrete(labels = c("No", "Yes")) +
  scale_fill_manual(values = c(`TRUE` = "dodgerblue", `FALSE` = "grey"), drop = FALSE,
                    labels = c("Significant", "Not significant"),
                    limits = c("TRUE", "FALSE")) +
  place_legend("top") +
  theme(panel.spacing.x = unit(1.5, "lines")) +
  labs(x = "Is mutated?", y = "AID expression", fill = "Wilcoxon test")
```


Telomere length does not correlate with mutation burden.

```{r scatterplot_telcontent_mutfreq}
num_mutations <- maf[, .(num_mutations = .N), by = biospecimen_id]

telomeres <- inner_join(telomeres, num_mutations, by = "biospecimen_id")

telomeres %>%
  filter(!num_mutations %in% boxplot.stats(num_mutations, coef = 5)$out) %>%
  ggplot(aes(tumor, num_mutations, color = clinical_variant)) +
  geom_point() +
  scale_colour_manual(values = colours$clinical_variant) +
  place_legend(c(0.5, 1)) +
  labs(x = "Tumor telomere content", y = "Mutation frequency",
       color = "Clinical variant")
```


Heatmap showing association with covariates.

```{r heatmap_mutstatus_patient_smg}
smgs_status %>%
  .[, status := as.integer(status)] %>%
  dcast.data.table(gene ~ sample, value.var = "status") %>%
  df_to_matrix() %>% {
      cv <- annotations$genome[colnames(.), "clinical_variant"]
      ebv <- annotations$genome[colnames(.), "ebv_type"]
      .[, order(ebv, cv)] } %>%
  pheatmap::pheatmap(
    annotation_col = annotations$genome, annotation_colors = colours,
    cluster_cols = FALSE, cluster_rows = FALSE, color = c("grey80", "grey50"),
    show_colnames = FALSE, drop_levels = TRUE)
```


No association between the centroblast score and mutations in SMGs.

```{r boxplot_mutationstatus_centroblastscore_gene}
smgs_status %>%
  inner_join(rownames_to_column(coo$ps, "sample"), by = "sample") %>%
  filter(!is.na(ps)) %>%
  group_by(gene) %>%
  mutate(pval = wilcox.test(ps[status == 1], ps[status == 0])$p.value) %>%
  ungroup() %>%
  mutate(gene = fct_reorder(gene, pval, median),
         signif = pval < (0.05 / n_distinct(gene))) %>%
  filter(dense_rank(pval) <= 5) %>%
  ggplot(aes(as.factor(status), ps, fill = signif)) +
  geom_boxplot() +
  facet_grid(~ gene) +
  scale_fill_manual(values = c(`TRUE` = "dodgerblue", `FALSE` = "grey"), drop = FALSE,
                    labels = c("Significant", "Not significant"),
                    limits = c("TRUE", "FALSE")) +
  theme(panel.spacing.x = unit(1.5, "lines")) +
  scale_x_discrete(labels = c("No", "Yes")) +
  place_legend("top") +
  labs(x = "Is mutated?", y = "Centroblast score", fill = "Wilcoxon test")
```


No association between centroblast score and the AID transcript ratio.

```{r scatterplot_centroblastscore_aidtxratio_clinicalvariant}
aid_tx_ratio <-
  (tx$abundance["ENST00000537228",] / tx$abundance["ENST00000229335",]) %>%
  enframe("sample", "aid_tx_ratio")

coo$ps %>%
  rownames_to_column("sample") %>%
  inner_join(aid_tx_ratio, by = "sample") %>%
  inner_join(metadata$mrna, by = "sample") %>%
  ggplot(aes(ps, aid_tx_ratio)) +
  geom_smooth(method = "lm") +
  geom_quasirandom(aes(color = clinical_variant), groupOnX = FALSE) +
  geom_vline(xintercept = 0, linetype = 2, color = "grey") +
  scale_colour_manual(values = colours$clinical_variant) +
  coord_cartesian(xlim = c(-60, 60)) +
  place_legend(c(0, 1)) +
  labs(x = "Centroblast score", y = "AID transcript ratio (short / long)",
       color = "Clinical variant")
```


No association between AID activity (measured by the number of mutation non-coding peaks) and the AID transcript ratio.

```{r scatterplot_aidactivity_aidtxratio_clinicalvariant}
colSums(wavelets_mat_bin) %>%
  enframe("sample", "aid_activity") %>%
  inner_join(aid_tx_ratio, by = "sample") %>%
  inner_join(metadata$mrna, by = "sample") %>%
  ggplot(aes(aid_activity, aid_tx_ratio)) +
  geom_smooth(method = "lm") +
  geom_quasirandom(aes(color = clinical_variant), groupOnX = FALSE) +
  scale_colour_manual(values = colours$clinical_variant) +
  place_legend(c(1, 1)) +
  labs(x = "AID activity", y = "AID transcript ratio (short / long)",
       color = "Clinical variant")
```


```{r}
maf_noncoding[, .SD[.N >= 10], by = id] %>%
  .[, id_symbol := paste0(symbol, " (", id, ")")] %>%
  .[!is.na(ig_overlap), id_symbol := ig_overlap] %>%
  ggplot(aes(fct_reorder(id_symbol, t_vaf), t_vaf, fill = ebv_status)) +
  geom_boxplot() +
  rotate_x_text() +
  labs(x = "Nearest gene", y = "VAF") +
  scale_fill_manual(values = colours$ebv_status)
```
