GenoGAM 2.0: scalable and efficient implementation of genome-wide generalized additive models for gigabase-scale genomes

Background GenoGAM (Genome-wide generalized additive models) is a powerful statistical modeling tool for the analysis of ChIP-Seq data with flexible factorial design experiments. However large runtime and memory requirements of its current implementation prohibit its application to gigabase-scale genomes such as mammalian genomes. Results Here we present GenoGAM 2.0, a scalable and efficient implementation that is 2 to 3 orders of magnitude faster than the previous version. This is achieved by exploiting the sparsity of the model using the SuperLU direct solver for parameter fitting, and sparse Cholesky factorization together with the sparse inverse subset algorithm for computing standard errors. Furthermore the HDF5 library is employed to store data efficiently on hard drive, reducing memory footprint while keeping I/O low. Whole-genome fits for human ChIP-seq datasets (ca. 300 million parameters) could be obtained in less than 9 hours on a standard 60-core server. GenoGAM 2.0 is implemented as an open source R package and currently available on GitHub. A Bioconductor release of the new version is in preparation. Conclusions We have vastly improved the performance of the GenoGAM framework, opening up its application to all types of organisms. Moreover, our algorithmic improvements for fitting large GAMs could be of interest to the statistical community beyond the genomics field.

of basis functions, distributions and further features for variable and smoothness selection. In its general form, the implementation for fitting a GAM minimizes a cost function using iterations whose time complexity are quadratic in the number of parameters. Moreover, the time complexity of the implementation for estimating the standard errors of the parameters, which are required for any statistical significance assessment, is cubic in the number of parameters. To allow the fitting of GAMs on complete genomes, which involves millions of parameters, we had proceeded with a tiling approach [3]. Genome-wide fits were obtained by fitting models on tiles, defined as overlapping genomic intervals of a tractable size, and joining together tile fits at overlap midpoints. With long enough overlaps, this approximation yielded computation times linear in the number of parameters at no practical precision cost. Furthermore, it allowed for parallelization, with speed-ups being linear in the number of cores.
Nonetheless, application of the current implementation remains limited in practice to small genomes organisms such as yeast or bacteria, or to selected subsets of larger genomes. A genome-wide fit for the yeast genome (ca. 1 million parameters) took 20 hours on a 60-core server. Fits for the human genome could only be done for chromosome 22, the smallest human chromosome.
Here we introduce a new implementation of GenoGAM that is 2 to 3 orders of magnitude faster. This is achieved by exploiting the sparsity of the model and by using out-ofcore data processing. The computing time for parameter and standard error estimation, as well as the memory footprint, is now linear in the number of parameters per tile. The same genome-wide fit for yeast is now obtained in 13 min on a standard 8-core desktop machine. Wholegenome fits for human datasets (ca. 300 million parameters each) are obtained in less than 9 hours on the same 60-core server.
Before describing the new implementation and results, we provide some necessary mathematical background.

GenoGAM models
In a GenoGAM model, we assume ChIP-Seq read counts y i at genomic position x i in the ChIP-Seq sample j i to follow a negative binomial distribution with mean μ i and dispersion parameter θ: where the logarithm of the mean μ i is the sum of an offset o i and one or more smooth functions f k of the genomic position x i : The offsets o i are predefined data-point specific constants that account for sequencing depth variations. The elements z j i ,k of the experimental design matrix Z is 1 if smooth function f k contributes to the mean counts of sample j i and 0 otherwise. A typical application is the comparison of treatment versus control samples, for which a GenoGAM model would read: where z j i = 0 for all control sample data points and z j i = 1 for all treatment sample data points. The quantity of interest in such a scenario is the log fold-change of treatment versus control at every genomic position The smooth functions f k are piecewise polynomials consisting of a linear combination of basis functions b r and the respective coefficients β (k) r : where b r are cubic B-splines, which are bell-shaped cubic polynomials over a finite local support [8]. The column of the n×p k matrix X k , where p k is the number of basis functions in smooth f k , represents a basis function b r evaluated at each position x i . Typically all smooth functions have the same bases and knot positioning, implying that all X k are equal to each other. Consequently, the complete design matrix X is the Kronecker product of the experimental design matrix Z and X k .
where X = Z ⊗ X k and the vector β is the concatenation of all β (k) . The fitting of the parameters β is carried out by maximizing the negative binomial log-likelihood plus a penalty function: where S is a symmetric positive matrix that approximately penalizes the second order derivatives of the smooth functions. This approach is called penalized B-splines or P-splines [9]. The I term adds regularization on the squared values of the β's, which is particularly useful for regions with many zero counts. The smoothing parameter λ controls the amount of regularization. Both the smoothing parameter λ and the dispersion parameter θ are considered as hyperparameters that are estimated by cross-validation [3]. Newton-Raphson methods are used to maximize Eq. (6). The idea is to iteratively maximize quadratic approximations of the objective function around the current estimate. The current parameter vector β t is updated as: where the negative inverse Hessian H −1 (β t ) captures the local curvature of the objective function, and the gradient vector f (β t ) captures the local slope. The iteration stops when the change in the log-likelihood or the norm of the gradient of the log-likelihood falls below a specified convergence threshold. Because the negative binomial distribution with known dispersion parameter θ is part of the exponential family, the penalized log-likelihood is convex and thus convergence is guaranteed.

