Multi-resolution modeling of a discrete stochastic process identifies causes of cancer
Authors: Adam Uri Yaari
Maxwell Sherman
Oliver Clarke Priebe
Po-Ru Loh
Boris Katz
Andrei Barbu
Bonnie Berger
Venue: ICLR 2021 Poster (2021)
Keywords: Computational Biology
non-stationary stochastic processes
cancer research
deep learning
probabelistic models
graphical models

Abstract: Detection of cancer-causing mutations within the vast and mostly unexplored human genome is a major challenge. Doing so requires modeling the background mutation rate, a highly non-stationary stochastic process, across regions of interest varying in size from one to millions of positions. Here, we present the split-Poisson-Gamma (SPG) distribution, an extension of the classical Poisson-Gamma formulation, to model a discrete stochastic process at multiple resolutions. We demonstrate that the probability model has a closed-form posterior, enabling efficient and accurate linear-time prediction over any length scale after the parameters of the model have been inferred a single time. We apply our framework to model mutation rates in tumors and show that model parameters can be accurately inferred from high-dimensional epigenetic data using a convolutional neural network, Gaussian process, and maximum-likelihood estimation. Our method is both more accurate and more efficient than existing models over a large range of length scales. We demonstrate the usefulness of multi-resolution modeling by detecting genomic elements that drive tumor emergence and are of vastly differing sizes.

Published as a conference paper at ICLR 2021
MULTI-RESOLUTION MODELING OF A DISCRETE STOCHASTIC PROCESS IDENTIFIES CAUSES OF CANCER
Adam Ya. Maxwell Sherman 13355 11335 Oliver Priebe
Po-Ru Loh Boris Katz! Andrei Barbul.2 Bonnie Berger
MIT CSAIL, 2 MIT CBMM, 3 Broad Institute of MIT and Harvard, 4 MIT Department of Mathematics,
5 Division of Genetics, Brigham and Women's 's Hospital Department of Physics, University of Pennsylvania
{ Iyarri, maxas, pri ebeo, boris, abarbu, bab @essaill 1 mit edu, poruloh@broadinst orr..e...ttttttttoooooo...g
ABSTRACT
Detection of cancer-causing mutations within the vast and mostly unexplored hu- man genome is a major challenge. Doing so requires modeling the background mutation rate, a highly non-stationary stochastic process, across regions of interest varying in size from one to millions of positions. Here, present the split- Doisson-Gamma (SPG) distribution, an extension of the classical Daammnn formulation, to model discrete stochastic process at multiple resolutions. We demonstrate that the probability model has a closed-form posterior, enabling effi- cient and accurate linear- time prediction over any length scale after the parameters of the model have been inferred a single time. We apply our framework to model mutation rates in tumors and show that model parameters can be accurately in- ferred from high- dimensioini epigenetic data using a convolutional neural network, Gaussian process, and maximum-likeline accurate and more efficient than existing models over a large range of length scales. We demonstrate the usefulness of multi-resolution modeling by detecting genomic elements that drive tumor emergence and are of vastly differing sizes.
estimation. Our method is both more
INTRODUCTION
Numerous domains involve modeling highly non-stationary discrrttetime and integer-valued stochas- tic processes where event counts vary dramatically over time or space. An important open problem of this nature in biology is understanding the stochastic process by which mutations arise across the genome. This is central to identifying mutations that drive cancer emergence (Lawrence et al. 2013).
Tumor drivers provide a cellular growth advantage to cells by altering the function of genomic element such as a gene or regulatory feature (e.g. promoter). Drivers are identifiable because they reoccur across tumors, but there are two major challenges to detecting such recurrence. First, driver mutations are rare and their signal is hidden by the thousands of passenger mutations that passively and stochastically accumulate in tumors (Stratton et al.) 2009: Martincorenaa & Campbell 2015). Second, because functional elements vary dramatically elements: -00-103 bases; and single positions), driver mutations accumulate across regions that vary many orders of magnitude. Accurately predicting the stochastic accumulation of passenger mutations at multiple scales is necessary to reveal the subtle recurrence of driver mutations across the genome.
size (genes: 10- bases; regulatory
Here, we introduce the split-Poisson Gamma (SPG) process, an extension of the Poisson-Gamma distribution, to efficiently model a non-stationary discrete stochastic process at numerous length scales. The model first approximates quasi-stationary regional rate parameters within small windows; it then projects these estimates arbitrary regions in linear time (10-15 minutes for genome wide inference). This approach is in contrast to existing efforts that model fixed regions and require computationally expensive retraining (e.g. over 5 hours) to predict over multiple scales of interest. (Nik-Zainal et al. (2016] Martincorena et al. (2017). We apply our framework to model cancer-specific mutation patterns (fig. D.. We perform data-driven training of our model's parameters and show that it more accurately captures mutation patterns than existing methods on simulated and real data. We demonstrate the power of our multi-resolution approach by identifying drivers across functional
Authors contrihuted equally to this work.

---

Published as a conference paper at ICLR 2021
elements: genes, regulatory features, and single base mutations. Despite the method having knowledge of genome structure, it detects nearly all gene drivers present in over 5% of samples while making no false discoveries and detects all previously characterized regulatory drivers. Detected events also include novel candidate drivers, providing promising targets for future investigation.
Regional mucatin
Gase pal-resslldion dfiver
Crosed drromaon
e
mln
Caophaggu
9kin
OUT
Regional
sedtet &+
CNN
GP
X5
GM
Figure 1: Non-stationary stochastic process modeling predicts mutation patterns and identifies cancer-specific driver mutations. Biological processes are shown in bluc, data processing is shown orange. a. Areas of the genome have varying epigenetic states (c.g. accessibility for transcription) depending on the tissue type. b. These epigenetic states set different mutation rates in different tissues. c. Our model takes these epigenetic tracks as input to estimate the regional mutation density across the genome (95% confidence interval in orange). d. Regional parameters and sequence context integrated via the provide arbitrary resolution mutation count estimates. Deviations between the estimated and observed mutation rates identify mutations that are associated with cancers in different tissues. e. The split-Poisson GGamma (SPG) model plate diagram (squares: inferred parameters; grey: observed input data).
Poisson-Gammer (SPG) distribution
1.1 PREVIOUS WORK
Numerous methods exist for modeling stationary stochastic processes (Lindsey) 2004). Far fewer exist for non-stationary processes because they are difficult to capture with the covariance functions of parametric models (Risser 2016). Non-stationary kernels have been introduced for Gaussian processes (Paciorek & Schervish 2004). but these may not be tractable on large datasets due to their computational complexity. More recently, there has been work developing Poisson-game sassonngammaa models for dynamical systems (Schein et al.] 2016) Guo et al. (018), but these methods have focused on learning relationships between count variables, not predicting counts based on continuous covariates.
In the particular case of modeling mutation patterns across the cancer genome, numerous compu- tational methods exist to model mutation rates within well-understood genomic contexts such as genes (Lawrence et al. 2013 Martincorena et al. 2017 Wadi al. [017 Mularoni et al. 2016
Juul et al.) These models account for < 4% of the genome (Rheinbay et al. (020). They are not applicable in non-coding regions, where the majority of mutations occur (Gloss & Dinger 2018). A
handful of methods to model genome-wide mutation rates have been introduced (Polak et al. 20155]
Nik-Zainal et al., 2016 Bertl et al., 2018). However, they operate on single length-scale of regions and require computationally expensive retraining to predict over each new length-scall Several methods rely on Poisson or binomial regression; however, previous work has extensively documented that mutation counts data are over-dispersed, leading these models to underestimate variance and yield numerous false-positive driver predictions (Lochovsky et al. 2015 Martincorena et al., 1007,Juu et al.] (211)). Negative binomial regression has recently been used to account for over-dispersion (Nik-Zainal et al. (016) and perform genome-wide mutation modeling and driver detection. However, resolution was coarse, and it only found few, highly recurrent driver mutations.
set
1.2 OUR CONTRIBUTIONS
This work makes three key contributions: 1) we introduce extension of the distribution to model non-stationary discrete stochastic processes at any arbitrary length scale without retraining; 2) we apply the framework capture cancer-specific mutation rates with unprecedented accuracy, resolution, and efficiency; and 3) we perform a multi-scale search for cancer driver mutations genome-wide, including the first-ever base-resolution scan of the whole genome. This search yields
2

---

