Population sequencing data reveal a compendium of mutational processes in the human germ line

Abstract

Biological mechanisms underlying human germline mutations remain largely unknown. We statistically decompose variation in the rate and spectra of mutations along the genome using volume-regularized nonnegative matrix factorization. The analysis of a sequencing dataset (TOPMed) reveals nine processes that explain the variation in mutation properties between loci. We provide a biological interpretation for seven of these processes. We associate one process with bulky DNA lesions that are resolved asymmetrically with respect to transcription and replication. Two processes track direction of replication fork and replication timing, respectively. We identify a mutagenic effect of active demethylation primarily acting in regulatory regions and a mutagenic effect of long interspersed nuclear elements. We localize a mutagenic process specific to oocytes from population sequencing data. This process appears transcriptionally asymmetric.

Fig 1.. Inference of spatially-varying mutational processes in germline.

(A) Observed spatial variability of mutational spectrum is modeled as a number of mutational processes with specific spectra and spatially-varying intensities. (B) Strand-independent mutational processes have equal rates of complementary mutations at each locus. Strand-dependent mutational processes produce two unequal patterns of complementary mutation rates at loci depending on the strand orientation of a genomic feature. (C) Example of a predicted strand-independent process. Loadings of complementary mutations of mutational component 9 are highly similar. (D) Example of a predicted strand-dependent process. Loadings of complementary mutations between two mutational components 1 and 2 are highly similar and characterize a single mutational process (left). In contrast to strand-independent process (see Fig. 1C), loadings of complementary mutations of mutational component 1 (upper right) and component 2 (lower right) are almost uncorrelated. (E) Theoretical scale-loading limitations for detection of mutational processes shows potentially high range of processes that can be recovered with the proposed approach. Simulations of mutational processes at different scales (quantified as half-life of simulated Ornstein-Uhlenbeck process) and spatial loadings (fraction of spatially-varying mutations of a process among total mutations, scheme) were based on parameters from TOPmed dataset. Quality of recovery was assessed using maximum absolute correlation between spectra of each simulated component and reconstructed components. (F) Reflection matrix reveals strand-dependency of processes and separates biological signal from noise. Correlation of spectrum of one mutational component with reverse complementary spectrum of another demonstrates clear separation into self-correlated components (5,8, 9,14) and pairs of mutually correlated components (1/2, 3/4, 6/7, 10/11, 12/13).

Fig 2.. Mutational processes are associated with distinct genomic features.

(A) Heatmap of correlations of intensities with genome features shows diverse modes of associations (left). For strand-dependent processes two spatial characteristics were considered: intensity (int.), estimated as the sum of intensities of two components, and asymmetry (as.), estimated as the difference between intensities of two components. Fraction of mutational variance explained by each process (middle) and scale, estimated as the half-life of the autoregressive model (right) are shown. (B) (top) The spectrum of one of the two components comprising process 1/2; (bottom) An example of intensities of components 1 and 2, associated with non-transcribed strands on chromosome 1. The bars on the bottom of the panel depict gene bodies (colors: cyan if transcribed strand is the reference strand and orange otherwise). (C) (top) The spectrum of one of the two components comprising process 3/4; (bottom) The association between the asymmetry of process 3/4 (component 3 – component 4) and the direction of the replication fork measured as a gradient of replication timing. (D) The spectrum of component 5 (top), and its association with replication timing (bottom). (E) (top) The spectrum of component 10; (bottom) Intensity of component 10 among 13 consecutive 10 KB-long regions adjacent to the transcription end site (TES). Box plot shows component 10 intensity. Mean intensity of component 10 in each region shown as a red point. Mean fraction of the transcribed nucleotides per region shown on the bottom.

Fig 3.. Oocyte-specific mutational process.

(A) The spectrum of one of the two components comprising process 6/7. (B and C) Examples of two loci with high intensity of process 6/7, estimated as the sum of intensities of component 6 and component 7. Black dots on top of the panels mark windows of high intensity that we call “maternal regions” (see Methods). Red dots show de novo maternal clustered mutations from Halldorsson et al. (18). (D) Enrichment of maternal clustered de novo mutations from Halldorsson et al. in maternal regions. The fraction of each chromosome that is attributed to “maternal regions” (high intensity of process 6/7) is shown in black. The fraction of maternal clustered mutations located within maternal regions on each chromosome is shown in red; the difference in size between the red and black bars indicates enrichment of clustered mutations within “maternal regions”. (E and F) Zoom in view of process 6/7 intensity spikes around FHIT and CSMD1 genes on non-transcribed strands. Bars on the bottom depict gene bodies (colors: cyan if transcribed strand is the reference strand and orange otherwise). (G) Difference between C>G mutation rate on transcribed or non-transcribed strand of a gene compared to a 100 KB region flanking the gene. Red dots correspond to genes within maternal regions and black dots corresponds to genes outside of maternal regions. Density plots on the top and right summarize the distributions on the X and Y axes. (H) Ratio of parent-specific de novo mutation rates between first and last parent age quartiles is shown, estimated independently for “maternal regions” and for the rest of the genome. The error bars show the 95% confidence interval (95% CI) for the ratio of two binomial proportions test.

Fig 4.. Cytosine deamination and cytosine demethylation.

(A and D) The mutational spectrum of component 8 is dominated by CpG>TpG (top), while component 9 is dominated by CpG>GpG and CpG>ApG (bottom). (B, C, E and F) Association between the intensity of components 8 and 9 and cytosine methylation or cytosine hydroxymethylation. (G) Process 8 is inversely correlated with the density of CpG islands, while process 9 is positively correlated. Blue dots represent density of CpG islands across a 50 MB long region of chromosome 5. (H) Effect of CpG islands on mutation rate in CpG context and in cytosines outside of CpG context. (I) Mutation rate in CpG context at transcription factor binding sites located outside of annotated CpG islands, as determined by ChIP-seq peaks (see Methods), normalized to the genome average mutation rate. 95% binomial proportion confidence intervals are displayed in transparent lines. Higher levels of demethylation at these sites lead to the accelerated rate of CpG transversions. (J) Illustration of the biochemical mechanisms suggested for processes 8 and 9. Enzymatic oxidation of methylcytosine (5-mC) leads to hydroxymethylcytosine (5-hmC) (35), which after additional steps of oxidation should be removed by glycosylase, leaving an abasic site (AP). During DNA replication, AP sites will frequently be converted to CpG>GpG mutations and more rarely to CpG>ApG mutations, matching the spectra of process 9 (36). Alternatively, successful repair of AP sites creates non-methylated cytosines. Alternately, spontaneous deamination of methylcytosine creates a T to G mismatch, enhancing the rate of CpG>TpG mutations. While deamination should be prevalent in CpG sites with high methylation levels, the mutagenic effect of demethylation should be prominent in CpG islands.