Standard error computation
For the purpose of statistical testing, variance of the smooth estimates are also needed. These are of the form: where H k is the Hessian with respect to the parameters β (k) and can be simply extracted from H.

Remarks on sparsity
The Hessian H is computed as: with W a diagonal matrix [6]. However, the number of nonzeros for each row of the design matrix X is at most 5 times the number of smooth functions because every genomic position x i is overlapped by 5 cubic B-splines b r only. Moreover, the penalization matrix S only has 5 nonzeros per row, as it encodes the second-order difference penalties between coefficients of neighboring splines [9]. Hence, the matrices X and S, and therefore H, which appears in the majority of the computations via Eqs. (7) and (8), are very sparse. Here we make use of the sparsity of these matrices to drastically speed up the fitting of the parameters.

Workflow
Data preprocessing consists of reading raw read alignments from BAM files, centering the fragments, computing the coverage y, and splitting the data by genomic tiles (Fig. 1). Afterwards, normalization factors for sequencing depth variation are computed using DESeq2 [10]. In the new version of GenoGAM we store the preprocessed data in HDF5 files [11] through the R packages HDF5Array [12] and rhdf5 [13]. This allows writing in parallel as the data is being preprocessed, which reduces the memory footprint of this step. For all subsequent matrix operations the Matrix package is used, which implements routines for storage, manipulation and operations on sparse matrices [14].
Fitting GenoGAM models on tiles is achieved by the Newton-Raphson algorithm (Eq. 7). This is done on few representative tiles during cross-validation in order to identify optimal hyperparameters λ and θ, and subsequently when fitting the model on the full dataset.
The variance of the smooth estimates (Eq. 8) is obtained using the sparse inverse subset algorithm as detailed in a subsection below. The implementation is based on the R package sparseinv [15], which wraps relevant code from the SuiteSparse software [16]. As in the previous GenoGAM model [3], fitting on different tiles is conducted in parallel. The result objects for the fits, variances and parameters are initialized prior to fitting on hard drive. This allows the processes to write results in parallel on the fly, ensuring fast computation and low memory footprint. The HDF5 storage is further optimized for reading time by adjusting HDF5 chunk size to the size of the tiles (for preprocessed count data) and chunks (for fits and variance). As HDF5 is not process-safe on R level, writing is serialized by a queuing mechanism.
The parallelization backend is provided by the R package BiocParallel. It offers an interface to a variety of backends and can be registered independently of GenoGAM. Parallelization is performed over chromosomes during the read-in process. Over tuples of folds and tiles during cross-validation process and over tiles during fitting process. Because some backends have a particular long start-up time, the use of many processes might end up dominating computation time. Specifically during cross-validation on small and limited number of regions, this might pose a problem. Therefore an optimal number of workers is automatically obtained and registered by the cross-validation function and reset on exit.

Newton-Raphson implementation for sparse matrices
We estimate the parameters β by maximizing the penalized log-likelihood using the Newton-Raphson iteration (Eq. 7). Due to the sparsity of the matrices X, D and S, H is sparse and cheap to compute. The inverse is never explicitly formed. Instead the linear system is solved by a direct solver using the SuperLU library [17]. Furthermore all matrices are stored in a sparse format, avoiding redundant storage of zeros.
Our new fitting algorithm differs from the one of mgcv in two ways. First, mgcv uses Iteratively Reweighted Least Squares, a Newton-Raphson method that employs the Fisher information matrix I, defined as the negative expectation of the Hessian H, instead of the Hessian in the iteration (Eq. 7). However, this did not lead to any measurable differences in the fitted parameters. Second, mgcv uses QR decomposition of the design matrix X [6]. However, general QR decomposition destroys the sparse Fig. 1 Schematic overview highlighting the difference between GenoGAM 1.0 and GenoGAM 2.0: Raw BAM Files are read-in, pre-processed normalized and written to hard drive in HDF5 format. Moreover, normalization factors for sequencing depth variation are computed using DESeq2 [10]. The resulting object is the dataset upon which fitting is done. Then global hyperparameters are estimated by cross-validation and for each tile coefficients are estimated via Newton-Raphson and standard errors via sparse inverse subset algorithm. The final model is written as a new object to hard drive in HDF5 format. Note, that the schematic view is a simplification: The pre-processed dataset and the fitted model are not generated in memory and written to HDF5 in the end. Instead, all HDF5 matrices are initialized on hard drive directly and the writing is done on the fly. Blue (GenoGAM 1.0) and orange colors (GenoGAM 2.0) mark differences between both GenoGAM versions, simultaneously displaying the content of this paper structure of X. We have investigated the use of sparse QR decompositions but this was less efficient than our final implementation.