Published as a conference paper at ICLR 2021
several new candidate driver events in the largely unexplored non-coding genome, which working on validating with experimental collaborators. Crucially, our approach allows fast, efficient, and accurate searches for driver elements and mutations anywhere in the genome without requiring arduous retraining of a model, a feat which is not possible with existing approaches.
MULTI-RESOLUTION MODELING OF A NON-STATIONARY DISCRETE STOCHASTIC PROCESS
We consider a non-stationary discrete stochastic process (M;; 1 1.2.... where M, is the integer- valued event count position i. Associated with each position is a real-valued, L-dimensional feature vector 77 that determines the instantaneous event rate 24 via an unknown function. Thus region {i, 1, ....i + N} of N contiguous positions characterized by L N feature matrix nR and an event count XR LieR-M M,. As training data, 7R. Xg. and M, are observed for some set of regions {R E T }. Then given set of feature matrices from unobserved regions { R H). ). the challenge is to predict the distribution of event counts over any arbitrary set of unseen positions that may or may not be contiguous. Real- woaldw examples include traders in a stock market, packets delivered to routers in network, and mutations accumulating at positions in the genome.
2.1 THE OPLIT-POISSON-GOE
PROCESS
We assume that the process is near-stationary within a small enough region R= (i,i+1, .... ... and that the L X N covariate matrix aa is observed. Thus the rate of events de within R is approximately constant and associated with 17R. albeit in an unknown way. A number events (Xx) may occur within R dependent dR and are then stochastically distributed individual positions within R, implying hierarchical factorization of the scalar random variables dg. XR. and M, (fig. [T) as
Pr(M: XR, AR: πR) Pr(M; k k|xx;
ARi MR Pr
XR and de are unknown nuisance variables and are marginalized in general as
Pr(M: = k(n))= 1/' Pr(Ax:TR) nR)
Σ Xg
Pr(M;
(R)
Since applications often require many posterior predictions over regions of varying sizes, we propose a prior parameterization that builds on the success and flexibility of the classical distribution while ensuring the marginalization has an easy-to-compute posterior distribution:
åR Gamma (1R, 0x)
(3)
XR Poissoní AR)
(4)
M; Binomial (XR.
(5)
where and Og are shape and scale parameters dependent TR. Pi is the time-averaged probability of an event at and ýú EjEppj* the normalized probability within R. plate diagram of the hierarchical model is presented in fig- De.
The above formulation provides simple, closed form solution to eq. distribution (See Appendix for details):
P as a negative binomial (NB)
Pr(M; = k| Cp.. 6R. PR.p::7 = NB CRR I O2 p.
(6)
Eq. implicitly assumes that events distributed independently units within R. Exploiting this assumption, eq. (6) immediately generalizes to consider any set of units / C R
Pr Σ M: k OR, (iilie TR NB k: ClR: + OR Eier Pi
The above formulation is an extension of the classical Poisson-Gamma distribution whereby the Poisson is randomly split by binomial. We term this split-Poisson-Gaammee (SPG) process. While
3

---

Published as a conference paper at ICLR 2021
the derivation of the SPG solution makes simplifying assumptions, the benefit is that the parameters OR and 6R need to be estimated only once for each overlappag region R. Estimates for a region of any other size can then be computed in constant time from eq. (7). If a new region R' is larger than R, we approximate the gamma distribution in a super-region containing R' as a superposition of the previously inferred parameters of each region of size R within the super-region (see section 2.2
2.2 INFERRING REGIONAL RATE PARAMETERS
The statistical power of SPG depends on accurate estimation of the regional gamma rate parameters Og and Og- We propose variational approach to enable flexible, accurate, and non-linear inference of these parameters from a set of covariates, Let G( 0) be a gamma distribution. By the central limit theorem, limg- G(a, 0) = N(H, o²) where μ = ab and o² = az².. We thus use a Gaussian process (GP) to non-linearly map covariates to regional estimates for μe and oz The variational estimates for the gamma parameters are then
0g = MR OR For a super-region R' = =R; + Ry. H4 = μRR + μR, and of o + oz.
Og = / o
limitation of this approach is that GPs can only operate on vectors of covariates. Thus a dimension- ality reduction method must be applied to the input matrix TR.- In cases where 77x includes spatial relationships, a convolutional neural network can be a powerful approach to dimension-reduucionn; however, other approaches are feasible (see section 3.2 and section 5.1
2.3 INFERRING AVERAGEE EVENT PROBABILITIES
The time-averaged parameters {pr i= ...) must also be inferred. Crucially, as seen in eq. ( these parameters are never used directly; instead, they are always renormalized sum to one within a region of interest. Thus, estimates do not need to reflect the absolute probability of an event at i but merely the relative rate of events between positions. Indeed, because of the renormalization procedure, the estimates need not even be a true probability distribution. Estimating can thus be accomplished by clustering units with similar relative rates of events. How this clustering should be performed will depend on the application of interest (see section 3.3 for concrete example).
FITTING PARAMETERS TO PREDICT CANCER MUTATION PATTERNS
We obtained publicly available mutation counts from four cancer cohorts previously characterized by the Pan-Cancer Analysis of Whole Genomes Consortium (PCAWG) (Campbell et al., 2020): esophageal adenocarcinoma (N 98 tumors; 2.7M mutations), skin melanoma (N = 70 tumors; n s 7.8M mutations), stomach adenocarcinoma (N = 37 tumors; n ss 480k mutations), and liver hepatocellular carcinoma (N 264 tumors; 3.3M mutations). Crucially, these data contain only the total number of mutations at each position in the genome. We do not know a priori which mutations are background mutations and which are driver mutations. We also do not know the true mean and variance of the underlying mutation rate in any region.
We do know that the mutation rate is highly associated with chemical modifications of the DNA that set the way it is processed in a cell, collectively termed the epigenome (Schuster-Böcce) Böckler & Lehner al. 2015). We obtained 733 datasets characterizing the patterns of these chemical
2012; Polak Polak modifications in 111 human tissues from Roadmap Epigenomics (Roadmap Epigenomics Consortium al., 2015). These data are the largest compendium of uniformly processed human epigenome sequencing currently available. Each track provides the -log P-value that particular modification present at each location of the genome
given tissue type. We additionally created two tracks that provide the average nucleotide and GC content in region based on the human reference genome GRCh37. See Appendix and supplementary data for additional information on the epigenetic tracks. The input matrix for each region l7 thus has 735 rows. We fixed the number of columns to be 100 irrespective of the size of R, where each column is the mean across R/ 100 adjacent positions.
3.1 ARTIFICALL DATASET
In order to evaluate the ability of SPG and other models to estimate the unknown mean and variance of regional rates, created simulated datasets with known mean and variance parameters dependent
4

---

Published as a conference paper at ICLR 2021
b
true mean R Osserve #muratons R?
Esophageal
Melanoma
Stomach
Data simulation
Parameter interence
CNN
5m true varainee 87
6R
OffeGG
NOW
Figure 2: Data simulation and regional parameters inference accuracy across methods. a. Simulated data experiment, Data simulation shown in green, parameter inference is shown in orange. Simulated jup and oR are computed as the mean and variance of sample's 's KNN cluster. The model is trained over randomly sampled event counts; the parameter estimates fäg and og are then compared to their true values. b. Pearson R² to the observed mutation count in the true (unsimulated) data (top) simulated mean (middle) and simulated variance (bottom) for the CNN+GP parameter estimation strategy (CNN+GP), random forest (RF), negative binomial regression (NBR) and binomial regression (BR). Results for additional estimation techniques are in Appendix.
on the observed input matrix (fig.Pp) Db) We created input matrices of size 735 x 100 from the epigenetic tracks (described above) for non-overlapping regions of 50,000 positions. To define non-stationary mean and variance of mutation rate dependent on the each region's 's input matrix, we reduced Te to an feature vector of size 735 by taking the mean across columns and used a k-nearest-neiihhors (KNN) strategy to identify 500 regions with similar epigenetic feature vectors; we then defined μl and o for each region as the mean and variance of the observed event counts across its 500 neighbouring regions. The number of observed events for that region was then randomly drawn from a negative binomial distribution defined by those parameters (full technical details in Appendix). Models were trained on the randomly drawn counts and evaluated on their ability to accurately infer the true mean and variance. We simulated 50kb regions following previous work (Rheinbay et al. 2020).
3.2 ESTIMATING DYNAMIC REGIONAL RATES WITH UNCERTAINTY
The input matrices Me E R 105 required significant dimension reduction before we could employ our GP-based variational strategy to infer SPG regional rate parameters. Columns encode the high- resolution spatial organization
the epigenome which have recently been shown to be important determinants of local mutation rate (Gonzale7-Perez et al., 2019 Akkemir et al. 2020). Therefore, we hypothesized that a convolutional neural network (CNN) would provide a powerful approach to produce -indeenniiooal embedding that retrains information about this local structure; the supervised nature a CNN further enables the resulting embedding to be optimized for the cancer of interest, which is crucial to performance since the epigenetic determinants of mutation rate vary drastically between cancer types (Polak et al.) 2015). We constructed a ID CNN model with ceattiionnrat--sscciate local epigenetic patterns
residual blocks and 3 fully-connected layers to map to regional mutation rates. The CNN non-linearly reduces πR E R 735 100 to a 16 dimensional feature vector its last feature layer. The CNN was trained to minimize the mean squared error between observed and predicted mutation counts. Due to the interchangeable nature of the rows, the ID kernels allow the network to identify arbitrary inter-track interactions. The final 16-dimension feature vector was then passed as input to sparse GP (Titsias), fit to maximize the likelihood of the observed mutation counts using 2000 inducing points and a radial basis function kernel (fig. [T).. We found that results were robust to the particular choice of kernel and hyperpriors placed over kernel parameters. While end-to-end training is possible (Bradshaw et al. 2017). we did not find necessary to achieve high accuracy in this particular application. A CNN is not the only method available to reduce dimensionality prior to GP inference; we investigated numerous other methods, but found the CNN+GP to produce the most accurate results (see Appendix).
5

---

Published as a conference paper at ICLR 2021
3.3 ESTIMATING TIME- AVERAGED EVENT PROBABILITIES
In the case of cancer mutation patterns, previous work showed that the mutation rate at any position i is heavily influenced by the nucleotide at and the two nucleotides directly adjacent to i; positions with this same "trinucleotide context" will have similar mutation patterns (Alexandrov et al. (013)
Following previous works (Mularoni et al., 2016 Wadi al. 2017 Martincorena et al. 2017
Weghorn & Sunyaev 2017). we used trinucleotide context to estimate Let ntn' be the trinuclcotide context centered position i. We estimate the probability that is mutated using the ensemble maximum-like likelihoo estimate of its cluster
Vnna
Pi Pna.n Nuzaa
where Nusw is the number of ntn' trinucleotides in the genome and Van is the number of times mutated within nt'. This approach alone explains little variance in sub-megabase regions (see Appendix) because it does not account for regional mutation rates.
3.3.1 COMPARING TO BENCHMARK MODELS
We compared SPG to three alternative approaches that have been previously used to learn both the mean and variance of regional mutation patterns genome-wide. The alternative models random forest (RF) regression (Polak et al. 2013), binomial regression (BR) (Bertl et al.] (2018). and negative binomial regression (NBR) (Nik-Zainal et al. 2016 Martincorena et al. (001)). For the RF, we used the Jackknife method (Wager et al.) to estimate the variance; this method requires (n)) trees where is the number of samples in the training set. BR and NBR directly specify the
variance as function of the mean: BR as OR μ 222 /n and NBR as o # Buue), where is an overdispersion parameter. Benchmarking comparisons were performed on the skin melanoma, esophageal adenocarcinoma, and stomach adenocareinoma cohorts.
3.4 MODEL TRAINING
For every region of size R, epigenetic features were extracted into matrices of size 735 tracks by 100 binned position columns, where each column was the mean across R 100 adjacent base-pairs. Regions with highly repetitive DNA sequence (<70% of 36mer sub-sequences being unique) where excluded from the training set to ensure high data quality as in previous analyses (Polak et al. 2015). Before training, high-quality data regions were strictly split into train (64%), validation (16%) and test (20%) sets. Predictions for excluded regions and held-out test sets were obtained after model training. Genome-wide predictions were generated using 5-fold coossvvlidattion The CNN received the full 735 x 100 matrices as input. Vector-based methods (RF, NBR, BR) received the 735-dimension vector of epigenetic values averaged across position columns. Following previous work, also included the expected number of mutations based on the trinuclcotide composition of region as an offset term in NBR and BR when predicting mutation counts NNi--aiial et al. 2016). Additional details on training (e.g. number of epochs) are in Appendix.
IDENTIFYING GENETIC DRIVERS OF CANCER
Because cancer drivers reoccur across tumors, driver elements (genes, regulatory structures, and individual base-pairs) will contain an excess of mutations relative to the of expected background mutations. The SPG model provides simple, efficient, and accurate method to search for this recurrence. We first estimate mean and variance of the background mutation rate using the CNN+GP estimation method. We then apply eq. (7) to search for statistical evidence that the number of observed mutations, k, exceeds expectation within every gene, known regulatory structure, and 50 bp window in the genome by changing the set of tested positions I. For gene, is the number of observed missense or nonsense mutations and is the set of all possible mutations in the gene. For both a regulatory element and window of fixed size, is the number of mutations observed in the element / window and / is the set of all positions within the element/ / window. If an element overlaps multiple 10kb regions, we merge the mean and variance estimates for overlapped regions as described in section 2 To maintain strict train-test separation, both the rate parameters and pı are estimated excluding the element being tested. We controlled family-wise error rate at the a = 0.05 level using Bonferroni correction for the total number of tests in genes, regulatory elements, 50bp windows.

---

Published as a conference paper at ICLR 2021
Esophagea
3.8
Vartance aaplareed
0.8
0.0
NBR
0.6
0.6
0.6
0.4-
0.4
0.4
0.2
0.2
0.2
10bp 20b0p 1kb 13ke 100ko
1M1
1000
10000 seb
1Deb 10oks 19bb
1000
100ep 14b
10ke 1ocks INb
Window size
Window size
Window size
Osserved
Predicted
Mutat
par Mb
Centromere
M
wWW
kenbification TERT driver melanoma
-egg0 1 observed preeluu
50-
Gene searth: Emppraaped
IFIAA prom.
Muttatoons
1066
50-
40-
Regulatory search: Liver
TP53
WDR74
prom.
30-
1.5
2.0
SMADD
INCRNA
2.0
MAP2
KRAS
DAK Skp
TERT driver detected
1.5
1og10 expected p-value
Position chromosome IMB)
Figure 3: SPG accurately models mutation density and detects driver events, Variance explained (Pearson R²) the observed mutation count by our SPG, binomial regression (BR) and negative binomial regression (NBR) across length-scales. b. Observed (dashed blue) and predicted (solid orange) mutation density from our GM at 1Mb regions across chromosome in melanoma. (c) Quantile quaantlee plots of expected and observed P-values for gene driver detection in esophageal adenocarcinoma (purple) and driver regulatory element detection in liver cancer (pink). Esophageal and liver were chosen only for the sake of readability: qq plots are similar for all cancers. (d) Model detection of a well-known non-coding driver in the TERT promoter in melanoma at 1kb resolution. Black dashed lines: cnfferroni-correcter genome-wide significance thresholds.
Gene information was obtained from Martincorena et al. (2117) and regulatory element information from Rheinbay et al. 22220). Driver detection was performed in all four cancer cohorts.
RESULTS
5.1 ACCURACY OF REGIONAL RATE PARAMETER ESTIMATION
We first evaluated various methods' abilities infer regional rate parameters, considering both new (CNN+GP) and existing (RF, NBR, BR) methods. We assessed each method's ability to learn the expected mutation rate by directly assessing the amount of variance (Pearson R²) it explained over observed mutation counts in 50kb real data windows (fig. 2 top), and found the CNN+GP estimation method performed the best, although random forest was a close second (results were similar when estimating the mean in simulated data; fig. (2p middle). We then evaluated each method's ''s ability to capture the variance of in the simulated data, quantified as the Pearson R² to the true variance. The CNN+GP method again outperformed the others (fig. 2 bottom). Notably, RF was unable to infer the variance beyond chance level, and thus we did not consider this method further because its inability to infer variance precludes accurate driver detection.
also considered other dimensionality reduction techniques including both non-neural and neural approaches well as supervised and unsupervised approaches, other approach achieved accuracy comparable to the CNN+GP over both mean and variance (see Appendix). Moreover, both parameters and found it significantly reduced model performance (7% decrease over mutation counts and 13% over of within 10kb windows in melanoma).
alternative to the CNN;
validated the necessity of the GP by directly optimizing the CNN predict
5.2 ACCURACY AND EFFICIENCY OF MUTATION RATE PREDICTION
To further compare the SPG performance to existing methods, we evaluated the accuracy and efficiency of each method over length scales ranging over 5 orders of magnitude (10-10$ positions). To evaluate SPG, we estimated the background mutation rate parameters, M4 and o in 10kb regions genome- wide using the CNN+GP estimation strategy; we then applied the SPG distribution to estimate mutation count distributions over all other region sizes. The existing methods with
7

---

Published as a conference paper at ICLR 2021
Method
100bp
1kb
10kb
100kb
1Mb
Multi-Scale
SPG
4m Is
3m333
36m35s
19s
28
42m40s
NBR
1h30m
7m3s
43s
6s
4s
1h37m56s
BR
44m36s
3m8s
158
5s
4s
48m8s
RF
>15h
>15h
14h2m
5m24s
28s
>15h
Table I: Run times for SPG, NBR, BR, and RF for five region sizes and multi-scale search. Reported times are for single rrin-validation--est split per model over 8 CPUs and 1 GPU machine. For SPG, parameters were inferred using the CNN+GP estimation method at 10kb, running the CNN and GP one time each; hence SPG's increased run time at 10kb relative to other region sizes. Bolded is the best multi-scale search time. Presented RF times are high due to the need for O(n) trees to estimate variance using the Jackknife method.
reasonable performance on both mean and variance prediction (BR and NBR) were trained to directly predict the count distribution in each region for each length scale genome-wide.
Across all tested window sizes and cancers, SPG outperformed existing methods. with performance particularly improved in esophageal adenocarcinoma and skin melanoma (fig. Bb)), crucial for high- accuracy driver detection downstream. Across 1 Mb windows, SPG explains 95% of the variance in mutation density across all three cancers (fig. Bb,b); this is >15% more variance than both existing methods (Fig. 3h) highlighting the ability of SPG to accurately capture regional distribution parameters and project them upwards. The decrease in variance explained in smaller window sizes is expected because observed mutation counts become increasingly stochastic relative to the expected number of mutations predicted by each method. The theoretical foundations of negative binomial regression and SPG are similar, both built upon the classical Poisson Poisson-gamma gamma model. SPG differs from NBR in three key ways that help explain its improved performance: 1) SPG models mutation patterns over arbitrary sets of positions enabling it to dynamically pool information across positions after a single training; in contrast, NBR operates on fixed regions and must be retrained for every new region size. 2) SPG's s variational inference method estimates the gamma parameters for each region independently; NBR estimates only the shape parameter independently for each window and uses a single scale parameter for all windows. 3) SPG's CNN data reduction enables non- non-linear mapping of spatial covariate information to mutation rate, whereas NBR can perform only linear inference and disregards the spatial organization of the genome.
SPG is also the most efficient method for multi-resolution search (appendix D.3). Initial training of parameters using the CNN+GP method for one fold of 10kb regions required 36 minutes using GPU. Projection to each additional scale using 8 CPUs required at most 4 minutes (table []. In contrast, training time for BR and NBR increases considerably the resolution decreases. Performing a search across resolutions of 50bp, 100bp. 500bp. 1kb, and 10kb would require >5h for negative binomial, >2h for BR, and only 52 minutes for SPG (Appendix). We have also found that parameter estimation on windows as large as 100kb does not significantly reduce accuracy across scale (Appendix). allowing SPG parameter estimation in considerably shorter time (e.g. only 8 minutes for 50kb).
5.3 IDENTIFICATION OF CANCER DRIVER MUTATIONS
We leveraged SPG's ability to model multiple resolutions to search the whole genome of each of the four cancer cohorts for gene drivers, non-coding regulatory drivers, and 50bp windows that may harbor a driver mutation. All significant results are provided as supplementary data tables. We compared our results those obtained from a previous comprehensive characterization of these cohorts by [ampbell et al. (2020), who used 13 different methods to identify drivers. Our model did not have access information about gene structure or function unlike the methods used in the previous characterization. Nonetheless, the model's 's p-values were well calibrated (fig. (3)), and we identified 19 genes with significant excess of missense or nonsense mutations. All T9 genes were previously reported as drivers by Campbell et al. (2020). We failed to detect only two known driver genes present in >5% of samples. This performance is on par with sfate-of-the-an methods specifically designed for driver gene identification (Rheinbay et al. 2020).
When analyzing non-coding regulatory elements, SPG's valuee were again well calibrated (fig. Bp), and it identified all non-coding drivers (n=11) identified by [Cambell et al. (2020). Moreover, SPG implicated several additional putative non-coding driver elements that had not been previously reported. Examples include 1) the promoter of the gene MTERFDI in esophageal cancer (P 3.1 x 10 10-8), whose over expression has been observed in numerous cancers, has been shown to promote cell growth, and decrease clinical survival (Zhang et al. 2014); 2) enhancer of DHX33 in
8

---

Published as a conference paper at ICLR 2021
liver cancer (P = 4.8 10-!1), whose over-expression has been shown to promote cancer development (Wang et al. 2019); and 3) the 5' UTR of ERNI in melanoma, which has been linked to cancer therapy resistance (SuSiić et al., 2018).
Finally, we performed the first, to our knowledge, genome-wide search for individual driver mutations. All significant genic hits fell within known driver genes whose functions have been experimentally validated including 7PF3, BRAF, KRAS, PIK3CA, and CTNNBI. In addition, SPG identified two recurrent mutations in the genes GPR98 and KLB that had not been previously identified in Campbell et al. (2020'ss 's analysis of the data. These mutations are listed as driver mutations in the Catalogue of Somatic Mutations in Cancer Tate et al. 2018). SPG implicated numerous hotspots in the mostly unexplored non-coding genome, including the well-known TERT promoter mutation (fig- (t). These results are promising targets for future studies of non-coding drivers in cancer cell lines and organoids.
DISCUSSION
We introduced an extension of the Poisson-Gamma distribution to model discrete-time, integer- valued stochastic processes at multiple scales. The split-Poisso-Gaamm (SPG) model makes several simplifying assumptions including: 1) that the process quasi-stationary in a small enough region; 2) events are distributed among the discrete units approximately independently; and 3) the behavior of the random variables can be captured by particular parametric distributions. The assumptions are necessary to derive a closed-form posterior distribution. This enables efficient prediction over multiple length-scales without having to re-estimate the model parameters. We additionally proposed variational inference strategy to reduce input dimensionality and estimate the parameters of the model using a CNN coupled with a GP. Indeed, the use of a CNN+GP for a distribution variational inference may be of use well beyond the SPG framework and discrete stochastic process modeling.
To demonstrate the utility of the SPG, we applied it to model mutation rates in cancer and identify genomic elements that drive tumor emergence. In the case of this application, previous work has established the validity of the above assumptions, demonstrating that the mutation rate is approxi- mately constant within 50kb regions (Rheinbay et al.) 2020) and that mutations occur approximately independently given each position's is trinucloeotide context (Martincorena et al. (017). We demon- strated that the approach is more accurate than other methods on both real and synthetic data. We also demonstrated that multi-resolution prediction enables identification of both known and novel putative drivers of cancer, including in the non-coding genome, a crucial open problem in genomics (Khurana et al., 2016 Rheinbay et al. 2020
SPG is also applicable to discrete stochastic challenges in other domains, particularly when anomaly detection is the goal. For example, cybersecurity is often interested in detecting malicious network activity that may occur over seconds, hours, or weeks. Such activity ought to appear as anomalous relative to the expected network traffic. However, similarly to cancer drivers, detecting such anomalies confounded by the fact that expected network traffic can vary dramatically over time. Thus detecting malicious activity requires modeling non-stationary event rates and searching for anomalous activity across multiple resolutions. SPG is highly suited for efficient execution of this task. While the details of parameter estimation will depend upon the application, we expect the variational Gaussian process approach will be broadly applicable and that, in the case of high-dimensional matrix input, a CNN will provide powerful tool reduce the data to an informative feature vector.
Another timely use-case is identifying infectious disease outbreaks. The task is to determine when a new infection hotspot is developing to implement containment measures. This task is challenging because infection rates vary by geography and by individuals (e.g. young vs. elderly). SPG provides a framework to identify infection hotspots while accounting for geographic and demographic risk. In this application, regional rate parameters would reflect geographic infection rates while individual people would be the "positions" of the stochastic process. The task would be to identify groups of people (e.g. schools, neighborhoods, cities, etc.) with more cases than expected. Local information (e.g. medical treatment availability, number of cases, population density, etc.) could serve as predictors to infer regional infection rates while an individual's demographics could provide the clustering criteria to estimate pl. Such an application could be particularly useful for early identification of hotspot outbreaks of COV These examples highlight the diverse situations to which SPG could be applied, and we expect that the number of applications will continue to grow as collections of time- and spaci-varying data grow increasingly large.
9