Variance computation using the sparse inverse subset algorithm
The Hessian H is sparse, but its inverse, the covariance matrix H −1 , usually is not. However, the variances of interest (Eq. 8) can be computed using only a subset of the elements of the inverse H −1 . Specifically, denoting for any matrix A: On the other hand, Eq. 9 implies that NZ(H) = NZ(X T WX) ∪ NZ(S) ∪ NZ(I). Since it follows that: Moreover, the nonzeros of the identity matrix I is a subset of the nonzeros of the second-order differences penalization matrix S [9]. Furthermore, the nonzeros of the second-order differences penalization matrix S, which penalizes differences between triplets of consecutive splines, is a subset of the nonzeros of X T X, since genomic positions overlap five consecutive splines when using cubic B-splines. Hence, NZ(H) = (l, j), ∃i, (l, j) ∈ C 2 i (X) . Together with Eq. 10, this proves the result.
Using only the elements of H −1 that are in NZ(H) applies to computing the variance of any linear combinations of the β based on the same sparse structure of X or a subset of it. Hence, it applies to computing the variance of the predicted value for any smooth function f k (x) or computing the variance of the derivatives of any order r of any smooth To obtain the elements of H −1 that are in NZ(H), we used the sparse inverse subset algorithm [18]. Given a sparse Cholesky decomposition of symmetric matrix A = LL T , the sparse inverse subset algorithm returns the values of the inverse A −1 that are nonzero in the Cholesky factor L. Since nonzero in the lower triangle of A are nonzeros in the Cholesky factor L [19], the sparse inverse subset algorithm provides the required elements of H −1 when applied to a sparse Cholesky decomposition of H. See also Rue [20] for similar ideas for Gaussian Markov Fields. To perform the sparse inverse subset algorithm, we used the R package sparseinv [15], itself a wrapper of relevant code from the SuiteSparse software [16].
Once the sparse inverse subset of the Hessian is obtained, σ 2 i can be computed according to Eq. (10) with a slight improvement: Because only the diagonal is needed from the final matrix product, the implementation does not perform two matrix multiplications. Instead, only the first product is computed, then multiplied element-wise with X T k and summed over the columns.

Results
Leveraging the sparse data structure allows for faster parameter estimation

Exact σ 2 computation by the sparse inverse subset algorithm
Alternatively to the direct computation of the inverse Hessian with consecutive computation of variance vector σ 2 , it is also possible to directly compute σ 2 . Here and hereafter the smooth function specific index k is dropped for simplicity. In a comment to the paper of Lee and Wand [21], a direct way to compute σ 2 without inverting H was proposed by Simon Wood [22]. The comment states, that in general, if y = Xβ, then Where P is the permutation matrix and L −1 is the inverted lower triangular matrix resulting from Cholesky decomposition of X T H −1 X. Figure 3 shows the comparison of both methods in time and memory on a single core, with the above proposed method depicted as "indirect" (blue). While both methods have linear memory footprint, the slope of the indirect method is around four times higher. The computation time is significantly in favor of the sparse inverse algorithm. This is because for every σ 2 i a triangular system has to be solved to obtain (XP T ) i L −1 . Although solving the complete system at once is faster, it had a high memory consumption when it came to increased number of parameters in our implementation. Thus the performance presented is based on batches of σ 2 i to obtain a fair trade-off between runtime and memory. Nevertheless, the difference remains around 2 orders of magnitude. Moreover runtime goes now linearly in practice for the sparse inverse subset algorithm compared to quadratically for the indirect method (dashed black lines fitted to the performance data).