---

Published as a conference paper at ICLR 2021
ACKNOWLEDGEMENTS
We acknowledge the contributions of the many clinical networks across ICGC and TCGA who provided samples and data to the PCAWG Consortium and the contributions of the Technical Working Group for collation, realignment, and harmonized variant calling of the cancer genomes used in this study. We thank the patients and their families for their participation in the individual ICGC and TCGA projects.
REFERENCES
N. Abdennur. Python bindings to UCSC BigWig and BigBed library. Contribute nvictus/pybbi development by creating an account GitHub, November 2018. URL htttps://iihhu
com/ http///ccttsyypypeeee original-date: 2016-05-16719:11:5585
K. C. Akdemir al. Somatic mutation distributions in cancer genomes vary with three-dimensional chromatin structure. Nature Genetics, 11:1118778-118 November 2020. ISSN 1546- 1718. doi: 10. 100-88/415880000---00
URL https: /www. htt.://ttue.eeeem////////// naturatuuee.ee wtt.coom//ra.. com/articles/
641588--20-0708-11 Number: 11 Publisher: Nature Publishing Group.
L. B. Alexandrov et al. Signatures of mutational processes in human cancer. Nature, 500(7463): 415-421, August 2013. ISSN 1476-4687. doi: 10.10388nattuee122 2477. URL https: //ww/ /www nature.
J. Bertl, Guo, M. Juul, S. Besenbacher, M. M. Nielsen, H. Hornshøj. J. S. Pedersen, and A. Hobolth. A site specific model and analysis of the neutral somatic mutation rate in whole-genome cancer data. BMC Bioinformatics, 19(1):447, April 2018. ISSN 1471-2105, doi: 10.1186/s12859-018 018-2141-2.
URL httpst:s:: peletps://doi.oom http:://dooooogg////
86 2211-889--222 8-21
J. Bradshaw, A. G. d. G. Matthews, and Z. Ghahramani. Adversarial Examples, Uncertainty, and Transfer Testing Robustness in Gaussian Process Hybrid Deep Networks. arXiv:1707.0244 [stat], July 2017. URL
707 01t./717777777666
arXiv: 1707.02476.
P. J. Campbell et al. Pan-cancer analysis of whole genomes. Nature, 578(7793):82-93, Februa 2020. ISSN 1476-4687.
doi: 10.1038//41586-0222 1969-6. URL http:::/tt//w nat ure. com/ articles articles/s41586-020-- 020- 1969-6
Number: 7793 Publisher: Nature Publishing Group.
R. G. Gallager. Stochastic Processes: Theory for Applications. Cambridge University Press, December 2013. ISBN 978-1-107-43531-5 Google-Books-ID. CGFBAgAAQBAJ.
J. R. Gardner, G. Pleiss, D. Bindel, K. Q Weinberger, and A. G. Wilson. GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration. arXiv:1809, 11165 [cs, stat], January 2019. WtlLattp://aaiiy/ http::///axxxv..r/a////////s////
111111899111111.. 11165
arXiv: 1809.11165.
B. S. Gloss and M. E. Dinger. Realizing the significance of noncoding functionality in clinical genomics. Experimental Molecular Medicine, 50(8), August 2018. ISSN 1226-3613. doi: 10.
01515/12276-118-088888 URL htt.s: www. ncbi. wwwmmnnn.....m.mm nih. hwt.m.iin....i.vo..... wtt.g//www/./cccleeeee/ PMC60828317
A. Gonzalez-Perez, Sabarinathan, and N. Lopez-Bigas. Local Determinants of the Mutational
Landscape of the Human Genome. Cell, 101(1:101-114, March 2019. ISSN 0092-8674. doi: 10.
1016// 2019 02.051. URL www. www.::sciiincddiee...... com wttp//////eeeeeeeeieeieeeee pii 679223333377
D. Guo, B. Chen, H. Zhang, and M. Zhou. Deep poisson gamma dynamical sys- tems. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems 31, pp. 8442- 8452. Curran Associates, Inc., 2018. URL http:://pppet
httipsccc/pape..
8064 -
gamna dyr.miccall seallssstemss. p
M. Juul. J. Bertl, Q. Guo, M. M. Nielsen, M. Šwitmicki, H. Hornshøj. T. Madsen, A. Hobolth, and
J. S. Pedersen. Non-coding cancer driver candidates identified with a sample- and position-specific model of the somatic mutation rate. eLife, 6. ISSN 2050-084X. doi: 10.7554/cLife.21777 URL https: / wwwww ncbi i nlm. m.m n wtt..h.www.e/.eeeeeooocooooe
/ 69/
10

---

Published as a conference paper at ICLR 2021
M. Juul, T. Madsen, Q. Guo, J. Bertl, A. Hobolth, M. Kellis, and J. S. Pedersen. nedDetect2: improved models of the site-specific mutation rate in cancer and driver detection with robust significance evaluation. Bioinformatics, 35(2):189-199, January 2019. ISSN 1367-4803. doi: 10. 10.1093biionnfommaaiccsr11111..1 URL https: //ac///accaamaaa c. www. http://///oooonoooooooo/oe att.iclee/35/2/1875545333
Publisher: Oxford Academic.
E. Khurana, Y. Fu, D. Chakravarty, F. Demichelis, M. A. Rubin, and M. Gerstein. Role of non-coding sequence variants in cancer. Nature Reviews Genetics, (7(2):33-188 February 2016. ISSN 1471-0064. doi: 10.1038/mr.20055.17 2015.17
na http:://www.ntuuurr..cc///rrrrrr
M. S. Lawrence et Mutational heterogeneity in cancer and the search for new cancer-associated genes. Nature, 499(7457):214-218 July 2013. ISSN 1476-4687. doi: 10.1038/nature1220 URL
atttps://wwx www.natturee. com htt.aaaaww.....caaaaaaaaaaaattta
J. K. Lindsey. Statistical Analysis of Stochastic Processes in Time. Cambridge University Press, August 2004. ISBN 978-1-139-45451-3 Google-Books-ID> podDRPIOFTQC.
L. Lochovsky, J. Zhang, Y. Fu, E. Khurana, and M. Gerstein. LARVA: an integrative frame-
work for large-scale analysis of recurrent variants in noncoding annotations. Nucleic Acids
Research, 43(17):8123-8144 September 2015. ISSN 0305-1048. doi: 10.10933ma//gkv803.
URL https: academic. oup. httppp..oo//a/atciceeeeeeee httpa////.////////aa//////////// lisher: Oxford Academic.
Pub-
I. Martincorena and P. J. Campbell. Somatie mutation in cancer and normal cells. Science, 349(6255): 1483-1489, September 2015. ISSN 0036-8075, 1095-9203, doi: http///cciccccaaaa00.5.... URL https:// //tc:///ccinnceeeee sciencemag org/ wtt.a/tent/34942//eeeee can Association for the Advancement of Science Section: Review.
Publisher: Ameri-
I. Martincorena, K. M. Raine, M. Gerstung, K. J. Dawson, K. Haase, P. Van Loo, H. Davies, M. R. Stratton, and P. J. Campbell. Universal Patterns of Selection in Cancer and Somatic Tissues, Cell, 17165):1029-1041. November 2017. ISSN 00928674, doi: cel.0001/..
URL
www.uub.ellviiiii com http://w.....................///
L. Mclnnes, J. Healy, N. Saul, and L. Grossberger. Umap: Uniform manifold approximation and projection. The Journal of Open Source Software, 3(29):861, 2018
L. Mularoni, R. Sabarinathan, J. Deu-Pons, A. Gonzalez-Perez, and N. López-Bigas. OncodriveFML: general framework to identify coding and non-coding regions with cancer driver mutations.
Genome Biology, 17(1):128, June 2016. ISSN 1474-760X. doi: 10. 10.1186/11059--11--099
S. Nik-Zainal et al. Landscape of somatic mutations in 560 breast cancer whole-genome sequences. doi: 10. 10.1008//ature17661 URL Number: 7605 Publisher: Na-
Nature, 534(7605):47-54, June 2016. ISSN 1476-4687.
https: www. www.uatuue... htt.://eccoooerrt...es////ttttte
ture Publishing Group.
J. Paciorek and M. J. Schervish. Nonstationary Covariance Functions for Gaussian Process Regres- sion. In Thrun, L. K. Saul, and B. Schölkopf (eds.), Advances in Neural Information Processing Systems /6, pp. 273-280. MIT Press, 2004. URL. http::ttt///. tttttp://papprrs.. wwwp.cc.ppppp......eeppppppppppe
2350 1rss0=oonntatinnary--oooo wtw.aaiiooncccyvaaaaaaaaa pdf
functions for for- oaaaee-iiiiiiiio-oooooooaaaaaaae www.eeseeeeeeeesiiooo.o..o
A. Paszke, S. Gross, S. Chintala, G. Chanan, Yang, VeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer. Automatic differentiation in pytorch. 2017.
Polak al. Cell-of-origin chromatin organization shapes the mutational landscape of cancer. Nature. 518(7539):300-366 February 2015. ISSN 1476-4687. doi: 10.1038/nature 14221. URL https: www. nature. com www./atticee//eaeeeeeeooeoeeeeee wtw.:/iiiilee/nauuee22222
K. Polimis, A. Rokem, and B. Hazelton. Confidence intervals for random forests in python. Journal of Open Source Software, 2(1), 2017.
11

---

Published as a conference paper at ICLR 2021
E. Rheinbay al. Analyses of non-coding somatic drivers in 2,658 cancer whole genomes. Nature, 178(7793):102-111 February 2020. ISSN 1476-4687, doi: 10. 10.10338/41188-022 1965-x. URL
https: www nature. com httpe///me.aaaeeee//1115------00
M. D. Risser. Review: Nonstationary Spatial Modeling, with Emphasis on Process Convolution and Dovariate-Drivens Approaches. arXiv:1610.m 02447 [stat], October 2016. URL htttp://arxiv
org/abs/1610 02447
arXiv: 1610.02447.
Roadmap Epigenomics Consortium et al. Integrative analysis of 111 reference human epigenomes. Nature, 518(7539):317334) February 2015, ISSN 1476-4687, doi: 10. 1038/atree 14248. URL
https: www nature www.eurree....me.. http://arccceseettteeeeeee
A. Schein, H. Wallach, and M. Zhou. Poisson-ganas gamma dynamical systems. D. Lee, M. Sugiyama, U. V. Luxburg. Guyon, and R. Garnett (eds.), Advances in Neural Information Processing Systems 29, pp. 5005-5013. Curran Associates, Inc., 2016. URL | http:
http:://ppprrs...// ni nips ccpeaper/
6083 5
gamma -amma-ddnammaaaa..mm.mmaamm ca systems pdf
B. Schuster-Böckler and Lehner. Chromatin organization is major influence on regional mutation
rates in human cancer cells. Nature, 488(7412):504 507, August 2012. ISSN 1476-4687. doi: 10. 10.1038/nature1 1273. URL http:tt::/// //ww.. nature com htt.:///rtt..ceee/aateeeeeeel//
S. Seabold and Perktold, statsmodels: Econometric and statistical modeling with python. In 9th Python in Science Conference, 2010.
M. R. Stratton, J. Campbell, and P. A. Futreal. The cancer genome. Nature, 458(7239):719-724 April 2009. ISSN 1476-4687. doi: 10. 1038/nature07943. URL (tLLhttps::// www. nature..ree....... com/ artic het.eeeeeattrrr/4444
J. G. Tate al. COSMIC: the Catalogue Of Somatic Mutations Cancer. Nucleic Acids Research, 47(D1):D941-D947, 10 2018. ISSN 0305- 1048. doi: 10. 10.1093/nar/gky1015 1015.
URL https:
///ba..org/10.109//rr//////
M. K. Titsias. Variational Learning of Inducing Variables in Sparse Gaussian Processes. pp.
L. Wadi Candidate cancer driver mutations in superenhancers and long-range chromatin interaction networks. bioRxiv, pp. 236802, December 2017. doi: 10.1101/236802. URLhttps Publisher: Cold Spring Harbor
w/www. www.:iii.iiv..... org/
Laboratory Section: New Results.
S. Wager, T. Hastie, and B. Efron. Confidence Intervals for Random Forests: The Jackknife and the Infinitesimal Jackknife. pp. 27.
J. Wang, W. Feng, 7. Yuan, J. D. Weber, and Y. Zhang. Dhx33 interacts with ap-2ß to regulate bel-2 gene expression and promote cancer cell survival. Molecular and Cellular Biology, 39(77), 2019. ISSN 0270-7306. doi: N0.112//CB.000177
BRLhttts:: //cc/////.///////////////cc../at asm. oog.//ntett./t./
39/17/e00017-D
D. Weghorn and S. Sunyaev. Bayesian inference negative and positive selection in human cancers. Nature Genetics, 49(12): December 2017. ISSN 1061-4036, 1546-1718. doi: 10.1038/ng.3987. URL www.ttttp:// www nature. caammaarticom wtt.:/re.....oaart.l/s/////////e
C. Zhang et al. Mterfdl functions as an oncogene. Oncotarget, 500), 2014. ISSN 1949-2553. doi: http:///....og/0.0066662222ooeae article/111er
140. URL https: wwwp///www/////////////////wo www..oncotarge..... com/
Šuštić, S. van Wageningen, E. Bosdriesz, R. J. D. Reid, J. Dittmar, C. Lieftink, R. L. Beijersbergen, L. F. A. Wessels, R. Rothstein, and R. Bernards. A role for the unfolded protein response stress sensor ERNI in regulating the response to MEK inhibitors KRAS mutant colon cancers. Genome Medicine, 10(1):90. November 2018 ISSN 1756 994X. doi: 10. 10.11866 1011188//13073000 011-00000-- URL
12

---

Published as a conference paper at ICLR 2021
A APPENDIX
In this appendix, we provide detailed information on:
1. The data used in this work including its origin and all preprocessing steps.
2. Additional method details including:
A derivation the closed-form marginal distribution of the graphical model presented in main text.
Architecture and training details of all models.
How the genome-wide search for driver mutations was performed.
The appendix also includes analysis of the sensitivity of negative binomial regression to detect well-known drivers genome-wide wiide and additional figures that provide context to results presented in the main paper.
B DATA
B. EPIGENETIC TRACKS
We obtained 733 1ooguo(Pvvalue) chromatin tracks representing the epigenetic organization of 111 human tissues from Roadmap Epigenomies Roadmap Epigenomics Consortium al. (2015) (see Appendix table "predictor_ ppredictr_ttackk.... seteee_etiiiosscsv//// ccs"". These tracks measure the abundance of particular chromatin mark genome-wide, with smaller (more significant) p-values reflecting a greater abundance of the chromatin mark genomic position. Chromatin marks chemical modifications of histones, the proteins used to package DNA within a cell. We additionally obtained 10 replication timing tracks from the ENCODE consortium. Replication timing assays measure the relative time which each position in the genome is replicated during cell division. For non-overlapping regions of predefined size and location (see main text for more details), we extracted the signal for each epigenetic track using 100 bins per region with pybbi Abdennur (2018). We additionally calculated the average nucleotide content each window by assigning each nucleotide a numeric value between I and 5 and taking the average across bin (N [unspecified nucleotide] I, A = 2, C G = 4, T 5), and we calculated the GC content as the percent of and C nucleotides in a bin, resulting in total of 735 epigenome tracks region. The mean values for each region were calculated as the mean chromatin signal for each track in the region.
B.2 MUTATION COUNT DATA
We downloaded somatic single-base substitution mutations identified in the ICGC subset of the Pan-Cancer Analysis of Whole Genomes Consortium cohorts of esophageal adenocarcinoma, skin melanoma, stomach adenocarcinoma, and liver hepatocellular carcinoma. These data are freely available for download from the International Cancer Genomics Consortium data portal (see fig. 4 We excluded mutations on sex chromosomes (X Y) because males and females carry different sets of these chromosomes, leading to differential mutation patterns. We summarized the data as mutation counts per window for window sizes of 50bp. 100bp. 500bp, 1kb, 5kb, 10kb, 25kb, 50kb, 100kb, and 1Mb.
B.3 RESTRICTION TO REGIONS OF HIGH MAPPARILITY
High-throughput genome sequencing works by randomly reading millions of short sequences of nucleotides 156-150 bases in length) from a target genome. These "reads" are then mapped to the human reference genome to reconstruct the target. A challenge is that short sequences of k nucleotides (kmers) can occur multiple times in the genome. This results in ambiguous mappings for some reads and thus degradation of data quality in regions composed of many kmers that occur multiple times across the genome. Following previous work Polak et al. (0015), we removed regions the genome with low quality data by calculating mappability score for each region. Mappability scores reflect how many times a particular kmer occurs in the genome and have been pre-computed for the human reference genome GRCh37. We required that region's 's average mappability score based 36mers
13

---

Published as a conference paper at ICLR 2021
(e.g. average across all sequences of 36 nucleotides in the region) be >70% reflecting that all 36mers in the region be >70% unique. The majority of the genome passed this threshold; for 10kb regions, for example, >75% of the genome passed this threshold. We chose to measure mappability with 36mers because this was the length of read used to generate the Roadmap Epigenomics sequencing data.
Esophageal Adenocarcinoma
Melanoma
Stomach Adenocarcinoma
8000
6000
5000
6000
Counts
4000
4000
3000
2000
1000
1000
0
100
200
300
200
400
600
800
20
Mutations/ / /5kbb
Mutations/ / 50kb
Mutations / 50kbb
Figure 4: Distribution of mutation counts in 50kb windows tiled across the genome with 36mer uniqueness 770% (see section B.3) for esophageal adenocarcinoma, skin melanoma, and stomach adenocarcinoma. Of note: esophageal adenocarcinoma has highly skewed distribution, skin melanoma has high mutation counts relative the other cancers, and stomach adenocarcinoma has low mutation counts relative to other cancers.
B.4 SYNTHETIC DATA SIMULATION
We generated synthetic datasets for each of the cancers in order to have datasets with known mean and variance rate parameters. To generate the datasets, we used a senearest-neighhee strategy identify the 500 nearest neighbors for each region. The mean and variance for that region were then taken to be the empirical mean and variance calculated from the 500 nearest neighbors. The number of "observed" mutations was then randomly sampled from a binomial defined by the mean and variance parameters. It is important to note that these datasets are purely derived for the purpose of comparing methods over datasets with a known ground-truth. They do not reflect mutation patterns in the real datasets. The specific steps to generate the simulated data were:
1. Generate vectors of the values for each of the 735 tracks (733 epigenetic tracks, GC content track, and average nucleotide content track) in 50kb regions of the genome with 36mer uniqueness >70%.
2. Perform ordinary least-squares (OLS) regression of the mean vectors against the observed number of mutations in each 50kb window for that cancer.
3. Scale each value in the feature vectors by its corresponding coefficient from OLS and compress the weighted mean vectors to 50 components using Principal Components Analysis (capturing >94% of the variance for each cancer).
4. For each region R, perform k-nearest-neighbor clustering with Euclidean distance to identify its 500 nearest neighbors in the PC space. Define the mean μı and variance o of the mutation rate in R to be the mean and variance of the KNN cluster.
5. For region R, randomly draw new "observed" number of mutations from negative binomial distribution defined using the associated mean and variance. Specifically, X ~ NB(a, 10 1)) where /2/ o and 0 = of/ /uR
We created two versions of the simulated data, in which all regions in the genome were used to estimate the rate parameters and one in which rate parameters were estimated separately within independent train and test subsets, Results were qualitatively indistinguishable.
14