Performance on human and yeast ChIP-Seq datasets
The previous version of GenoGAM could only be partially applied genome-wide for megabase-scale genomes such as the yeast genome and was impractical for gigabasescale genomes such as the human genome. A genomewide model fit with two conditions and two replicates each took approximately 20 h on 60 cores [3]. With computational and numerical improvement on one side and a data model largely stored on hard drive on the other side, runtime and memory requirements have dropped significantly. Figure 4 shows the runtime performance on seven human ChIP-Seq datasets with two replicates for the IP and one or two replicates for the control. The analysis was performed with 60 cores on a cluster, the memory usage never exceeded 1.5 GB per core and was mostly significantly lower. The overall results show that around 20 min are spent with pre-processing the data, which is largely occupied by writing the data to HDF5  files. One hour of cross-validation, to find the right hyperparameters and around 7 to 8 h of fitting, amounting to a total runtime of 8 to 9 h. It is also notable, that the transcription factors NRF1, MNT and FOXA1 include two controls instead of one, thus efficiently increasing the amount of data to fit by a third, but the runtime by around 40 min.
Additionally, the same yeast analysis is shown by running on a laptop with 8 cores for comparison to the previous version. The total runtime is around 13 min with the cross-validation significantly dominating both other steps (around 9 min). This is due to the fact, that the number of regions used is fixed at 20, resulting in 200 model fitting runs for one 10-fold cross-validation iteration. Hence, for a small genome like the yeast genome, hyperparameter optimization may take more time than the actual model fitting. Note, that during cross-validation the only difference between human and yeast analysis is the underlying data and the parallel backend. However the runtime on yeast is only 1/6 of the runtime in human. Both factors play a role in this: First, the parallel backend in the yeast run uses the Multicore backend, allowing for shared memory on one machine. While the human run uses the Snow (simple network of workstations) backend, which needs to initiate the workers and copy the needed data first, resulting in an overall greater overhead. Second, convergence on yeast data is generally faster due to higher coverage resulting not only in less iterations by the Newton-Raphson, but also during cross-validation.

Replication of previous benchmark analyses show equivalent biological accuracy
To demonstrate that GenoGAM 2.0 leads to the same results than GenoGAM 1.0 we have re-generated benchmark analyses of the first paper [3]. The first benchmark is a differential occupancy application that demonstrates that GenoGAM has greater sensitivity for same specificities than alternative methods ( Fig. 5a-b). The second benchmark shows that GenoGAM is on par with alternative methods to infer peak summit positions in ChIP-Seq data of transcription factors (Fig. 5c). Consistently, with the fact that GenoGAM 2.0 fits the same function than GenoGAM 1.0, the performance on these two benchmarks matched.
These improvements have required us to re-implement the fitting of generalized additive models, since GenoGAM 1.0 was based on an generic R package for fitting generalized additive models. We have restricted the implementation so far to the negative binomial distribution. Therefore, application to methylation data, which requires the quasi-binomial distribution, is not yet supported.

Conclusion
We have significantly improved the implementation of GenoGAM [3] on three main aspects: Data storage, coefficient estimation and standard error computation. We showed its runtime and memory footprint to scale linearly with the number of parameters per tiles. As a result, GenoGAM can be applied overnight to gigabase-scale Replication from our previous study [3] with GenoGAM 2.0. a Replication of figure 3A from our previous study [3]. ROC curve based on a quantile cutoff of 0.15 (15% of the genes are assumed to be true negatives). GenoGAM (orange and blue) has a constantly higher recall with a lower false positive rate. b Replication of Fig. 3b from our previous study [3]. Area under the curve (AUC) for all possible quantile cutoffs from 0 to 1 in steps of 0.01. GenoGAM 1.0 (blue) and GenoGAM 2.0 (orange) are almost identical and are thus largely overlapping. Up to a cutoff of 0.6, GenoGAM (orange and blue) performs consistently better than all competitor methods by around 0.03-0.04 points above the second best method (csaw and DESeq2, pink and green, respectively). The entire range of quantile cutoffs is shown out of completeness, reasonable values are between 0.15 and 0.25. c Replication of supplementary figure S9C from our previous study [3]. Proportion of significant peaks within 30 bp of motif center and 95% bootstrap confidence interval (error bars) for six ENCODE transcription factors (CEBPB, CTCF, USF1, MAX, PAX5, YY1) on chromosome 22 and for the yeast TFIIB dataset genome datasets on a typical lab server. Runtime for mega-base genomes like the yeast genome is within minutes on a standard PC. Finally, our algorithmic improvements apply to GAMs of long longitudinal data and can therefore be relevant for a broader community beyond the field of genomics.