---

Published as a conference paper at ICLR 2021
C METHODS
C.I GRAPHICAL MODEL
Here we derive the closed form negative binomial distribution presented the main text as the graphical model marginal distribution events at some unit i in region R. We use the following notation:
Mj # mutations observed at pos (observed)
pt genome-wide probability of observing mutation at the nucleotide context of i (inferred)
iy normalized probability of observing mutation at i in region R (inferred)
2p the background mutation rate in region R (unobserved)
Xg: # background mutations in region R (unobserved)
Mgg the expected background mutation rate in region R (inferred)
o². the variance of background mutation rate in region R (inferred).
Me covariates associated with the behavior of the stochastic process within R (observed)
As presented in the main text and main Figure 1, the graphical model implies the factorization
Pr,M, ARl ag. OR. Pi; TR) Pr(M k|XR. ii MR) Pr(Xx x ARi 7R) (ie)-Pr(Ag) Og., 0g: nR) (10)
where
aR μ o
0₂ o HR.
Since πR is given in each equation, we suppress it for notational ease. To marginalize out XR, note that
Pr(M: = k(2R) Σ Pr(M, k[XR pi) Pr(Xx = x(2R)
is equivalent to split Poisson process GGallager 2013). Thus
Pr(M: = k Aa)
Possion (M; = k: p.
(11)
We now marginalize out the unknown rate parameter 2qp-
1 P(M: k|ARR: P caR OR d22x
PMM: = klp ORR, OR)
1/a0 (piAa)
BiAR
A AR. dÂR
k!
[(aRR) 0
2
- 1/ dÂR.
kII (LR)
0
Making the substitution , 2,, 1/Op) and noting that the resulting integrand is unnormalized gamma distribution, we have:
Op
P(M = klp;, CRR 0R)
I(cRR
kII OLR 0TH
11 + 1/ OR
I(aR +k)
piOR
RR
(.T((RR)
5.0p +1
5² +
NB
M, = k: CRR
p,0p +
15

---

Published as a conference paper at ICLR 2021
C.2 OVERVIEW OF PARAMETER ESTIMATION PROCEDURE
Estimation of regional rate parameters: As training data, we use a set of imput matrices {πR; R E } and associated mutation counts (XR; R J.. First, a CNN is trained to take ΠдR as imput and predict XR as output, using mean squared error loss. The final 16-dimension feature vector of the trained CNN then used as input to train Gaussian process to predict the mutation count XR and the associated estimation uncertainty by maximizing the likelihood of the observed data. The mean and variance output by the GP were used as estimates for μRR and oz.
Estimation of time-averaged event probabilities: the time-average probability of an event at Pi was estimated based on it's trinucleotide composition, t.', where n is the nucleotide at i 1- 1, is the nucleotide at and n' is the nucleotide at i 1 in the reference genome. We first counted every occurrence of n,t. n in the human genome and then counted the number of times the middle nucleotide of the 3mer was mutated across the genome. The maximum likelihood estimate of pi is then the ratio of the number of observed mutations of the 3mer divided by the total occurrences of the 3mer.
C.3 REGIONAL PARAMETERS ESTIMATION METHODS
To compute a model's R² accuracy to UR and of for regions R of size S, the genome was divided into non-overlapping contiguous segments of size S. To assure high data quality, any region with mappability score < 70% was excluded from further analysis. The remaining windows (accounting for more than 75% of the genome) were randomly divided into train and test sets in 80-20 split respectively. The test set was held-out and served solely for evaluation purposes. The train set was then divided into train and validation sets by another 80-20 split respectively (train set 64%, validation = 16%, and test = 20% of the considered regions with mappability score < 70%, see appendix B.3
C.3.1 GAUSSIAN PROCESS FEATURE VECTOR GENERATION
All networks were independently trained for 20 epochs with a batch size of 128 samples and using the Adam optimizer to minimize mean squared error loss to either the true mutation count (CNN and FCNN) or input tensor (AF). After training the model parameters using the train set, predictions over the held-out test set were computed by 1) extracting the last 16-dimensional feature layer (middle feature layer for AE) for all sets over the best performing model over the validation set across all epochs (according to the validation accuracy); 2) training multiple GPs (typically 10) to predict mutation counts using the 16 dimension feature vectors of the train set as input (see appendix C.3.2 for details); 3) taking the mean UR and of of all 10 runs over the test set as the ensemble prediction of the model. All neural network models were implemented in Pytorch Paszke et al. (2017).
1. Convolutional neural network (CNN): The CNN contains 4 convolutional blocks with 2 batch normalized convolutional layers and ReLU activation. The first block transformed the input tesor from 735 x 100 to 256 x 50 with 256 channels and a double stride. The other blocks are ResNet-style residual blocks that maintain their input dimension to facilitate residual connections, with 256, 512, and 1024 channels respectively. Between each of the 3 residual blocks there is a double stride (ReLU activated and batch normalized) convolutional layer, which divides the tensor length by two and doubles its height with additional channels. The output of the last residual block is flattened and passed through fully- fully-connected layers. The first two are ReLU activated and reduce the dimensionality of the tensor to 128 and 16 dimensions respectively. The last uses linear functions to reduce the tensor to a single cell holding the output of the regression. This forces linear relation between the regression output and the last feature layer, thus simplifying the function the GP needs to learn, which we found empirically improves the GP's accuracy.
2. Fully- neural network (FCNN): The FCNN has an architecture similar to the CNN's 3 fully-connected layers but with an input space of the mean epigenetic vector (735 dimensions). Thus, the FCNN is computationally similar to the CNN, but operates on the mean vector instead of the full matrix an input. The FCNN is designed to demonstrate maximum performance possible when reducing the input tensor to an averaged feature vector.
16

---

Published as a conference paper at ICLR 2021
3. Autoencoder neural network (AE): The encoder of the AE used the same architecture as the CNN, excluding the last linear fully connected layer. The decoder has a mirror architecture with the same number of parameters but differs in the internal design of the convolutional blocks. Convolutional layers were replaced by 1D transpose convolutional layers with no batch normalization and residual connections. The AE was designed to demonstrate the predictive power of feature embedding that was not optimized to specific task but produced in a way comparable the CNN.
4. Other dimensionality reduction methods: PCA was computed using the Python Scikit- learn package with default settings and UMAP was computed via Python's 's umap-learn package Mclnnes et al. (2018) with 20 nearest neighbours and Euclidean distance. Both methods were computed over the entire training set (80%) with no validation set and reduced the mean epigenetic vector dimensionality (735 dimensions) to 16, just like all other models. Prior to processing, we log- ongtrranforreed the epigenetic data as we found this improved prediction accuracy downstream.
C.3.2 GAUSSIAN PROCESS
We implemented a sparse, inducing-poinn Gaussian process Titsias with a radial basis function kernel using Python's 's GPy Teychrch package mackaee/Gardner Gardner al. (2019). The GP was optimized with 2000 inducing points using the Adam optimizer for 100 steps. All features were mean-centered and standardized to unit variance prior to training. For each dataset, we ran the GP ten independent times and calculated the ensemble mean of the mean and variance predictions from each of the individual runs. We took these ensemble predictions as the mean and variance for each region.
C.3.3 ALTERNATIVE MODELS
We implemented previously proposed alternative methods Polak et al. (2015); Nik- et al. (2016); Martincorena al. (2017) for the estimation of μR and oğ without the use of GP. These methods use the mean epigenetic vector as an input.
1. Random forest (RF); RF regression was implemented via the Ensemble Methods module in the Python Scikit-learn package, with maximum depth of 50 trees. Since RF does not directly compute variance, we implemented the Jackknife method as described in Wager al. (we have compared our implementation Polimis et al. (2017) and found them highly correlated). Wager et al. suggests that the number of estimators, i.e., trees, must be linearly related to the number of samples obtain reasonable estimates of the variance. We chose to have one tenth as many estimators as samples in an attempt to keep running time within reasonable limit for datasets of smaller region sizes. Even so, for 10kb regions (containing approximately 300K regions), RF required >2 hours to train.
2. Negative binomial regression (NBR): : As described in section 3.3.2 of the main text, NBR directly specifics the variance as o = μ (1 BURR), where B is an overdispersion parameter. When B 0 NBR reduces to Poisson regression, also widely used in the community. NBR was implemented via the discrete module in the Python statsmodels package [Sabold & Perktold (2010) with the Broyden-Fletcher-
Fletcher Glodcher-Collfarb Goldearbbllabbbirooee Shanno optimization algorithm and IK maximum iterations. Epigenetic predictors were fog-transformed and reduced to 20 principle components, following the Midi-staadrrdMMaatiocooooo
al. (2017) in both train and test sets. When used to compare against the GM we also included the expected number of mutations based on the sequence context model (see main paper section 3.2) as an exposure term in the model as in previous work Nik- et al. (2016); Martincorena et al. (007)).
3. Binomial regression (BR): Following a previous study Berti et al. (2018) that suggested multinomial regression to model multiple types of mutations, we also considered binomial regression (as the binary version of multinomial regression applicable to our simple counts data) as method to model mutation rates at high resolution. BR was implemented via the generalized linear module in the Python statsmodels package Seabold & Perktold (2010). As in previous work Nik-Zainal et al. (016);; Martincorena et al. (2017), we included the expected number of mutations based on the sequence context model (see main paper section 3.2) as an exposure term in the model. As with NBR, the epigenetic predictors
17

---

Published as a conference paper at ICLR 2021
were log- transsransfor and reduced to 20 principle components for both train and test sets following tatte-of-theaat recommendations Martincorena et al. (2017).
C.4 EMPIRICAL VARIANCE ESTIMATION
For real data, the true variance in mutation counts of a region unknown. Thus to estimate variance empirically for given model, we used the following approach:
For a region in the test set, perform k-nearest neighbors clustering with Euclidean distance to identify the 500 regions in the train set that are most similar to the region of interest based on the model's feature embedding. For all models, a feature embedding of 16 dimensions was used.
2. Calculate the empirical variance as the variance of the KNN cluster.
Since feature embeddings are model-specific, we calculated an empirical variance estimate per model. The feature-vecte embeddings for models specified in section C.3.1 were the feature vectors used as input to the GP. Models specified section C.3.3 do not create or require comparable feature vectors and therefore were not considered in the main paper results. However, to measure the ability of these methods to estimate empirical variance (Fig- 7 we computed their feature vectors by 1) taking the dot product of the model parameters and the input data mean vectors and 2) reduced these scaled vectors to 16 dimensions via a PCA reduction (cxplaining 80%-95% of the variance across the different region scales). For RF, we took the model parameters to be the feature importance weights derived from the trained forest and for NBR, we used the model coefficients as the parameters.
C.5 PERFORMING A GENOME- WIDE SEARCH FOR CANCER DRIVER MUTATIONS
For each cancer, the background mutation rate parameters were estimated across the genome using 5-fold cross validation in 10kb, 25kb and 50kb regions. While the model is robust to choice of 10kb, 25kb or 50kb region size (fig. 5). the 25kb and 50kb models include some additional regions of the genome due to the the mappability threshold (see section B.3). To analyze the largest possible subset of the genome, we performed our analysis iteratively: we first searched for drivers using regions accessible via the 10kb model; then searched additional regions not accessible by the 10kb model in the 25kb model and then in the 50kb model. To search for drivers, we applied our probabilistic model estimate the mutation count distributions in 50bp regions across the genome, and we then searched for 50bp regions with significantly more observed mutations than expected under the null distribution of our model. We controlled false-discovery rate at the 0.05 level using Bonferroni conferroni-correteed corrected p-value threshold of P<le-9.
To compare our hits with known cancer drivers. we tabulated the recurrent driver mutations reported by PCAWG that were present in our dataset, including in the TERT promoter, a well known non- coding driver. While most recurrent driver mutations are activating mutations (e.g. cause a gain of cellular function), we also found recurrent mutations in the tumor suppressor genes TP53 and SMAD4. Recurrent mutations a single position are far less likely in tumor suppressor genes because any deleterious mutation can act as potential cancer-causing mutation. For example, TP53 had 6 genome-widd -ide significant 50bp regions, consistent with its status as a crucial tumor suppressor that can be knocked-out with many different mutations (see table 2). Methods specialized to discover driver genes are necessary to find tumor suppressor genes in general Lawrence et al. (2013): Mularoni et al. (0066) Martincorena et al. (007)
18

---

Published as a conference paper at ICLR 2021
Esophageal Adenocarcinoma
Melanoma
Stormach Adenocarcinoma
1.0
1.0
10kb
10kb
Variance explained
0.8
25kb
Variance explained
0.8
Variance explained
0.8
25kb
50kb
50kb
0.6
100kb
0.6
0.6
100kb
0.4
0.4
10kb
0.4
25kb
0.2
0.2-
50kb
0.2
100kb
0.0
0.0
100bp
1kp
10kb 100kb
1Mb
1006p 1kp
10kb
100kb
1Mb
100bp
1kp 10kb 100kb 1Mb
Window size
Window size
Window size
Figure 5: Model robustness region size. We tested the robustness of our GM estimates to the choicc of the scale of region R over which Hy and og were inferred with the CNN+GP. Here show GM's Pearson R² accuracy the observed number mutations over a range sizes different choices of initial region size Melanoma shows a slight decrease in performance at larger scales, suggesting local chromatin structure more strongly influences mutation rates in this cancer.
C.6 ENVIRONMENT AND COMPUTE TIME
benchmark rum 10kb scale with 10 GP reruns takes 2-3 hours on a single 24 Gb Nvidia RTX GPU, with 8 CPU cores and 756GB RAM. Thus, a full 5-fold of the entire genome takes 10-15 hours. Due to the model's robustness to scale, this time may significantly reduced without drastic loss of accuracy by using larger region scales (e.g. only 30-40 minutes for 50Kb regions, fig. 5 Importantly, after completing the CNN+GP training, projections to lower or higher scales via the GM require no additional training.
D APPENDIX RESULTS
D.1 NEGATIVE BINOMIAL REGRESSION DOES NOT DETECT WELL- -NWWN DRIVERS GENOME- WIDE
Negative binomial regression is the only other method that has been used perform unbiased genome-wide search for driver mutations Nik-Zainal et al. (2016): [Rheinbay evaluated how the sensitivity of NBR to detect driver mutations genome-wide compared with the sensitivity of our method. While all known melanoma drivers present in >3 samples were found by the GM by projecting down to only 1kb scale, NBR at 1kb fails to detect TERT, the only known common non-coding driver mutation, yielding a p- value that was an order of magnitude less significant than the genome-wide significance for this scale. Similarly, while the GM detects all known esophageal adenocarcinoma drivers by projecting down to 100bp. NBR over 100bp fails to detect KRAS, an important genic driver of esophageal cancers, again yielding a p-value that was an order of magnitude less significant than the genome-widd significance threshold for 100bp. Note: presented results at 50bp in the text to highlight our model's 's ability to search in arbitrarily small regions, but all known drivers for esophageal adenocarcinoma also detected in search over regions of 100bp.
(Rhheiinnaa et al. (2020). We thus
D.2 CONVOLUTIONAL NEURAL NETWORK OUTPERFORMS OTHER DIMENSIONALITY REDUCTION ALTERNATIVES FOR GAUSSIAN PROCESS
We first evaluated the methods for regional rate first and second moment inference, HR and o using simulated datasets. We calculated accuracy as the Pearson R² of the estimated mean and variance to the simulated ground-truth mean and variance. CNN+GP, FCNN+GP, NBR and RF accurately inferred μg, with Ryy 0.95 for all three datasets (fig. 64). However, PCA+GP, UMAP+GP, and AE+GP consistently under- dedro-perfored (fig. (6 left), suggesting supervision when creating feature vectors is critical for the GP downstream performance.
The CNN+GP and FCNN+GP outperformed the other models when estimating the simulated variance (fig. 66. right), suggesting the ability to represent arbitrary functions is important for learning
19

---

Published as a conference paper at ICLR 2021
true mean (50kb)
true variance (50kh)
Empirical variance
0.8-
Esophageal
0.6
0.4-
0.0
1,0
0.8
Neloroma R?
0.6
0.4
0.2
0.0
1.0
0.8
Stonsm
0.6-
0.0
65
KBR
PCA+oP
UNA op
AEvap
NBR
PCA+GP
LMAP-00
AE+0P
FOWR-GP Sp
CWN op
10c
2508
500D
7500
1000
10kG
73KB
100G
ROW
OWN
GP
UMAP-GP
AE+
FC
ONN-GP
UMAP
Figure 6: and of estimation accuracy three cancer types: esophageal adenocarcinoma (top row), skin melanoma (middle row). and stomach adenocarcinoma (bottom row). a. R accuracy of all models with respect to simulated He (left) and dg (right) 50kb. b. R² accuracy of GP-based models to observed number of mutations (left) and empirical variance (right) across scales in real data.
uncertainty in complex dataset. This conclusion is strengthened by the observation that UMAP and AE enabled relatively accurate variance estimation despite mediocre performance over the mean. Importantly, the clusters used for the simulated data were computed from mean epigenetic vectors; thus our CNN architecture (receiving an input in matrix form) was at disadvantage. Nonetheless, the CNN+GP most accurately learned both μą and of across all three simulated datasets (Fig 64). with slight improvement over the FCNN+GP.
To further compare the approaches, we applied the GP coupled models to estimate real mutation counts from the three cancers on multiple scales. Models were compared by their R² to the observed mutations over the test set and to an empirical variance based on the model's own feature vectors (fig. 6)) (see Appendix). The CNN+GP outperformed the FCNN+GP model over observed mutation counts and empirical variance estimation for all three cancer types. Additionally, the performance advantage of the CNN appeared to grow window size and observed mutation counts increased. This suggests that local epigenetic patterns play an appreciable role in setting mutational processes and indicates that model is well-designed to leverage the recent growth in genomics corpus sizes.
20

---

Published as a conference paper at ICLR 2021
Melanoma
Esophageal
Stomach
Observed #mutations R²
1.0
0.8
0.6
0.4
0.2
0.0
1.0
Empirical variance R2
0.8
0.6
0.4
0.2
0.0
100
25
50
100
50
75
100
cnn
random, forest
negative_ bitivmeabnnoialll
Figure 7: NBR and RF R² accuracy (with CNN+GP reference) to observed number of mutations (top row) and empirical variance (bottom row) in real data and across multiple scales for each cancer type: melanoma (left), esophageal adenocareinoma (middle) stomach adenocareinoma (right). Due to Jackknife method requirement number estimators be linear with respect to number of samples, estimating RF variance at scale <50kb was computationally infeasible (with >8,000 estimators).
Melanoma
Esophageal
Stomach
250
Predicted varience
8000
0000
20000
150-
25000-
105-
10000
100
200
700
100
100
350
200
13
Preciited mean
Predicted mean
Predicted mean
Figure Mean (UR) YS variance (o for the ground-truth simulated data (blue) and predictions each model across cancer types: melanoma (left), esophageal adenocarcinoma (middle), stomach adenocareinoma (right). NBR significantly over estimates of in high mutation count regions hecause of its strict quadratic relation the predicted mean. RF consistently underestimates of FCNN+GP accurate low to medium mutation count windows, but overestimates of with respect the CNN+GP high mutation count regions.
21

---

Published as a conference paper at ICLR 2021
Esophageal Adenacarcinoma
Melanoma
Stomach Adenocarcinoma
Expected mutations 50kk
150-
30-
300
100
200
100
200
# obs. mutations
# obs. mutations
abs. mutations
30
40-
Expected mutations 1kkb
8-
30
20-
6-
154
20-
10-
0-
0-
0
20 # mutations
10
30
mutations
Figure 9: Number of observed mutations number of expected mutations based the sequence context model alone in 50kb and 1kb regions with mappability >70% across three cancers. Sequence context explains <10% of variance 50kb and <1 % of variance kk scales for cancers.
D.3 EXISTING WHOLE-GENOME REGRESSION MODELS ARE TIME INEFFICIENT AT MULTI-RESOLUTION SEARCH
All existing regression models (RF, NBR, BR) require retraining for each desired scale. A requirement that becomes computationally challenging at finer resolutions (c.g. >1.5h for NBR at 100bp). To provide an estimate of the differences between existing methods and our SPG, we performed multi-scale time analysis presented in However, it does not include scales <100bp, such as 50bp used this work to detect driver hot-spots. log-log transform the scale against the run-time exposes a polynomial relation between the the window size and time (for small enough scales where the compute power is not governed by the machine's 's memory and system operations). Extending this relation to scale small as 50bp run-time is high 1.5h for BR and 2.5h for NBR. Making the overall run-time for a typical multi-resol sean of 50bp. 100bp, 500bp. 1kb, 10kb over 2h for BR and over 4h for NBR, while the SPG run-time remains under 1h.
22

---

Published as a conference paper at ICLR 2021
Log-log models runtime across scale
NBR
10a
Time (s)
10a
10
111
102
10a
10
10
10
Scale (BP)
Figure 10: Log-log run-time of current whole-genome regression methods that require retraining per desired scale. Run-time increases polynomially with scale beyond a threshold (where memory operations dominate the computation of the method).
Melanoma
Esophogeal
Stomach
500p Regions
50bp Regions
500p Regions
Observed logio Droalum
Conerved segi0 prestues
Oboerved logi0 pviiuves
Theoretical 10g10 p-yluess
Theoretical log10 p-values
Theoretical log10 paraues
14t Regions
11b Regions
1kb Regions
Osserved log10 povitues
Observed 1og values
1og10
20.0-
7.5
5.0
+
2.5
Theoretical log10
Theoretical log10
Theoretical log10
Figure 11: -logio(P-value) quantile-quantile (q9) plots for expected observed number of mutations in 50bp and 1kb windows using our graphical model with parameters estimated 10kb regions for each cancer. Under properly calibrated null model. p-values generated from null distribution are uniformly distributed between zero and one. QQ plots thus provide qualitative assessment of the accuracy of model's null distribution: the observed p-values should closely match the expected p-values from uniform distribution (red line) except at extremely small valuuu where observations from the alternate model should be found. The step-like nature of the qq-plot for stomach adenocarcinoma in 50bp regions is because the null distribution is discrete (negative binomial) and the dataset has relatively few mutations; thus cach 50bp bin can have only one of few possible mutation counts (typically between and 5).
23

---

Published as a conference paper at ICLR 2021
Chrom
Start
End
Observed
Expected
p-value
17
7577500
7577549
13
0,0141
1.75 x 10-30 30
17
7577100
7577149
10
0.0153
7.43 x 10-23 23
17
7577550
7577599
8
0.00856
3.72
10
17
7578400
7578449
8
0.0147
2.76 x 10 18
17
7578500
7578549
6
0.0129
6.16 x 10 14
17
7578200
7578249
6
0.0146
1.28 x 10 13
Table All 50bp windows significant recurrent mulations in the TP53 gene from genome wide driver search in esophageal adenocarcinoma.
24


Reviews:
[{"review_id": "ozE6-kc6zY0", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary: The paper extends Poisson-Gamma models for non-stationary sequences, in a manner that allows partitioning the counts according to a binomial model to account for multiple resolutions. This generalisation is motivated well with a biological application of practical relevance, and the proposed method is particularly strong in enabling linear computational scaling required for analysis of large genome data.\n\nReasons for score: I am leaning towards rejection in the current form.  The contribution itself is worth publishing and the method is likely to be valuable for the application, but the presentation would need to be improved (especially regarding the GP+CNN part; see the detailed feedback) to better communicate the technical contributions for the ICLR audience. Conditional on improved presentation, I would be leaning towards acceptance.\n\nDetailed feedback: The proposed split Poisson-Gamma (SPG) process seems reasonable, but in technical terms is a fairly straightforward hierarchical structure and can be constructed using standard properties of Poisson-Gamma. It is well motivated by the resulting efficient inference and this specific application, and may find uses in other applications as well, but does not provide a very clear theoretical contribution that would open immediate follow-up research directions for more general modelling questions.\n\nMy main problem with the paper concerns the structure of the presentation. While the authors motivate the model well and provide very clear illustrations for the application, the method sections are disconnected and I had trouble following the technical contributions. In particular, the connection between Section 2 (SPG) and the technical algorithm required for using it (Section 3.3) is unclear. To me it seems the GP+CNN part is an integral part of the overall solution and a contribution in itself (and, in fact, an important one -- the SPG alone is not quite sufficient as theoretical contribution). It provides a concrete algorithm using SPG and is general, but now the description is provided only after talking about specific data and looks more like a minor technical detail with no proper theoretical justification. For me, the paper would be more natural if Section 3.3 (and possibly some other parts of Section 3) would be described after Section 2 as description of how SPG is used in practice, and if they would use shared notation and terminology. This would make the practical approach easier to follow and the contributions more clear.\n\nThe empirical experiments and illustrations for the application are well carried out, and serve as good demonstration of the method. However, a reader uninterested or uneducated in this specific application will have some trouble figuring out how well the method works; this could be improved by complementing the results with clear artificial data of slightly simpler nature.\n\nModifications after discussion:\nIncreased score by one as the presentation in the revised version has clearly improved, along the lines requested in the original review.", "extra_scores": {}, "extra_text": {}}, {"review_id": "AXEOYOlRZvm", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors present the split Poisson Gamma (SPG) distribution, an extension of the Poisson-Gamma distribution, to model a discrete non-stationary stochastic process. SPG has an analytical posterior allowing accurate prediction after the model parameters have been inferred a single time. The authors apply the SPG to model tumor mutation rates and show that model parameters can be accurately inferred from high-dimensional epigenetic data. This is achieved through a combination of CNNs, GPs and MLE. The results are promising in detecting tumor drivers such as genes, regulatory structures and base-pairs.\n\nThe paper is well written, with motivation described and prior literature being discussed. Figure 1 is well laid out to drive home the point. \n\nComments:\n\n* Section 2: are all the distributions univariate? If not, a table giving dimensions will be helpful. What is the dimension of the covariates (eta_R)? It would also help if a plate diagram of the generative model is given. \n\n*  Why was CNN chosen for dimensionality reduction?\n\n*  Were other non neural architectures used for dimensionality reduction?\n\n*  What is the significance of ‘735’ epigenetic tracks? \n\n*  How valid is the assumption that events are distributed independently in the mutation space? I think this is too restrictive but at the same time simplifies the problem. A discussion of this simplification will be essential for the medical/biology audience.\n\n*  Title of paper has 'causes'. Title on ICLR webpage has 'cusses'. Please correct if possible.\n\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "-P9YofAKHdT", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Short summary: This paper introduces a split-Poisson-Gamma model to capture discrete-time, integer-valued stochastic processes at multiple scales. Although it seems to be simple and incremental compared with Poisson gamma models, this novel method may has some impacts in modeling mutation rates and identifying genomic elements that drive tumor emergence. Thus, I vote for weak acceptance.\n\n-quality: the technique contribution appears to be simple, incremental but could be useful in real applications of detecting cancer-related mutations.\n\n-clarity: the most parts of the paper is clearly written. Nonetheless, I find the paper is improvable by clearly describing the input data in Sec.2.\n\n-originality: to my knowledge, it is the first attempt to develop split-Poisson gamma model for discrete-time, integer-valued stochastic processes at multiple scales. \n\n-significance: it appears to be an incremental contribution in the family of Poisson gamma models. I think the split-Poisson gamma model is a sensible and useful method in detecting cancer-related mutation rates.\n\nPros: a useful and very sensible model for detecting cancer related mutation rates.\n\nCons: A clear description of the input data is missing in the beginning of Sec.2. I think it could be better if the authors can position the split-Poisson gamma model among the closely related Poisson gamma models, to clearly distinguish the main differences and why the novel model outperforms the others. \n\nIn addition to estimating mu_R and sigma_R^2, how to place hyper priors over these parameters , and perform MAP estimations?\n\nI agree the detection of cancer-causing mutations is a significant application of the multi-resolution modeling of discrete stochastic processes. Nonetheless, I would suggest the authors to discuss the broad applications of the split-Poisson gamma (SPG) model for a venue like ICLR. \n\ncomments:\n-\"NB\" notations are used in eq.(6,7) but are not defined for negative binomial yet.\n-you miss a \"Pr\" for lambda_R in Eq.(9).\n-can you present the complete procedure  to perform parameter inference in the supplement?\n-In C.1 of pp14, it should be NB(M_i;alpha_R, 1/(p_i \\theta_R+1)).", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Reviewers agree that the paper excels in providing a principle pipeline that combines CNNs and GPs with a Poisson-Gamma distribution to provide a generic approach for multiresolution modelling of tumour mutation rates. As a whole such combination of techniques addresses a key challenge in computational biology that also scales to large datasets. 
Decision: Accept (Poster)