Transcription start sites experience a high influx of heritable variants fuelled by early development

Transcription start sites experience a high influx of heritable variants fuelled by early development | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Transcription start sites experience a high influx of heritable variants fuelled by early development Miguel Cortés Guzmán , David Castellano , Clàudia Serrano Colomé , Vladimir Seplyarskiy , View ORCID Profile Donate Weghorn doi: https://doi.org/10.1101/2025.02.04.635982 Miguel Cortés Guzmán 1 Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology , Dr. Aiguader 88, Barcelona 08003, Spain 2 Universitat Pompeu Fabra (UPF) , Barcelona, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site David Castellano 1 Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology , Dr. Aiguader 88, Barcelona 08003, Spain 5 University of Arizona , Tucson, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Clàudia Serrano Colomé 1 Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology , Dr. Aiguader 88, Barcelona 08003, Spain 2 Universitat Pompeu Fabra (UPF) , Barcelona, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site Vladimir Seplyarskiy 3 Division of Genetics, Brigham and Women’s Hospital, Harvard Medical School , Boston, MA, USA 4 Department of Biomedical Informatics, Harvard Medical School , Boston, MA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Donate Weghorn 1 Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology , Dr. Aiguader 88, Barcelona 08003, Spain 2 Universitat Pompeu Fabra (UPF) , Barcelona, Spain Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Donate Weghorn For correspondence: dweghorn{at}crg.eu Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Mutations drive evolution and genetic diversity, but the impact of transcription on germline mutagenesis remains poorly understood. Here, we identify a hypermutation phenomenon at transcription start sites in the human germline, spanning several hundred base pairs in both directions. We link this TSS mutational hotspot to divergent transcription, RNA polymerase II stalling, R-loops, and mitotic—but not meiotic—double-strand breaks, revealing a recombination-independent mechanism distinct from known processes. Notably, the hotspot is absent in de novo mutation data. We reconcile this by showing that TSS mutations are significantly enriched with early mosaic variants often filtered out in de novo mutation calls, indicating that the hotspot arises during early embryogenesis. Mutational signature analysis reinforces these findings and implicates alternative non-homologous end joining and maternal mutation clusters. Our study provides the first detailed description of a germline TSS mutation hotspot, with broad evolutionary and biomedical implications. Introduction The mutation rate in human cells is influenced by DNA sequence composition and epigenetic factors. 1 , 2 , 3 Although germ cells and somatic tissue cells share many of these factors, there are fundamental differences in the observed distribution and nature of mutations, particularly due to exogenous and disease-related mutational processes active in the soma. 4 In addition to replication time, which is one of the most important determinants of mutation rate, transcription and the directionality of its effects are also debated. 5 , 6 Transcription affects mutation rate in at least two important ways: (1) via the joint effects of transcription-associated mutagenesis (TAM) and transcription-coupled repair (TCR) on the transcribed sequence itself and (2) via the effects of transcription factor (TF) binding in the regulatory sequence elements. 7 The balance of TAM and TCR on expressed genes is a topic of great interest. In the soma, gene expression is generally thought to correlate negatively with observed mutation density, likely due to a dominant role of TCR, especially in tumours with a strong influx of mutation-inducing DNA damage. 2 , 8 In the germline, the consensus is not as clear. A recent study reported a largely negative correlation between transcription and mutation density, similar to that observed in the soma. 9 , 10 However, this effect was not observed in another study 4 and reversed when some potential confounding factors were taken into account, 11 consistent with previous reports. 12 , 13 Hence, it is currently an open question how large the actual impact of transcription on mutation rates in the germline is when all other covariates of mutation rate are taken into account. Regarding the influence of TFs, it has been shown that promoter-proximal TF binding sites in the soma have an increased mutation density, especially in melanomas and binding sites belonging to CTCF. 14 , 15 , 16 This effect is likely due to a combination of limited DNA repair and, in particular, increased DNA damage. 17 , 18 Interestingly, Perera & Poulos et al . (2016) found that TF binding alone cannot be the sole reason for the increased mutagenicity at DNAse I hypersensitive sites (DHSs) and speculated that transcription initiation must play an important role. 16 In the germline, an excess of mutations around testis-active promoter-proximal TF binding sites, especially for T>G mutations, was recently found. 19 , 20 To shed more light on these open questions, we investigated the effects of transcription on mutation rate by examining the variability of observed mutation density in transcribed regions and their genomic neighbourhood. Using extremely rare variants (ERVs) from two large human whole-genome sequencing cohorts, we uncovered a remarkable mutational hotspot of non(CpG>TpG) mutations near the transcription start site (TSS) in the human germline, extending several hundred base pairs in both directions. Surprisingly, the hotspot is not detectable in de novo mutations (DNMs) from family sequencing data. We show that this apparent discrepancy is resolved by a significant enrichment of the TSS with early mosaic mutations, which are largely removed in family sequencing mutation calling pipelines, unveiling early development as a key factor in TSS mutagenesis. Consistently, analysis of the mechanistic factors driving TSS mutagenesis showed no association with meiotic double-strand breaks, but instead revealed that the hotspot is associated with divergent transcription, stalling of RNA polymerase II (RNAP II), R-loops and mitotic double-strand breaks. Mutational signature decomposition further corroborates the role of mosaic variants and suggests involvement of non-canonical homologous recombination repair, transcription-dependent mutational processes, and a process typically associated with clustered mutations that arise in the female human germline. Results Transcription start sites of protein-coding genes show an excess of extremely rare polymorphisms To quantify the direct and indirect effects of transcription on the mutation rate in the human genome, we analysed four sets of mutations: Extremely rare variants that segregate in the human population obtained from (1) gnomAD 21 and (2) the UK Biobank (UKBB, www.ukbiobank.ac.uk , sample allele frequency AF < 0.01%), (3) de novo mutations from family sequencing data ( Methods ) and (4) somatic mutations from cancer genomes. 22 Like other studies, we considered ERVs as a proxy for neutral germline variations, as they are only slightly affected by the effects of selection. 23 , 24 We excluded hypermutable tumours from the pan-cancer cohort, as they exhibit greater sequence context dependencies than our model considers. 25 , 26 Furthermore, we subdivided the mutations into the two main classes of non(CpG>TpG) and CpG>TpG mutations, as the latter have a 1-2 orders of magnitude higher average mutation rate due to methylation-dependent mutagenicity. 27 For each set of mutations, we then calculated the relative mutation density in non-overlapping genomic windows of 100 or 1000 base pairs (bp) in length. We anchored our analysis at the transcription start site (TSS) and transcription termination site (TTS) of 14763 non-overlapping protein-coding genes and analysed up to 50 kilo-base-pairs (kb) upstream and downstream of the TSS and TTS. Our measure of relative mutation density per window is the ratio between the observed and expected number of mutations averaged over genes, µ. The expected number of mutations takes into account the differences in mutation propensity due to DNA sequence composition. It is based on a 5-mer sequence context model, which is calculated separately for the four different mutation data sets and the two mutation classes from the respective genome-wide average mutation probabilities. Because we treated complementary mutations separately and orientated them along the direction of transcription, this model also accounts for the effects of transcription strand bias. 24 Since we found that both gene length and distance to other genes are positively correlated with µ , likely due to their correlation with replication time, 28 we also defined a second measure of sequence-corrected mutation density, µ ′ , for visualisation ( Methods ). This measure was designed to capture the variability of mutation density at length scales ≤ 150 kb in and around transcribed regions that is neither caused by DNA sequence composition nor by the larger regional variability in mutation rate. The most remarkable feature of the relative mutation density in transcribed and neighbouring untranscribed regions that we found is a striking mutational hotspot at the TSS for non(CpG>TpG) ERVs, Figure 1a-b . This localised excess reaches up to 14% at the 1-kb scale. Zooming in to 100 bp, the relative density of ERVs in the first 100 transcribed sites is increased by up to 35%. There is also evidence of a hotspot of cancer mutations at the TSS, but it is much narrower (≈300 bp) and shifted upstream towards untranscribed regions, Figure 1c , consistent with previous results. 16 Stratification by cancer type showed that this signal is driven by several tumour types, including bladder, breast, oesophageal, head and neck, lung, lymphoid, ovarian, pancreatic and gastric cancers ( Extended Data Figures 1-2 ). In lymphoid tumours, the mutational excess is restricted to the first 2 kb downstream of the TSS, consistent with somatic hypermutation in the variable domains of immunoglobulin genes, 29 whereas in all other cancer types with a hotspot, the mutational excess has an upstream component. Notably, liver tumours show the opposite pattern, with a clear deficit of mutations, especially in the first 1-2 kb downstream. Download figure Open in new tab Figure 1: Extremely rare variant and pan-cancer mutation density around the transcription start and termination site. Average number of mutations across 14763 protein-coding genes divided by the expectation based on the 5-mer sequence context and by the gene-level mutation density, µ 0 , upstream and downstream of the TSS (orange and blue, respectively) and upstream and downstream of the TTS (green and yellow, respectively). Each dot represents µ 0 in a 1-kb window (left) or 100-bp window (right). Error bars represent the 90% confidence intervals across 100 bootstrap replicates. (a-c) non(CpG>TpG) mutations using (a) UKBB ERVs, (b) gnomAD ERVs and (c) pan-cancer somatic mutations. (d-f) CpG>TpG mutations using the same data sets. Previous work has shown that TF binding in cancer and in the germline can lead to locally increased mutation at promoter DHSs, suggesting that TSS hypermutability may be related to transcription. 16 , 20 Since transcription in the human genome is not restricted to protein-coding genes, we applied the same methodology as for protein-coding genes to three types of non-protein-coding genes: 3991 upstream antisense or “divergent” lncRNAs, 8454 intergenic lncRNAs and 1660 pseudogenes. The promoters of intergenic lncRNAs are usually embedded in enhancers, whereas the promoters of divergent lncRNAs are shared with the promoters of protein-coding genes and transcription occurs upstream and in the opposite direction. 30 Our analysis revealed an even larger TSS hotspot of non(CpG>TpG) mutations on divergent lncRNAs than for protein-coding genes, both for ERVs (47% excess for gnomAD and UKBB) and for somatic cancer mutations (49% excess) at 100 bp ( Extended Data Figures 3-8 ). Conversely, mutations on pseudogenes and on intergenic lncRNAs showed only a slight excess (<11%) near the TSS at the 100-bp length scale for ERVs and no significant excess for the pan-cancer data set. Although divergent lncRNAs are expected to be transcribed more frequently than intergenic lncRNAs, these results may suggest that, in addition to gene expression itself, bidirectional transcription is mechanistically related to the TSS-proximal mutation hotspot. In contrast to non(CpG>TpG) mutations, we found a deficit of CpG>TpG mutations in all mutation sets near the TSS at the 1-kb level, Figure 1d-f , as expected due to a lower mutation rate of unmethylated CpG islands in the vicinity of the TSS of expressed genes. 31 Interestingly, against the background of this general deficit, ERVs show a small relative excess of CpG>TpG variants in the first 100 bp downstream of the TSS. Finally, for all mutation sets and classes, we found on average a lower µ ′ in transcribed than in neighbouring untranscribed regions (see Extended Data Figure 9 for results showing µ ), which can be explained in large part by the strong entanglement of transcription and replication landscapes, especially in the soma 32 ( Extended Data Figure 10 ). To test the effects of negative selection on our results, we repeated the same analyses without exons and conserved non-coding sequence elements. If strong negative selection removed mutations at the TSS, we would expect to observe an even more pronounced hotspot when we exclude conserved sites. However, the TSS hotspot for ERVs is reduced after removal of exons and conserved sequences ( Extended Data Figures 11-12 ), consistent with the mutation rate being primarily a function of distance from the TSS. Early mosaic mutations are significantly overrepresented near the transcription start site A more direct indicator for the germline mutation rate than ERVs are de novo mutations (DNMs). Surprisingly, however, DNMs did not show a significant excess of non(CpG>TpG) mutations near the TSS, which could not only be explained by the smaller sample size, especially downstream of the TSS, Figure 2a-b . We hypothesised that a possible explanation could lie in the filtering applied to de novo mutation calls from family sequencing data. First, mutations that match common variants that segregate in the human population are often filtered out, 34 but this should remove a negligible number of variants. Second, by construction early mosaic variants that occurred in the parents are removed during mutation calling from family sequencing data, Figure 2c . Mosaic variants have been found to contribute to several of the germline mutational signature components identified in Seplyarskiy & Soldatov et al . (2021), 24 which represent independent mutational processes involved in ERV mutagenesis. Therefore, we next decomposed the ERV mutational profile near the TSSs of protein-coding genes into these germline signature components, Figure 2d . We found that the relative enrichment of a component near the TSS compared to the gene body correlates significantly with the enrichment of that component with mosaic variants compared to de novo variants, as identified in Seplyarskiy & Soldatov et al . ( p = 0.006, Pearson), Figure 2e . Accordingly, three of the four components significantly enriched in mosaic variants are increased 1.5- to 4.8-fold near the TSS (components 5/6, 7 and 11), supporting the hypothesis that the TSS is influenced by somatic mutagenesis in early development. The two components most enriched in the comparison of mosaic against de novo variants (components 7 and 11) are also globally most enriched in a comparison of ERVs and DNMs, although the correlation across all signatures is not significant. Download figure Open in new tab Figure 2: Early mosaic non(CpG>TpG) variants are enriched at the TSS. (a-b) Relative mutation density µ ′ upstream (orange) and downstream (blue) of the TSS for gnomAD ERVs (top panels) and DNMs (bottom panels) in 1-kb windows (left) and 100-bp windows (right): (a) non(CpG>TpG) mutations, (b) CpG>TpG mutations. (c) Illustration of mosaic mutation timing and filtering in standard family sequencing datasets, inspired by Fig. 1 in Jonsson et al . (2021). 33 PGCS=primordial germ cell specification. (d) Decomposition of ERVs in 100-bp bins into germline mutational signature components. 24 (e) Enrichment of mutational signature components with mosaic mutations (defined in 24 ) as a function of the ratio of the component weight in the first and the 15th 1-kb bin downstream of the TSS (left) and the ratio of weights estimated from [− 10, 10] kb around the TSS in ERVs and DNMs (right, ϵ = 0.005). (f) Ratio of observed and expected number of non(CpG>TpG) early mosaic (maroon), late mosaic (turquoise), de novo (yellow) and extremely rare (lilac) variants around TSSs of protein-coding genes, averaged in bin brackets ± [1], ± [2,5] and ± [6,50] kb. Expected number based on a mononucleotide model. Confidence intervals derived from a Poisson model fit. ** p = 0.005. Given these compelling indications that mosaic variants might explain the TSS mutational hotspot, we compiled a data set of mosaic mutations from 11 published studies 33 - 43 , separating into early and late mosaic variants. Remarkably, we found a highly significant increase of 52% in the density of early mosaic mutations immediately downstream of the TSS ( p = 0.005), Figure 2f . This strongly suggests that the difference between the 14% excess in ERV density and the 4% excess observed in DNMs 1 kb downstream of the TSS is indeed due to mosaic variants. At the same time, the observed ERV excess in the 1 kb upstream of the TSS (7%) is within the bootstrap intervals of the DNM data, indicating that the two data sets are compatible in this genomic region, consistent with no significant deviation of mosaic variant density from the expectation in this bin. Taken together, these results suggest that removal of mosaic mutations in DNM data sets explains the absence of the mutation hotspot downstream of the TSS when using DNMs instead of ERVs as a proxy for germline mutagenesis. The TSS germline mutational hotspot is associated with mitotic double-strand breaks, divergent transcription, RNA polymerase II stalling and R-loops Next, we aimed to identify genomic and epigenetic features that affect the hypermutation phenomenon at the TSS of protein-coding genes. First, we performed a multiple negative binomial regression analysis of the sequence-corrected ERV mutation density, µ tb , across all protein-coding genes t at the 1-kb window level using 38 feature variables. The dependence of each feature on position relative to the TSS (i.e., on the bin ID b) was modelled with interaction terms ( Methods ). To extract the differential effect of each feature on mutation density at the TSS, we compared the coefficient of the interaction term at the TSS with the one far away from it, taking into account the covariance structure between the coefficients. Two features were significantly associated (adjusted p-value q < 0.05) with increased ERV density in the first 1 kb on both sides of the TSS, which was replicated in both the gnomAD and UKBB ERV data sets: (1) annotated sites of somatic double-strand breaks (DSBs) 44 and (2) G/C content, Figure 3a . Five additional features correlated significantly positively with ERV mutagenicity 1 kb downstream of the TSS: (3) H3K27ac, a histone mark found at the TSSs of actively transcribed genes, 45 (4) TATA box promoter, (5) the ratio of PRO-seq read density at the TSS and further downstream, a measure of RNAP II stalling, (6) G/C skew, a measure of the probability of R-loop (i.e., transcription-associated DNA-RNA hybrid) formation, 46 and (7) divergent transcription, defined based on GRO-cap data. 47 This last finding was supported by the strongest predictor of mutational excess 1 kb upstream of the TSS: (8) co-localisation of the protein-coding gene with an upstream antisense lncRNA gene, an alternative measure of divergent transcription. Download figure Open in new tab Figure 3: Determinants of germline mutation density around the transcription start site. (a) Difference in standardised multiple negative binomial regression coefficient of regressor bin interaction between bin brackets [1] and [6,50] kb upstream (orange) and downstream (blue) of the TSS for non(CpG>TpG) mutations. All significant terms replicating across UKBB and gnomAD ERVs that also have a significant total effect size at the TSS in at least one of the two data sets are shown ( q ≤ 0.05). Confidence intervals derived from coefficient standard errors using single-step intervals correction. 48 dens.=density, div. prom.=divergent promoter. (b) Same as (a) , but using single negative binomial regression with interaction terms for early mosaic (top left), late mosaic (top right) and de novo (bottom) variants. The associations with H3K27ac and divergent transcription confirm the notion that both transcription itself and its bi-directionality independently increase germline mutation density at the TSS, consistent with the large mutational excess at divergent lncRNA promoters ( Extended Data Figures 3 ). R-loops and RNAP II stalling consolidate the link of the mutational hotspot to transcription, possibly capturing related phenomena. 49 , 50 Potentially the most remarkable finding is the strong association of the ERV hotspot in the first 1 kb on both sides of the TSS with double-strand breaks ( q ∈ [0, 10 - 2 ] across UKBB and gnomAD). Since our DSB estimates are based on measurements in a healthy somatic cell line 44 and a recent study found mutational hotspots in the vicinity of putative sites of meiotic DSBs, 51 we also included the annotations of those sites in our feature list (“PRDM9”). While we indeed found a significant positive genome-wide effect of PRDM9 on non(CpG>TpG) mutation probability in both ERV data sets across transcribed and untranscribed regions ( . Extended Data Figures 13-14, Supplementary Table 1 ), the signal was not statistically significantly associated with the TSS ( Supplementary Table 2 ), suggesting that the association of the TSS mutational hotspot with somatic DSBs has a different aetiology to PRDM9-mediated events. We next aimed to test if mosaic variants recapitulated the associations with genomic and epigenetic features observed for ERVs by regressing separately on each of the 38 feature variables and their interaction with bin ID. Remarkably, despite the much smaller sizes of these mutation data sets (early mosaic: 0.002%, late mosaic: 0.001%), this identified two of the seven features that were also significantly positively associated with ERV mutation density in the first 1 kb downstream of the TSS: somatic double-strand breaks (early mosaic) and PRO-seq read density ratio (late mosaic), Figure 3b . These results suggest a clear timeline for TSS mutagenesis, with error-prone double-strand break repair during early developmental somatic cell divisions followed by RNAP II stalling at developmental gene promoters driving mutation accumulation. DNMs confirm the significant association with RNAP II stalling (“PRO-seq ratio”, q < 0.05) downstream, which may in part be driven by late mosaic variants, as the majority of late mosaic variants in our DNM data set could not be removed due to the limited number of studies with families with multiple progeny. Finally, we recovered a significant association of DNMs with DSBs and divergent transcription upstream of the TSS, as we had already observed for ERVs, supporting an important role for late germline mutagenesis in the mutation excess in promoter regions. The mutation hotspot at the TSS is associated with extraordinary mutational processes Given its unusual nature, we hypothesised that the mutational hotspot near the TSS is due to mutational processes that diverge from background processes ubiquitously active in the human germline and soma. 4 Therefore, we decomposed the mutational profile in each window around the TSSs of protein-coding genes into COSMIC mutational signatures for single base substitutions (SBS), 52 Figure 4a-c . Although these mutational signatures were derived from cancer sequencing data, we hypothesised that there might be some overlap and informative similarities with germline mutational processes beyond the background signatures SBS1 and SBS5. We quantified both the relative contributions of mutational signatures ( Figure 4a ) and the absolute number of mutations attributable to each signature (corrected for the target size, Figure 4b ). Download figure Open in new tab Figure 4: Mutational signatures and transcription strand biases affecting the TSS mutational hotspot. (a) Mutational signature decomposition in 1-kb (left) and 100-bp (right) windows upstream and downstream of the TSS of protein-coding genes. Results are shown for ERVs (UKBB and gnomAD), DNMs and PCAWG pan-cancer mutations (top to bottom). Only bins with more than 1000 mutations pooled across genes were used (others are marked in grey). (b) Mutation density stratified by signature for mutations in bins upstream (“uTSS”) and downstream (“dTSS”) of the transcription start site, shown for the 1-kb bins directly neighbouring the TSS (“bin=1”) and for mutations pooled across bins 6 to 50 kb (“bin>5”). (c) Mutational profile in the five bins around the TSS [-2,3] kb after subtraction of SBS1, SBS5 and SBS40b,c. (d) Relative contribution (quantiles) of non-gene-level genomic and epigenetic features used in the regression analyses as a function of position relative to the TSS. (e) non(CpG>TpG) ERV mononucleotide mutation density, stratified by transcription strand (coding or template) and replication strand (leading or lagging) in 1-kb windows around the TSS. The decomposition of the ERV profiles revealed four mutational signatures that occur exclusively in the region around the TSS: SBS3, SBS11, SBS40b and SBS40c. SBS3 has been associated with deficiencies in homologous recombination repair (HRR) of double-strand breaks in cancer tumours. 53 A likely candidate for alternative repair of such DSBs is the polymerase theta-mediated end joining (TMEJ) pathway, 54 and we have recently shown that indels associated with this repair pathway are co-localised with SBS3 point mutations. 55 In addition, we have linked SBS40 (whose subcomponents are SBS40b and SBS40c) to TMEJ during aborted HR or stalled replication. 55 Thus, the mutational signatures strongly suggest non-canonical DSB repair in the TSS germline mutational hotspot, consistent with the significant association of TSS hypermutability with somatic DSBs ( Figure 3 ). In both ERV data sets, we also found several other non-background mutational signatures with increasing contribution closer to the TSS: SBS12, SBS16, SBS19, SBS30 and SBS39, all of which are also present near the TSS in the DNM data set. Of these, the largest proportion of mutations near the TSS originates from SBS39 (10%). This signature is dominated by C>G mutations and shows a high similarity to a germline mutational signature component associated with maternal de novo mutation clusters (“component 8/9” in Seplyarskiy & Soldatov et al . 24 ). To test the de facto association between SBS39 and component 8/9, we stratified the mutations by chromosomes enriched in maternal mutation clusters (“maternal” chromosomes 8, 9, 15 and 16 34 , 38 , 24 ) and found that SBS39 mutation density on these chromosomes is systematically increased in all three germline mutation data sets ( Extended Data Figure 15 ). While signature SBS39 is present on transcribed sequences further downstream (and more so than in untranscribed regions, in agreement with 24 ), the much larger number of mutations associated with this signature near the TSS suggests that the underlying mutational process is particularly influenced by TSS-proximal events. As it has been hypothesised that maternal clustered mutations are associated with meiotic crossover events, 34 , 38 the latter could be the driving mechanism behind SBS39 enrichment at the TSS, even though we found no interaction of PRDM9 or recombination rate with the TSS mutation excess ( Figure 3 ). However, both recombination rate and PRDM9 levels (measuring the DNA occupancy of the homology search protein DMC1 51 ) are relatively low in the 1-kb bins neighbouring the TSS compared to sites further upstream or downstream, Figure 4d . This suggests that meiotic double-strand breaks, in addition to their negligible effect on mutational excess near the TSS, also have a small effect on the relative weights of mutational signatures near the TSS. In addition to C>G, the TSS mutational profile showed pronounced contributions from C>T and T>C mutations after subtraction of background processes, Figure 4c . T>C variants are in agreement with the increased weight of signatures SBS12 and SBS16 near the TSS. In cancer, these signatures are primarily observed in liver tumours and co-occur so frequently that they contaminate each other as well as the background process SBS5, 56 suggesting a more fundamental role in mutagenesis. Notably, SBS5 is relatively less active near the TSS. Both SBS12 and SBS16 exhibit a large transcription strand bias, suggesting a differential effect of TAM and TCR. 57 , 58 To link this to germline mutations, we calculated the transcription and replication strand biases for all mononucleotide mutation types, indeed showing strong biases for T>C ERVs, Figure 4e . In contrast to T>C variants, the large contribution of C>T mutations to the background-corrected mutational profile near the TSS has no obvious counterparts among somatic mutational signatures. Instead, we hypothesise that this accounts for the identification of several signatures that together approximate the observed C>T pattern in Figure 4c , including SBS11, SBS19, SBS30 and SBS32. In cancer tumours, the characteristic APOBEC-induced C>G and C>T mutation patterns in TCA and TCT contexts (corresponding to SBS2 and SBS13) dominate the mutational spectrum around the TSS, which is consistent with previous reports. 59 Besides T>C, we also found large transcription strand biases for all other mutation types except T>G (see Extended Data Figure 16 for cancer mutations). As for the TSS mutational hotspot, the mutation rate is accelerated both upstream and downstream for almost all mutation types near the TSS (see Extended Data Figure 17 for the 100-bp scale). In a previous report, the large TSS-proximal excess of C>T mutations on the coding strand was associated with damaged ssDNA during R-loop formation. 57 However, we observe that in addition to an increase on the coding strand, the C>T mutation density is decreased on the template strand, which requires an alternative explanation. Finally, we also found an association with the direction of replication, particularly for C>T and T>C mutations ( Figure 4e ), implying that both transcription and replication co-operate to generate these mutation types. 60 Regarding the background process SBS1, this CpG>TpG-associated signature contributes between 10-20% of mutations in DNMs, while its weight in ERVs is reduced to a few percent, highlighting the recurrence of CpG>TpG mutations and their subsequent undercounting in the ERV class. 61 At the same time, all four data sets in Figure 4a show the expected depletion of SBS1 mutations in the 2-4 kbs around the TSS due to demethylation of CpG islands ( Figure 1d-f ), which coincide with the promoter regions of 60-70% of genes. Negative selection and biased gene conversion erode the mutational hotspot in the germline ERVs serve as proxies for germline mutational processes unaffected by selection. When repeating our analysis with higher-frequency variants, we found that the TSS hotspot of the evolutionarily younger mutations is not reflected in older polymorphisms, consistent with direct and background selection acting on it, 62 Figure 5a . This also explains why the hotspot was not detected in substitution data. 31 Quantifying direct selection yielded a dN/dS ratio of 0.35 for common exonic SNPs and a significant depletion of non-synonymous ERVs (dN/dS≈0.82), indicating that variants with allele frequencies up to 0.01% (i.e., <14 and <30 copies in gnomAD and UKBB, respectively) were already subject to measurable purifying selection, Figure 5b . In contrast, DNMs showed a slight dN/dS increase above 1, reflecting the overrepresentation of diseased individuals in this cohort. 63 Download figure Open in new tab Figure 5: Impact of selection and biased gene conversion on the TSS mutational hotspot. (a) Regional variability in µ ′ for non(CpG>TpG) mutations in 100-bp bins around the TSS for allele frequencies below 0.01% (ERVs), between 0.01% and 0.1%, between 0.1% and 1%, between 1% and 10%, between 10% and 90% (“common SNPs”), and above 90% (top to bottom, gnomAD). (b) From top to bottom: dN/dS for ERVs, common SNPs and DNMs for non(CpG>TpG); ratio of µ of common SNPs and ERVs for three different filtering criteria (all mappable sequences, only conserved elements and all except conserved elements); density of conserved elements; and mean recombination rate per site around the TSS. Note that a slightly elevated value of µ common / µ ERV on untranscribed sequences upstream of the TSS is likely spurious and due to the genome-wide normalisation of µ common . (c) Regional variability in µ ′ stratified by biased gene conversion (BGC) with favoured, indifferent and unfavoured mutations (left to right) and by allele frequency (top to bottom). Error bars indicate the 90% confidence intervals for the mean across 100 bootstrap replicates (when possible). To assess purifying selection including non-coding sequences, we calculated the ratio of mutation densities in common polymorphisms ( µ common ) and ERVs ( µ ERV ), yielding an average µ common / µ ERV of 0.79 for coding sequences, Figure 5c . The sharper reduction near TSSs compared to distal regions likely reflects the higher density of conserved elements, modulated by the effects of recombination disrupting genetic linkage. GC-biased gene conversion (BGC), which favours C/G over A/T nucleotides, 64 may also contribute to hotspot erosion. The TSS hotspot is strongest for mutations indifferent to BGC (A<>T, C<>G, 48% excess) or favoured by BGC (AT>CG, 39%), and weakest for BGC-unfavoured mutations (CG>AT, 24%), Figure 5d . Consistently, BGC-unfavoured mutations showed the greatest reduction at TSSs among common SNPs, supporting the roles of purifying selection and BGC in hotspot erosion during evolution. Gene expression correlates positively with mutation density in the soma and the germline when accounting for other factors Finally, we used our genome-wide window-based negative binomial regression model, controlling for 38 possible covariates of mutation rate, to test the global effects of expression on mutation density. A recent study found an overall negative correlation between transcription and germline mutagenesis. 9 , 10 However, a later re-analysis came to opposite conclusions, 11 which in turn are consistent with earlier observations. 12 , 13 In the soma, a significant negative correlation between transcription level and mutation rate is usually assumed. 2 Here, using the expression data from Xia et al . (2020) 9 we found a significant positive association between µ tb and gene expression for non(CpG>TpG) on transcribed sequences in both ERV data sets and a partially significant negative correlation for CpG>TpG ERVs , Extended Data Figures 13-14 . In cancer, replication time is the most important predictor of mutation density in all mutation types and regions ,which is consistent with previous reports. 2 , 13 Interestingly, correcting for this and the other features, we found a positive correlation with gene expression in our non-hypermutable cancer tumours for all mutation types .While this observation may be different for some cancer types when considered in isolation, it suggests that the reported apparent negative effect of gene expression on somatic mutation density is driven by other covariates of mutation rate, mainly replication time 32 , 65 ( Extended Data Figure 10 ). Discussion Transcription and its associated processes are pivotal in shaping mutation rates. In somatic cells, it has been well-established that transcription start sites exhibit elevated mutation rates in many cancers, primarily due to increased DNA damage during transcription factor binding. Our study demonstrates the existence of a distinct mutational hotspot in the human germline, spanning several hundred nucleotides upstream and downstream of the TSS. Despite variations in ancestry among cohorts, we replicated these findings in two independent datasets, comprising 70000 individuals (gnomAD) and approximately 150000 individuals (UKBB). The TSS hotspot is linked to several mechanisms, including mitotic double-strand breaks, RNAP-II stalling, divergent transcription, and R-loops. Mutational signature analysis further supported these mechanisms, highlighting four key mutational processes: (1) T>C mutations linked to transcription (SBS12 and SBS16, known from liver tumors); (2) C>G mutations associated with transcription and complex crossovers in the female germline (SBS39); (3) C>T mutations of unknown origin (SBS11 and others); and (4) defective double-strand break repair (SBS3, SBS40b, and SBS40c, all mutation types). Surprisingly, de novo mutations did not show a significant excess near the TSS. Beyond the smaller sample size and reduced statistical power of the DNM dataset, we attributed this finding to the exclusion of early mosaic variants from family sequencing data. Supporting this, we observed a significant 52% excess (p = 0.005) of early mosaic mutations within the first 1 kb downstream of the TSS, linked to double-strand breaks and RNAP II stalling. Mosaic variants have long been recognized as crucial contributors to germline mutational burden, with lower-bound estimates of their proportion ranging from 4 to 26%. 34 , 35 , 66 , 67 , 68 Our findings suggest that cohort studies relying on de novo mutations for disease associations may miss key regulatory and exonic mosaic variants at the TSS that play a significant role in disease phenotypes. 41 Neglecting the mutational hotspot at the TSS can also distort estimates of the strength of negative selection. A recent study using ERVs to estimate the fraction of ultra-selected variants encountered an unexpected and implausible negative fraction of strongly deleterious variants at the TSS, driven by a mutational model that did not account for the natural excess of mutations in this region. 69 Our findings provide a framework to improve mutational models, enabling more accurate inferences of negative selection in promoters, 5’-UTRs, and exonic regions. 70 , 71 Our results align with evidence of mutational excess at transcription factor binding sites of genes highly expressed in the male germline. 20 They also complement recent studies detailing mutational excess near recombination-associated PRDM9 binding sites 51 and constitutive replication origins, 72 both associated with double-strand breaks and TMEJ. Moreover, recent insights reveal connections between SBS3 and imperfect DSB repair with non-canonical recombination attempts, 38 while meiotic DSB hotspots lacking PRDM9 binding fail to undergo homologue-guided repair. 73 These findings underscore the critical yet underexplored role of non-canonical DSB repair in germline mutagenesis. Our study places this process, alongside transcription, at the centre of TSS mutagenesis, revealing a specific developmental timeline of events. Methods Definition of Genomic Windows Processing of Genomic Positions of Interest We obtained FANTOM5 gene annotations from the ZENBU genome browser ( https://fantom.gsc.riken.jp/zenbu/gLyphs/index.html#config=cZzb9dCAO5Ptk_A-QAjmGC ). 74 , 75 This resource provides robust locations of Transcription Start Sites (TSSs) since annotations were created by associating peaks of Cap Analysis of Gene Expression (CAGE) to 5’ ends of genes annotated across different genomic sources. 76 There are 21069 proteincoding genes, 8187 divergent long non-coding RNAs (lncRNAs), 13105 intergenic lncRNAs, 4510 pseudogenes, and 12239 genes of other types in the FANTOM5 robust gene track. Genomic positions of interest were defined as a maximum of 50 kilo-base-pairs (kb) upstream and 50 kb downstream of both the TSS and Transcription Termination Site (TTS). The analysis was limited to the gene length or the first/last 50 kb if the gene was shorter or longer than 50 kb, respectively. We binned these genomic positions of interest into non-overlapping windows of up to 1 kb or 100 bp, taking as a starting point the TSS or TTS of each gene in the autosomes. Transcription, Mappability, and Liftover Filters We excluded sites from the windows defined above that lay in untranscribed upstream TSS or downstream TTS regions that could actually be transcribed as part of other genes different from the one considered in each case. If untranscribed upstream or downstream positions from different neighbouring genes overlapped, we assigned the position to the nearest TSS or TTS. Hence, no position is represented more than once in our analysis and all positions are labelled as a function of distance to its closest TSS/TTS gene annotation. Sites transcribed as part of more than one gene, regardless of the transcription strand orientation (same or opposite), were excluded from the analysis. We subtracted hard-to-map regions from the genomic windows in favour of using only sites present in uniquely mappable regions, based on the CRG-36mers alignability filter obtained from https://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg19&g=wgEncodeMapabilit y. 77 These regions can promote the calling of artificial mutations and boost the counts of recurrent mutations across samples. 77 Since we worked with some mutation data sets originally provided in the hg38 genome assembly, in order to ensure a comparable window target size of mutation across data sets after using liftover to convert hg38 to hg19 coordinates, 78 we identified the sites in hg19 that can be uniquely mapped to hg38 and for which the identity of the reference nucleotide does not change after lifting over. Sites not meeting these two specifications were subtracted from the genomic windows of interest. After applying the filters described above, we were left with the following genes for downstream analyses involving the estimation of mutation density: 14763 protein-coding genes, 3991 divergent lncRNAs, 8454 intergenic lncRNAs, 1660 pseudogenes and 1479 genes of other types. Supplementary Table 3 contains the identifiers and FANTOM5 coordinates of all these genes together with their corresponding number of windows and sites per relative position to the TSS or TTS. Labeling Putative Regions under Selection Exons from GENCODE version 19 and conserved sites as predicted by GERP scores from transcribed and neighbouring untranscribed regions were identified as putatively selected regions. 79 , 80 We produced versions of the genomic windows where these putatively selected regions were subtracted from the set of analysed sites (windows without conserved elements) or where they were left as the only sites to analyse (windows with only conserved elements). Throughout the manuscript, we note when we have used windows with or without these additional modifications. Mutation data sets gnomAD and UK Biobank Germline variants were collected from gnomAD v3 ( https://gnomad.broadinstitute.org/ ) and from UK Biobank ( www.ukbiobank.ac.uk ). 21 , 81 We used a liftover to convert both gnomAD and UK Biobank data from hg38 to hg19 coordinates. The derived allele frequency spectrum from gnomAD was partitioned into six frequency intervals. Alleles occurring in 1-14 copies (or <0.01%) were considered extremely rare variants, followed by alleles between 15-140 (0.01%-0.1%), 141-1400 (0.1%-1%), and 1401-14000 (1%-10%). Common polymorphisms were defined as alleles occurring between 14001-126000 (10%-90%), and high-frequency polymorphisms were defined as alleles occurring in more than 126000 copies (top 10%). We downsampled the original variable sample size at each site to 140000 chromosomes and discarded sites with a sample size lower than 140000 chromosomes (those represented less than 5% of all segregating sites). For the UK Biobank data, we used approximately 300000 chromosomes and kept only variants with an allele frequency lower than 0.01% (ERVs). There are 414,433,805 gnomAD ERVs, 4,814,941 gnomAD common variants, and 520,004,351 UK Biobank ERVs before restricting to mutations in our genomic windows of interest. Cancer Data Somatic mutations were retrieved from the PCAWG (Pan-Cancer Analysis of Whole Genomes) database. 22 We removed hypermutable tumors/patients of skin, colorectal, and a few other tumor types described in. 52 The identifiers of the patients used in this study are listed in Supplementary Table 4 . There are 22,587,102 mutations across cancer types before restricting to mutations in our genomic windows of interest. De novo and mosaic mutation data sets Mutations from parent-offspring or multi-generation sequencing and somatic mosaic mutations from deep sequencing were collected across 13 studies. Mutations shared between studies or between supplementary tables of the same study were filtered out, such that all data sets contain only unique mutation calls. A total of 10856 early mosaic variants were collated from Supp. Table 3 of Ju et al . (2017), 36 containing low-VAF mutations from blood (163 unique mutations) the file gonosomal.dnms of Sasani et al . (2019) 35 (440 unique mutations) Supp. Table 1 of Jonsson et al . (2021), 33 containing near-constitutional (VAF>45%) mutations from blood that differed between monozygous twins, excluding mutations shared by the twins (874 unique mutations) Supp. Table 2 of Jonsson et al . (2021), 33 containing mutations detected in both soma and germline of one monozygotic twin, the latter evidenced by transmission to an offspring (3842 unique mutations) Supp. Table 3 of Jonsson et al . (2021), 33 containing mutations present in both twins with <25% VAF or significant VAF difference between twins (85 unique mutations) Supp. Table 2 of Rodin et al . (2021), 37 containing predicted non-constitutional SNVs from brain (2166 unique mutations). Supp. Table 3 of Maury et al . (2024), 43 containing SNVs with VAF<0.4 from brain, presumed to have arisen during prenatal neurogenesis (3286 unique mutations). A total of 4764 late mosaic variants, i.e., variants classified as de novo variants in the following studies after filtering out variants present in parental blood, but shared between siblings, were collated from Supp. Table 15 of Goldmann et al . (2016) 39 (26 unique mutations) Supp. Table 3 of Yuen et al . (2017) 40 (2555 unique mutations) Supp. Table 4 of Jonsson et al . (2017) 34 (84 unique mutations after removal of 207 variants shared with Halldorsson et al ., 2019) Supp. Table 2 of An et al . (2018) 41 (1257 unique mutations) Goldmann et al . (2018), 38 obtained directly from the authors, but also available at dbGaP accession phs001522.v1.p1 (9 unique mutations). the file post-pgcs.dnms of Sasani et al . (2019) 35 (303 unique mutations) the file aau1043 datas5 revision1.tsv of Halldorsson et al . (2019) 42 (455 unique mutations) Supp. Table 2 of Jonsson et al . (2021) 33 (75 unique mutations) A total of 774930 de novo variants (“DNMs”), i.e., variants identified in the following studies that were not shared between siblings, were collated from www.nlgenome.nl of Francioli et al . (2015) 82 (11016 unique mutations) Supp. Table 3 of Goldmann et al . (2016) 39 (35730 unique mutations) Supp. Table 3 of Yuen et al . (2017) 40 (113960 unique mutations) Supp. Table 4 of Jonsson et al . (2017) 34 (9468 unique mutations, after removal of 88721 variants shared with Halldorsson et al ., 2019) Supp. Table 2 of An et al . (2018) 41 (231362 unique mutations) the file aau1043_datas5_revision1.tsv of Halldorsson et al . (2019) 42 (180287 unique mutations) Goldmann et al . (2018), 38 obtained directly from the authors, but also available at dbGaP accession phs001522.v1.p1 (73737 unique mutations) the files second_gen.dnms and third gen.dnms of Sasani et al . (2019) 35 (27684 unique mutations) Supp. Tables 2 and 3 of Richter et al . (2020) 83 (72096 unique mutations, after removal of 88279 variants shared with An et al ., 2018) Supp. Table 2 of Jonsson et al . (2021) 33 (19590 unique mutations, after removal of 369 variants shared with Halldorsson et al ., 2019). Mutation Density across Genomic Windows Computation of µ and µ ′ The observable µ measures the observed mutation probability of a given genomic locus, corrected for the sequence composition of the locus. It is equivalent to the observed number of mutations in a locus divided by the number expected under a uniform distribution. The expected number of mutations is derived from the mutational matrix which represents the global average pentanucleotide-sequence-context dependence of all active mutational processes. That is, the matrix element m ij quantifies the probability of the central base in the pentanucleotide of type i = 1,. .., 1024 (corresponding to AAAAA, AAAAC, .. . , TTTTT) on the coding (i.e., non-transcribed) strand to mutate to state j = 1,. .., 4 (corresponding to A, C, G and T): where is the number of mutations of type ( i, j ) at site l = 1,. .., L and L is the total sequence length. The abundance of the pentanucleotide context i among all L sites is given by ,with denoting the indicator function and c l the pentanucleotide at site l, such that .The summation is taken over all L genomic sites considered in a given analysis. By default, this includes all sites in the (up to) 50 kb untranscribed sequence upstream of the TSS, the (up to) 50 kb transcribed sequence downstream of the TSS, the (up to) 50 kb transcribed sequence upstream of the TTS and the (up to) 50 kb untranscribed sequence downstream of the TTS of each gene. All 4 gene categories (14763 protein-coding, 3991 divergent lncRNAs, 8454 intergenic lncRNAs, and 1660 pseudogenes) were processed together. When a gene was less than 100 kb in length, double-counting of mutations on the transcribed region was avoided by assigning each mutation to whichever element it was closer to out of TSS and TTS. When conserved sequences were excluded from the analysis, the matrix was correspondingly recalculated from only those parts of the sequence that remained. When using a window size of 100 bp (instead of 1 kb), the considered regions were correspondingly reduced to 5 kb (instead of 50 kb). For each data set (UKB ERVs, gnomAD ERVs, gnomAD common variants, DNMs, and pan-cancer), a separate mutation probability matrix was derived. To compute µ tb in a given genomic window b of a given gene t , we divided the observed mutation count in the window by the expected, as computed from the matrix: where W ∈ {100, 1000} is the size of the window. In other words, all observed mutations in the window were counted and this count was divided by the number of mutations expected in the window based on the genome-wide average mutation probabilities. We then took the average of µ tb across all genes to arrive at our observable µ b ( Extended Data Figure 9 ). Since we separated mutations into two categories (non(CpG>TpG) and CpG>TpG), removed conserved sequences (including exons) for some analyses and a gene’s TTS may not coincide with the 3’ boundary of the window it is located in, we accounted for the increased amount of noise in µ tb coming from genes that have fewer than W sites in a given window by weighting each value with the fraction of sites considered coming from the gene, w tb : To increase the information content of the figures, we also derived a second observable, µ ′ . Beyond the pentanucleotide sequence context, µ ′ also corrects for the regional background mutation probability of the genomic locus in which the gene is embedded since the latter correlates with gene length and regional gene density. To this end, for each window in a given gene t , we divide that window’s µ tb by the ratio g t of the observed and expected mutation count of the entire gene, where ℓ t denotes the total length of regions considered from the given gene, including up-stream TSS untranscribed, downstream TSS transcribed, upstream TTS transcribed, and downstream TTS untranscribed. Error Bars for µ and µ ′ To estimate the sampling variance around our estimates of µ and µ ′ , we derived simple percentile confidence intervals from bootstrap samples. 84 Briefly, for a given gene type, we sampled with replacement as many genes as analysed together with their corresponding observed mutations and repeated this 100 times. µ and µ ′ were calculated after each resampling iteration as described in the previous section. We found the values representing the percentile 5 and 95 of the resulting bootstrap distributions for each window and defined the 90% confidence interval as the range of values between them. Error bars in Figure 1 and in similar plots represent these intervals and points are the means of the bootstrap distributions. Mutation Density Estimates for Mosaic Variants To study the mosaic mutation density at different genomic windows, we introduced a variant of the analysis consisting of estimating a mutational matrix with only 12 mutation types representing the possible mononucleotide substitutions without considering sequence contexts across sites. This matrix is used in place of the pentanucleotide context matrix described in previous sections to compute the expected number of mutations per window. We do this in the context of Figure 3e due to the low number of mutations in mosaic data sets which would result in an extremely sparse pentanucleotide matrix. The glm function from the R stats package with family = Poisson and all other parameters left as defaults was used to automate the estimation of mutation density and significance testing via the summary method. 85 An interceptless model was fit with the mutation count as the response variable, the window position relative to the TSS as the only explanatory variable, and the logarithm of the expected mutations as an offset. The resulting coefficient for each set of windows is equivalent to the logarithm of the sum of observed mutations across the windows of the set divided by the corresponding sum of expected mutations. A 2-sided z-test with the null hypothesis that the model’s coefficients are equal to 0 was performed (p-values reported in the main text). 95% confidence intervals were derived assuming normality of the sampling distribution of the coefficients using the confint method. Coefficients and interval limits were exponentiated back to the original scale of the observed-to-expected fold change of mutations for interpretability purposes in Figure 3 . Signature Decomposition Analyses To identify genomic window-specific candidate mutational processes contributing to the observed mutations, we ran a signature decomposition analysis using non-negative least squares (NNLS) optimisation as implemented in the nnls R package. 86 , 87 Pooled mutation counts across genes were calculated according to the positions of the windows of origin of each mutation relative to the TSS. We leveraged two different sets of signatures: COSMIC v3.4 56 and the germline components from Seplyarskiy et al . 24 COSMIC mutational signatures were extracted using whole-genome (WG) cancer data and were not normalised by the trinucleotide abundance of the genome. Since our analysis is restricted to genomic windows, we need to re-scale our data to account for their different trinucleotide composition. To do so, we divided each trinucleotide mutation type by its window-specific abundance and multiplied it by the WG abundances. This re-scaled input is the one we used for the NNLS decomposition. Signatures with an assigned weight lower than a cutoff (between 0.01 and 0.1, depending on the study and the bin) were not plotted, and their weights were aggregated into a category called “Others”. Since genomic windows with very few mutations cannot be reliably decomposed using NNLS, we restrict the decomposition to windows with at least 1000 mutations. Moreover, in order to account for the effect of the number of mutations on the quality of the decomposition when comparing different windows, we downsample the number of mutations in each window to the number of mutations in the window with the lowest count. Finally, to compare the mutation load associated with the signatures across windows, we multiply the signature weights of each window by the total number of mutations in the window and divide by the total number of analysed sites in the window, giving the mutation density. The germline components were extracted from extremely rare variants. 24 These components were normalised by the trinucleotide abundance and are available in two different formats: extracted from data normalised by the standard deviation and without normalising. We used the non-normalised version for our analyses. Moreover, these signatures have 192 mutation categories, compared to the standard 96 used in the COSMIC signatures. This is because mutation types were not collapsed to the pyrimidine bases in order to analyse strand biases. Therefore, when decomposing samples with these components, we do not collapse the mutation types and we normalise the input by the trinucleotide abundance of the windows. To compare the component weights in ERVs and DNMs, we pooled mutations from 20 kb around the TSS (10 kb to either side) and decomposed all these mutations together, for each data set separately. Genomic Features Below are listed the 38 genomic features used throughout this work together with their sources. We aimed to assign one value or label per feature to each of the genomic windows of up to 1 kb (only with the transcription filter defined in the “Transcription, Mappability and Liftover Filters” section). When required, we used the Big Binary Indexed suite of tools to convert files to BED format. 88 Other tools used for data parsing and genomic coordinate intersection were wig2bed and bedtools . 89 , 90 When necessary, we lifted genomic coordinates from hg38 to hg19 using liftover . 78 Most of the source files we use can be accessed through the University of California Santa Cruz (UCSC) genome browser or Gene Expression Omnibus (GEO) resources. 91 , 92 Unless otherwise specified, for non-categorical features based on various signal intensities, we intersected the reported regions in the source files with our genomic windows and computed the per-base signal density in each window assuming a signal of 0 when a region was not reported in the files. Features denoting properties of genes had their values or labels propagated to all windows associated with a given gene. Huvec and H1hesc histone marks: we averaged the signal files across both cell lines of 11 different histone marks measured by the EN-CODE project. 93 These include H3K9me3, H3K4me1, H3K4me2, H3K4me3, H3K9ac, H3K9me1, H3K27ac, H3K27me3, H3K36me3, H3K79me2, and H4K20me1. Signal file contents and download locations are documented at http://genome.ucsc.edu/cgi-bin/hgFileUi?db=hg19&g=wgEncodeBroadHistone . CpG methylation: since we found poor methylome coverage in somatic tissue-specific data sets such as ENCODE (data not shown), we used instead the high-coverage methy-lome from the iMETHYL database ( http://imethyl.iwate-megabank.org/downloads.html ). 94 The methylome of human blood cells (monocytes, CD4T, and neutrophils) was downloaded from the iMETHYL database and used as a proxy for the methylation status of CpG dinucleotides in the human genome across tissue types. DNase I Hypersensitivity Sites (DHS): DHS files were downloaded for the tissue of origin of each tumor type. As a proxy for germline DHS, we use the H1 human embryonic stem cell line. Data originate from the ENCODE project and are deposited as the file wgEncodeRegDnaseClusteredV3.bed.gz at http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeRegDnaseClustered/?C=M;O=D . 93 To calculate the pan-cancer DHS signal, we took the weighted mean signal across tissues of origin according to the number of tumors from each tissue type in PCAWG (see Supplementary Table 5 ). Genomic Evolutionary Rate Profiling (GERP) score: conserved elements according to the GERP score were retrieved from the Sidow Lab at Stanford University ( http://mendel.stanford.edu/sidowlab/downloads/gerp/index.html ). 80 Somatic double-strand breaks (DSB): we obtained a high-sensitivity single-nucleotide resolution somatic DSB map for NHEK cells from the work of Lensing et al ( https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=G SE7817 2, file BREAK primary n1.NOdups.nomdel. default_peaks.narrowPeak.gz). 44 Huvec cells H2A.Z histone mark: we downloaded the ENCODE signal file documented at https://genome.ucsc.edu/cgi-bin/hgTables?db=hg19&hgta_group=regulation&hgta_track=wgEncodeBroadHistone&hgta_table=wgEncodeBroadHistoneHuvecH2azSig&hgta_doSchema=describe+table+schema . 93 Nucleosome density: we retrieved high-resolution lymphoblastoid B-cell MNase-seq stable nucleosome regions from Gaffney et al . 95 We focused on the normalised nucleosome occupancy signal provided by NucPosDB (ht tp s: //g en er egu la ti on. or g/N GS /s tab le _n ucs /h g1 9/, file GSE36979_Gaffney2012_Bcells_MNase-seq_stable_100bp_hg19.bed.gz). 96 Meiotic recombination hotspots associated with PRDM9: we retrieved a signal of human testis meiotic double-strand break hotspots found through their relationship with the binding of DMC1 at PRDM9 binding motifs from Supp. Table 1 of Hinch et al . 51 We extended the coordinates of the hotspots 100 bp in both the 5’ and 3’ directions to capture the main peak of the mutation footprint around hotspots that is reported by Hinch et al . 51 The signal intensity reported in the data was used for all 201 bp associated with each hotspot after extending the annotations. Replication time: we computed the mean wavelet-smoothed signal across 6 cell lines measured by the ENCODE project. 93 The cell lines are: Gm12878, Helas3, Hepg2, Huvec, K562, and Mcf7. The data files and download locations are documented at https://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg19&g=wgEncodeUwRepliSeq . We intersected the reported regions with our genomic windows and computed the perbase signal density in each window. If a base had no reported replication time, the base was ignored altogether during the calculation. Recombination rate: the average genetic map of recombination rates computed from the paternal and maternal maps was retrieved from the results of Kong et al . documented at https://genome.ucsc.edu/cgi-bin/hgTables?db=hg19&hgta_group=map&hgta_track=decodeRmap&hgta_table=decodeSexAveraged&hgta_doSchema=describe+table+schema . 97 The per-base density of this feature was calculated for each genomic window similarly to that of replication time. Replication strand: we took the data from the Mcf7 cell line used for calculating the replication time feature and leveraged it to compute a replication direction map using the reptDir R package. 98 The method used by the package is related to the descriptions provided by previous work on replication direction, 99 , 8 and is documented in the cited package. We specified a minimum replication direction domain length of 250 kb and a minimum slope of 0.01 between individual regions in the domains. Then, we used the obtained direction map to determine if each genomic window was replicated in a single direction, in a mix of both, or if it had base pairs for which the direction was unreliable. For windows replicated in a single direction, we were able to determine the identities of the leading and lagging strands of replication. As mutations were counted in the coding strand of transcription, we created a categorical variable with 3 possible labels indicating if the coding strand is also the leading strand, the lagging strand, or if such thing is unknown as is the case for windows with base pairs replicated in an unreliable direction or with a mix of base pairs replicated in both directions. G/C skew: we calculated the difference between the number of guanines and cytosines for each genomic window and divided it by the sum of both nucleotide counts. The computation is applied to the template strand of transcription according to the annotations of the particular gene associated with each window. CpG content: derived from the number of CpG sites in the genome ( https://doi.org/10.6084/m9.figshare.1415416.v1 ). 100 This is the count of sites that are found in each genomic window divided over the total bases in the window. G/C content: this is the simple sum of the number of guanines and cytosines in each genomic window divided by the number of total bases in the window. CpG island (CpGI) content: we defined CpGIs using the annotations file available at http://genome.ucsc.edu/cgi-bin/hgTrackUi?hgsid=1102630601_okTYJ1lN4a4kGQ6CPBEtw3g218bg&c=chr1&g=cpgIslandExtUnmasked . 101 We intersected the coordinates with our genomic windows and computed the fraction of CpGI bases in each window. Density of G-quadruplexes: genomic coordinates of G-quadruplexes in the human genome with experimental evidence were downloaded from the EndoQuad database ( http://chenzxlab.hzau.edu.cn/EndoQuad/#/download , file Human_eG4.txt). 102 We computed the number of bases that intersected each genomic window and divided it by the total number of bases in the window. Gene with CpGI promoter: promoters of each gene were defined as regions ranging from -1.2 kb to +0.3 kb relative to the TSS of BioMart Ensembl75 annotations. 103 The CpGI annotations cited for the previous feature were intersected with these promoters. A binary categorical variable was defined depending on whether or not at least 1 bp of the promoter intersected with a CpGI. Gene TATA box score: we retrieved the position frequency matrix for the TATA box motif from the JASPAR database ( https://jaspar.elixir.no/matrix/MA0108.3/ ). 104 Each matrix element was divided by the sum of elements in the column where it is located. The result was further divided by 0.25 which represents the background frequency for all 4 nucleotide types. Then, to obtain the Position Weight Matrix (PWM), 105 a pseudocount of 1 was added to all elements before taking their base 2 logarithms. We obtained the sequences in the genome found in sliding windows of the same size as the TATA box motif between -40 bp and -20 bp of the TSS of each gene. 106 The PWN was used to score each sliding window and their reverse complements. Finally, the TATA box score is the maximum score across all sliding windows of each gene. Gene with upstream divergent promoter: paired TSS annotations derived from Global Run-on cap (GRO-cap) sequencing data obtained from k562 cells were retrieved from the Supp. Data set of Core et al . 47 Each pair annotates a GRO-cap TSS in the plus strand that is located 150 bp or less from another GRO-cap TSS in the minus strand. We defined the putative divergent promoter region as the gap between the two GRO-cap TSS defined in each pair. In a minority of cases, the coordinates of the two GRO-cap TSS overlapped and we took the putative divergent promoter as the overlapping region. We then defined a binary categorical variable depending on whether or not putative divergent promoters overlapped with the first genomic window upstream of the TSS belonging to each gene. Gene with downstream divergent promoter: this feature is similar to the one defined above with the only difference being that putative divergent promoters are intersected with the first bin downstream of the TSS. Gene length: corresponds to the total length of the whole gene of origin of each genomic window. The lengths include intronic regions and are computed directly from the FANTOM gene models. 75 Gene expression: as a proxy of gene expression in the human germline, and given that more than 75% of DNMs are from paternal origin, 107 we used single-cell transcriptomic unique molecular identifier counts from male germline cells provided upon request by Xia et al . 9 For somatic expression, we considered data from the Genotype-Tissue Expression (GTEx) project version 4 ( https://dcc.icgc.org/releases/PCAWG/transcriptome/transcript_expression/ , file GTEX_v4.pcawg.transcripts.tpm.tsv.gz ). 108 , 109 We averaged the expression of each gene across samples of the same tissue using the metadata file GTEX_v4.metadata.tsv.gz available at https://dcc.icgc.org/releases/PCAWG/transcriptome/metadata . 109 Supplementary Table 6 shows the number of samples averaged for each tissue type and the correspondence with PCAWG cancer types. To obtain the pan-cancer expression, we applied the same procedure as described for the DHS feature using the mean expression across tissues as input for the weighted average. Gene is essential gene: essential genes found by Wang et al . and processed according to Weghorn and Sunyaev were retrieved from the supplementary files of the latter authors. 110 , 111 We defined a binary categorical variable denoting whether or not a gene is in this set of essential genes. Gene is Tumor Suppressor Gene (TSG): TSGs without a dual role as oncogenes were extracted from the Cancer Gene Census previously available at https://cancer.sanger.ac.uk/census 112 (560 genes total before filtering). Both tiers 1 and 2 were considered. We defined a binary categorical variable denoting whether or not a gene is in this set of TSGs. Gene is oncogene: similar to the previous feature but for oncogenes without a dual role as TSGs. Gene Precision Run-on Sequencing (PRO-seq) ratio: we obtained strand-specific maps of single-nucleotide resolution engagement of RNA polymerase II complexes in K562 cells from the file GSE60456_RAW.tar at https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE60456 . 113 The reads of these PRO-seq data originating from the transcribed strand of each gene were identified. We summed the reads counted in the first 500 bp of the first genomic window downstream of the TSS of each gene. We did the same for the reads counted in the remainder of the aforementioned window and the rest of the downstream TSS windows. We then divided these sums by the appropriate sum of sites to obtain the per-site PRO-seq signal density at the beginning and at the gene body of each gene. We added a pseudocount of 1 to these values before computing the difference between the base 10 logarithms of the beginning and gene body signals. The result is the final PRO-seq ratio feature for each gene. Gene from maternal, non-maternal, or other chromosomes: a categorical variable with three possible labels denoting the chromosome group where the gene is found was specified. The definition of the groups was guided by the results of Seplyarskiy et al . who showcase a mutational process acting on the human germline that is enriched with clustered mutations of maternal origin. 24 We defined the chromosome groups according to the fraction of maternal mutation clusters per chromosome that are shown in the cited work. We focused on the four most extreme chromosomes on both sides of the spectrum. Groups consist of maternal chromosomes (8, 9, 15, and 16), non-maternal chromosomes (1, 13, 18, and 20), and other chromosomes. Regression Analyses Modelling of Mutation Density We model the mutation density across genomic windows in the data sets of interest as a variable dependent on predictors derived from genomic features. This is done within a Generalised Linear Model (GLM) framework using a log link function: 114 Where y is a n -dimensional vector containing mutation counts for the different genomic windows (the numerator of Equation 2 ). X is a n -by-p-dimensional design matrix, β is a p -dimensional vector of regression coefficients which are to be estimated from the data and u is a n -dimensional vector of exposures (denominator of Equation 2 ). The ln( u ) term acts as an offset that adjusts the mutation counts exactly by the expected number of mutations in the corresponding window. 114 Under the assumption that the conditional errors around the model’s predictions are Poisson distributed, we estimated the regression coefficients using the glm function from the stats R package with family = poisson or using the glm.nb function from the MASS R package (Poisson and negative binomial regression, respectively). 85 , 115 Negative binomial regression is useful in the presence of overdispersion as it introduces additional variance by allowing extra gamma-distributed errors dependent on a single additional parameter that can be estimated together with the regression coefficients. 115 For all non(CpG>TpG) mutation data sets, negative binomial regression was used. For CpG>TpG mutation data sets, attempts to fit a negative binomial regression model tended to not converge due to a lack of overdispersion in these cases. 116 Therefore, Poisson regression was applied to these data sets. Parameters for the estimation routines were left as defaults except for the number of maximum iterations, which was increased to 1000. Design Matrix and Regression Coeffcient Interpretation In the context of the regression results shown in Figure 3 , we aimed to describe changes in mutation density across genomic windows as a function of both position-unspecific changes (main effects) and position-specific changes in genomic features (interaction effects). 117 To achieve this, the design matrix was defined with matrix augmentation as: Where 1 T is a column-vector of ones. A is a matrix composed of m = 1, …, M ( M = number of features) column-vectors containing the genomic feature information of all windows such that A = [ x 1 x 2 .. . x m ]. B is a matrix of k = 1, …, K ( K = number of position categories) column-vectors encoding the position category membership of any particular window such that B = [b 1 b 2 .. . b k ], and C is a matrix containing all Hadamard products between pairwise combinations of column-vectors in A and B such that C = [( x 1 ⊙ b 1 ) ( x 1 ⊙ b 2 ) .. . (x m ⊙ b k )]. The number of column-vectors in X is therefore p = 1 + m + k + mk . Categorical variables were included in A through dummy coding and position categories in B through sum coding. 118 Overall, this design allows us to make the following interpretations: Main effects: they are represented by the regression coefficients .These are the average changes in the log mutation density across position categories per unit change in a non-categorical genomic feature. For categorical features, the coefficient represents a similar change when switching between the reference and non-reference feature labels. For almost all of these features, the reference is “no” and the only coefficient associated with each feature reflects the change to “yes”. The only exceptions are maternal chromosome and replication strand status, where there are two possible non-reference labels and hence two coefficients for each of these features instead of one. The reference labels for these features are “other” and “unknown”, respectively (see the “Genomic Features” section). Interaction effects: they are comprised by regression coefficients . These are additional changes in log mutation density on top of the main effects that also depend on the set of features in A but only apply in the context of specific position categories. In typical parametric regression fashion, the interpretation of coefficients associated with a particular genomic feature is done considering all other features that have been adjusted for (i.e., held at some constant value). Note that by taking any coefficient and computing 100 × (1 − e β ), one can obtain the corresponding percentage change in the original unlogged mutation density scale. For regressions in Extended Data Figures 13-14 , we used a similar approach as described above but the C matrix component of the design matrix was omitted (i.e., models with no interactions), and the matrix B was instead reworked into a simple integer vector denoting distance from the TSS (predictor denoted as “bin”). For regressions in Figure 3f , including all genomic features in X causes convergence problems while fitting the models, due to the low number of mutations in some of the data sets considered in this analysis. Therefore, in this case, we opted to fit one regression per feature keeping the design with main and interaction effects. To simplify the coefficient interpretation of the numerous non-categorical features that we include in each regression, each with different units of measurement, standardised regression coefficients are provided for all analyses. 119 To achieve this, we simply mean-centered and variance-scaled the non-categorical predictor column vectors in A after filtering out windows with missing data for any of the features and before calculating C when this matrix was needed. Addressing Multicollinearity and Model Diagnosis For all regressions that considered more than one genomic feature at a time, we performed two procedures to rule out major problems confounding the estimation of regression coefficients and their associated inference procedures: Multicollinearity is associated with unstable estimation of coefficients and inflated standard errors that can obscure model interpretation and inference. 120 To prevent strong multicollinearity, we applied the following procedure: with the whole set of 38 features, a simplified model with no interactions was fitted ( C in Equation 6 omitted from the design matrix). Then, the Generalised Variance Inflation Factor (GVIF) was calculated for each feature using the car R package. 121 To make non-categorical and categorical features comparable, we took the GVIF to the power of the reciprocal of the degrees of freedom associated with each feature as the multicollinearity estimator. 122 The feature with the maximum value of this estimator was identified and if its value exceeded a predefined threshold, it was removed and a new reduced model was fitted. We repeated this procedure iteratively until no features exceeded the threshold. With the set of remaining features, we fitted the final models for each data set including at this point the C matrix in the design if necessary. We worked with a threshold of 5 as a middle ground between conservative values and more lax criteria. 123 Model mispecification can manifest as severe overdispersion or underdispersion of the data compared to the model’s predictions. 124 Particularly, overdispersion leads to underestimation of coefficient standard errors and therefore false positives when performing inferences. 125 After fitting each regression, we examined the residual dispersion statistic provided by the DHARMa R package. 126 We found that in most regressions, we observe only slight deviations from equidispersion. In the cases where stronger deviations exist, they trend towards underdispersion which has the effect of making the inference more conservative. 124 In Supplementary Tables 1-2 , we provide the list of genomic features that were filtered out due to the GVIF threshold for each regression shown. We also provide the dispersion statistics along with other summary statistics such as the number of mutations and number of windows involved. For all reported regression results, we ensured there were no warnings or errors raised by the algorithms estimating the regression coefficients. Regression Coeffcient Hypothesis Testing and Inference The aim of the regression analyses presented in Figure 3 is to compare the position-specific effects of genomic features on mutation density within 1 kb from the TSS and at distances 6 kb or farther from the TSS. As a result, vectors of matrix B in Equation 6 encode window membership to four position categories: downstream TSS, downstream far, upstream TSS, and upstream far. This entails that each predictor derived from the genomic features is associated with four position category interactions. After estimation of the regression coefficients, for each predictor, we computed the difference between TSS and far interaction effects associated with the predictor for downstream and upstream positions. This corresponds to comparing the effect on log mutation density that is exclusively attributable to changing the predictor at the TSS against the analogous effect far from the TSS. To complement these statistics, we also calculated the total effect of changing each predictor at the TSS by summing the corresponding interaction effect with the main effect of the predictor. For assessing the significance of these linear combinations of regression coefficients, we assumed normality of the sampling distribution of the estimates to leverage the multcomp R package. 48 The package’s routines allow for computing p-values based on z-tests under the null hypothesis that arbitrary sums or differences of coefficients equal 0 while taking into account the covariance structure between coefficients. Furthermore, p-values and confidence intervals can be adjusted for multiple testing to control the family-wise error rate within each regression. 48 We performed the adjustment using the default parameters of the summary and confint methods for glht objects except for test = adjusted(maxpts = 10 * 50000) , which was needed for successful calculation of p-values when testing many simultaneous hypotheses. For regressions in Extended Data Figures 13-14 , we focused on all 1-kb windows found within 50 kb upstream or downstream from the TSS. For each mutation data set, we fitted one separate regression for each of upstream and downstream windows without interactions. To assess the multiple testing-controlled significance of each predictor, we applied the same approach described above to test the null hypotheses that individual regression coefficients are equal to 0. Mutation Density Strand Bias across Genomic Windows We defined different mutation categories based on the type of mononucleotide substitution (A>G, T>G, A>T, G>T, C>G and C>T), on the strand where the mutation is counted relative to transcription (coding and template), and based on the replication strand annotations described in the “Genomic Features” section. Mutations of each type were summed across all genes according to the genomic window in which they appear. For each window, the resulting mutation counts were divided by the corresponding sum of the number of nucleotides that can lead to each type of mutation. The resulting mutation densities are grouped by mononucleotide substitution type and displayed in panels of Figure 3e (gnomAD data), Extended Data Figure 17 (gnomAD data), and Extended Data Figure 16 (PCAWG data). Evolution Analyses dN/dS Ratio across Genomic Windows In order to have a reading frame of reference for each gene, we used coding regions from GENCODE version 19 whose length is a multiple of 3 to determine the effect of mutations on amino acids. 79 For simplicity, the transcript used for each gene was simply the one specified in the gene expression data we used (see the “Genomic Features” section). Borrowing from the definitions of Equation 2 , we denote the number of mutations of type x in transcript t and genomic window b as where δ x ( l, j ) is a function that evaluates to 1 when the mutation to state j at site l produces a variant of type x in transcript t and to 0 otherwise. The two types of possible mutations x ∈ { z, s } are non-synonymous and synonymous, respectively. Similarly, we define the expected number of mutations of each type as . To obtain the dN/dS ratio for window b we evaluate the function: Potential sampling variance around the ratios was addressed through resampling. Briefly, for each transcript and window, we randomly sampled with replacement as many mutations as observed 100 times. For every resampling iteration, we calculated the dN/dS ratios as described above for each window. Similarly to what was described in the “Error Bars for µ and µ ′ ” section, resamples were used to calculate dN/dS confidence intervals and data points shown in Figure 5b (gnomAD data). For DNMs, due to low numbers of synonymous mutations in some windows, sometimes resampling of mutations would yield 0 synonymous mutations. In these cases, the dN/dS ratio cannot be computed. In windows where this happens, we plot only the observed dN/dS ratio estimate without error bars. Relative Comparison of Mean Mutation Density across Genes for ERVs and Common Variants at Different Genomic Windows We computed µ tb as defined in Equation 2 for the gnomAD ERVs and common variants. Then, we calculated a simple average of µ tb across genes for each genomic bin to obtain µ ERV and µ common , respectively. At each bin position, the ratio µ common / µ ERV was defined as a statistic to measure the relative change in average mutation density across genes due to the action of selection. Possible sampling variance of this statistic is addressed through resampling. We sampled with replacement as many mutations as observed across all genes and windows for each data set 100 times. We calculated µ common / µ ERV as explained above for each window in each resampling iteration. Similarly to what was described in the “Error Bars for µ and µ ′ ” section, resamples were used to calculate µ common / µ ERV confidence intervals and data points shown in Figure 5b (gnomAD data). Effect of GC-biased Gene Conversion on Mutation Density across Genomic Windows For each gene and genomic window, we calculated a value similar to µ tb from Equation 2 but that only takes into account subsets of mutation types. We defined 3 subsets according to the expected action of GC-biased Gene Conversion (BGC) on mutations: A<>T/C<>G substitutions (BGC-indifferent), AT>GC substitutions (BGC-favoured), and GC>AT (BGC-unfavoured). We continued to calculate the analogous version of µ b from Equation 3 for each of the mutation subsets keeping weights proportional to the total number of sites in each window. Then, we also obtained the similarly modified version of g t from Equation 4 in order to compute the mutation subset-specific of each of the 3 aforementioned mutation categories. Error bars and data points for these statistics are plotted in Figure 5c (gnomAD data) and were derived in the same way as described in the “Error Bars for µ and µ ′ ” section. Code availability The code we developed for the analyses presented in this manuscript are available at https://github.com/weghornlab/TSS_hypermutability . Author contributions Conceptualisation: D.C., D.W., data curation: M.C.G., D.C., C.S.C., D.W., formal analysis: M.C.G., D.C., C.S.C., V.S., D.W., methodology: M.C.G., D.C., C.S.C., V.S., D.W., supervision: D.W., writing: M.C.G., D.C., C.S.C., V.S., D.W. Funding We acknowledge support of the Spanish Ministry of Science and Innovation through the Centro de Excelencia Severo Ochoa (CEX2020-001049-S, MCIN/AEI /10.13039/501100011033), and the Generalitat de Catalunya through the CERCA programme. David Castellano has received funding from the European Union’s Horizon 2020 MSCA postdoctoral COFUND programme under grant agreement INTREPiD No. 754422. This work was also partially funded by the Spanish Ministry of Science and Innovation through grants PGC2018-100941-A-I00 and PID2021-128976NB-I00. Competing interests The authors declare no competing interests. Additional information Supplementary Information is available for this paper. Correspondence and requests for materials should be addressed to dweghorn{at}crg.eu . Acknowledgments We thank the members of the Weghorn Group for helpful discussions and suggestions, especially Miguel Rodriguez Galindo. We thank Bo Xia and Itai Yanai for providing access to the curated male germline expression patterns. This publication used data from the UK Biobank under the application ID 81963. We thank Silvia Bonas Guarch and Federico Billeci for their support with UK Biobank data analysis. References ↵ Schuster-Böckler , B. & Lehner , B. Chromatin organization is a major influence on regional mutation rates in human cancer cells . Nature 488 , 504 – 507 ( 2012 ). OpenUrl CrossRef PubMed Web of Science ↵ Lawrence , M. S. et al. Mutational heterogeneity in cancer and the search for new cancer-associated genes . Nature 499 , 214 – 218 ( 2013 ). OpenUrl CrossRef PubMed Web of Science ↵ Alexandrov , L. B. et al. Signatures of mutational processes in human cancer . Nature 500 , 415 ( 2013 ). OpenUrl CrossRef PubMed Web of Science ↵ Moore , L. et al. The mutational landscape of human somatic and germline cells . Nature 597 , 381 – 386 ( 2021 ). OpenUrl CrossRef PubMed ↵ Jinks-Robertson , S. & Bhagwat , A. S. Transcription-associated mutagenesis . Annual review of genetics 48 , 341 – 359 ( 2014 ). OpenUrl CrossRef PubMed ↵ Duan , M. , Speer , R. M. , Ulibarri , J. , Liu , K. J. & Mao , P. Transcription-coupled nucleotide excision repair: New insights revealed by genomic approaches . DNA repair 103 , 103126 ( 2021 ). OpenUrl CrossRef PubMed ↵ Thornlow , B. P. et al. Transfer rna genes experience exceptionally elevated mutation rates . Proceedings of the National Academy of Sciences 115 , 8996 – 9001 ( 2018 ). OpenUrl Abstract / FREE Full Text ↵ Haradhvala , N. J. et al. Mutational strand asymmetries in cancer genomes reveal mechanisms of dna damage and repair . Cell 164 , 538 – 549 ( 2016 ). OpenUrl CrossRef PubMed ↵ Xia , B. et al. Widespread transcriptional scanning in the testis modulates gene evolution rates . Cell 180 , 248 – 262 ( 2020 ). OpenUrl CrossRef PubMed ↵ Xia , B. & Yanai , I. Gene expression levels modulate germline mutation rates through the compound effects of transcription-coupled repair and damage . Human Genetics 141 , 1211 – 1222 ( 2022 ). OpenUrl CrossRef PubMed ↵ Liu , H. & Zhang , J. Higher germline mutagenesis of genes with stronger testis expressions refutes the transcriptional scanning hypothesis . Molecular Biology and Evolution 37 , 3225 – 3231 ( 2020 ). OpenUrl CrossRef PubMed ↵ Park , C. , Qian , W. & Zhang , J. Genomic evidence for elevated mutation rates in highly expressed genes . EMBO reports 13 , 1123 – 1129 ( 2012 ). OpenUrl Abstract / FREE Full Text ↵ Chen , C. , Qi , H. , Shen , Y. , Pickrell , J. & Przeworski , M. Contrasting determinants of mutation rates in germline and soma . Genetics 207 , 255 – 267 ( 2017 ). OpenUrl Abstract / FREE Full Text ↵ Poulos , R. C. et al. Functional mutations form at ctcf-cohesin binding sites in melanoma due to uneven nucleotide excision repair across the motif . Cell reports 17 , 2865 – 2872 ( 2016 ). OpenUrl CrossRef PubMed ↵ Sabarinathan , R. , Mularoni , L. , Deu-Pons , J. , Gonzalez-Perez , A. & López-Bigas , N. Nucleotide excision repair is impaired by binding of transcription factors to dna . Nature 532 , 264 ( 2016 ). OpenUrl CrossRef PubMed ↵ Perera , D. et al. Differential dna repair underlies mutation hotspots at active promoters in cancer genomes . Nature 532 , 259 ( 2016 ). OpenUrl CrossRef PubMed ↵ Mao , P. et al. Ets transcription factors induce a unique uv damage signature that drives recurrent mutagenesis in melanoma . Nature Communications 9 ( 2018 ). ↵ Elliott , K. et al. Elevated pyrimidine dimer formation at distinct genomic bases underlies promoter mutation hotspots in uv-exposed cancers . PLoS genetics 14 , e1007849 ( 2018 ). OpenUrl CrossRef ↵ Kaiser , V. B. et al. Mutational bias in spermatogonia impacts the anatomy of regulatory sites in the human genome . Genome Research 31 , 1994 – 2007 ( 2021 ). OpenUrl Abstract / FREE Full Text ↵ Seplyarskiy , V. et al. A mutation rate model at the basepair resolution identifies the mutagenic effect of polymerase iii transcription . Nature Genetics 1 – 8 ( 2023 ). ↵ Karczewski , K. J. et al. The mutational constraint spectrum quantified from variation in 141,456 humans . Nature 581 , 434 – 443 ( 2020 ). OpenUrl CrossRef PubMed ↵ Campbell , P. J. et al. Pan-cancer analysis of whole genomes . Nature 2020 578:7793 578 , 82 – 93 ( 2020 ). Publisher: Nature Publishing Group . OpenUrl CrossRef PubMed ↵ Carlson , J. et al. Extremely rare variants reveal patterns of germline mutation rate heterogeneity in humans . Nature communications 9 , 3753 ( 2018 ). OpenUrl CrossRef PubMed ↵ Seplyarskiy , V. B. et al. Population sequencing data reveal a compendium of mutational processes in the human germ line . Science 373 , 1030 – 1035 ( 2021 ). OpenUrl Abstract / FREE Full Text ↵ Pleasance , E. D. et al. A comprehensive catalogue of somatic mutations from a human cancer genome . Nature 463 , 191 – 196 ( 2010 ). OpenUrl CrossRef PubMed Web of Science ↵ Dietlein , F. et al. Identification of cancer driver genes based on nucleotide context . Nature Genetics 52 , 208 – 218 ( 2020 ). OpenUrl CrossRef PubMed ↵ Cooper , D. N. & Krawczak , M. Cytosine methylation and the fate of cpg dinucleotides in vertebrate genomes . Human genetics 83 , 181 – 188 ( 1989 ). OpenUrl CrossRef PubMed Web of Science ↵ Bickmore , W. A. & Van Steensel , B. Genome architecture: domain organization of interphase chromosomes . Cell 152 , 1270 – 1284 ( 2013 ). OpenUrl CrossRef PubMed Web of Science ↵ Odegard , V. H. & Schatz , D. G. Targeting of somatic hypermutation . Nature Reviews Immunology 6 , 573 – 583 ( 2006 ). OpenUrl CrossRef PubMed Web of Science ↵ Sigova , A. A. et al. Divergent transcription of long noncoding rna/mrna gene pairs in embryonic stem cells . Proceedings of the National Academy of Sciences 110 , 2876 – 2881 ( 2013 ). OpenUrl Abstract / FREE Full Text ↵ Hodgkinson , A. & Eyre-Walker , A. Variation in the mutation rate across mammalian genomes . Nature reviews genetics 12 , 756 – 766 ( 2011 ). OpenUrl CrossRef PubMed ↵ Petryk , N. et al. Replication landscape of the human genome . Nature communications 7 , 10208 ( 2016 ). OpenUrl CrossRef PubMed ↵ Jonsson , H. et al. Differences between germline genomes of monozygotic twins . Nature Genetics 53 , 27 – 34 ( 2021 ). OpenUrl CrossRef PubMed ↵ Jónsson , H. et al. Parental influence on human germline de novo mutations in 1,548 trios from iceland . Nature 549 , 519 – 522 ( 2017 ). OpenUrl CrossRef PubMed ↵ Sasani , T. A. et al. Large, three-generation human families reveal post-zygotic mosaicism and variability in germline mutation accumulation . Elife 8 , e46922 ( 2019 ). OpenUrl CrossRef PubMed ↵ Ju , Y. S. et al. Somatic mutations reveal asymmetric cellular dynamics in the early human embryo . Nature 543 , 714 – 718 ( 2017 ). OpenUrl CrossRef PubMed Web of Science ↵ Rodin , R. E. et al. The landscape of somatic mutation in cerebral cortex of autistic and neurotypical individuals revealed by ultra-deep whole-genome sequencing . Nature neuroscience 24 , 176 – 185 ( 2021 ). OpenUrl CrossRef PubMed ↵ Goldmann , J. M. et al. Germline de novo mutation clusters arise during oocyte aging in genomic regions with high double-strand-break incidence . Nature genetics 50 , 487 – 492 ( 2018 ). OpenUrl CrossRef PubMed ↵ Goldmann , J. M. et al. Parent-of-origin-specific signatures of de novo mutations . Nature genetics 48 , 935 – 939 ( 2016 ). OpenUrl CrossRef PubMed ↵ C Yuen , R. K. et al. Whole genome sequencing resource identifies 18 new candidate genes for autism spectrum disorder . Nature neuroscience 20 , 602 – 611 ( 2017 ). OpenUrl CrossRef PubMed ↵ An , J.-Y. et al. Genome-wide de novo risk score implicates promoter variation in autism spectrum disorder . Science 362 , eaat6576 ( 2018 ). OpenUrl Abstract / FREE Full Text ↵ Halldorsson , B. V. et al. Characterizing mutagenic effects of recombination through a sequence-level genetic map . Science 363 , eaau1043 ( 2019 ). OpenUrl Abstract / FREE Full Text ↵ Maury , E. A. et al. Somatic mosaicism in schizophrenia brains reveals prenatal mutational processes . Science 386 , 217 – 224 ( 2024 ). OpenUrl CrossRef PubMed ↵ Lensing , S. V. et al. Dsbcapture: in situ capture and sequencing of dna breaks . Nature methods 13 , 855 – 857 ( 2016 ). OpenUrl CrossRef PubMed ↵ Wang , Z. et al. Combinatorial patterns of histone acetylations and methylations in the human genome . Nature genetics 40 , 897 – 903 ( 2008 ). OpenUrl CrossRef PubMed Web of Science ↵ Ginno , P. A. , Lott , P. L. , Christensen , H. C. , Korf , I. & Chédin , F. R-loop formation is a distinctive characteristic of unmethylated human cpg island promoters . Molecular cell 45 , 814 – 825 ( 2012 ). OpenUrl CrossRef PubMed Web of Science ↵ Core , L. J. et al. Analysis of nascent rna identifies a unified architecture of initiation regions at mammalian promoters and enhancers . Nature Genetics 46 , 1311 – 1320 ( 2014 ). OpenUrl CrossRef PubMed ↵ Hothorn , T. , Bretz , F. & Westfall , P. Simultaneous inference in general parametric models . Biometrical Journal 50 ( 2008 ). ↵ Bhatia , V. et al. Brca2 prevents r-loop accumulation and associates with trex-2 mrna export factor pcid2 . Nature 511 , 362 – 365 ( 2014 ). OpenUrl CrossRef PubMed Web of Science ↵ Shivji , M. K. , Renaudin , X. , Williams , Ç. H. & Venkitaraman , A. R. Brca2 regulates transcription elongation by rna polymerase ii to prevent r-loop accumulation . Cell reports 22 , 1031 – 1039 ( 2018 ). OpenUrl CrossRef PubMed ↵ Hinch , R. , Donnelly , P. & Hinch , A. G. Meiotic dna breaks drive multifaceted mutagenesis in the human germ line . Science 382 , eadh2531 ( 2023 ). OpenUrl CrossRef PubMed ↵ Alexandrov , L. B. et al. The repertoire of mutational signatures in human cancer . Nature 578 , 94 – 101 ( 2020 ). Number: 7793 Publisher: Nature Publishing Group . OpenUrl CrossRef PubMed ↵ Nik-Zainal , S. et al. Landscape of somatic mutations in 560 breast cancer whole genome sequences . Nature 534 , 47 – 54 ( 2016 ). OpenUrl CrossRef PubMed ↵ Wyatt , D. W. et al. Essential roles for polymerase ✓-mediated end joining in the repair of chromosome breaks . Molecular cell 63 , 662 – 673 ( 2016 ). OpenUrl CrossRef PubMed ↵ Serrano Colome , C. , Canal Anton , O. , Seplyarskiy , V. & Weghorn , D. Mutational signature decomposition with deep neural networks reveals origins of clock-like processes and hypoxia dependencies . bioRxiv ( 2023 ). Publisher: Cold Spring Harbor Laboratory eprint : https://www.biorxiv.org/content/early/2023/12/08/2023.12.06.570467.full.pdf . ↵ Cosmic mutational signatures (v3.4 - october 2023 ). signatures/sbs/. Accessed: 2023-10-24 . https://cancer.sanger.ac.uk/ ↵ Polak , P. & Arndt , P. F. Transcription induces strand-specific mutations at the 5’ end of human genes . Genome Research 18 , 1216 – 1223 ( 2008 ). OpenUrl Abstract / FREE Full Text ↵ Seplyarskiy , V. B. & Sunyaev , S. The origin of human mutation in light of genomic data . Nature Reviews Genetics 1 – 15 ( 2021 ). ↵ Mas-Ponte , D. & Supek , F. Mutation rate heterogeneity at the sub-gene scale due to local dna hypomethylation . Nucleic Acids Research gkae252 ( 2024 ). ↵ Seplyarskiy , V. B. et al. Error-prone bypass of dna lesions during lagging-strand replication is a common source of germline and cancer mutations . Nature genetics 51 , 36 – 41 ( 2019 ). OpenUrl CrossRef PubMed ↵ Lek , M. et al. Analysis of protein-coding genetic variation in 60,706 humans . Nature 536 , 285 ( 2016 ). OpenUrl CrossRef PubMed Web of Science ↵ Charlesworth , D. , Charlesworth , B. & Morgan , M. The pattern of neutral molecular variation under the background selection model . Genetics 141 , 1619 – 1632 ( 1995 ). OpenUrl Abstract / FREE Full Text ↵ Rodriguez-Galindo , M. , Casillas , S. , Weghorn , D. & Barbadilla , A. Germline de novo mutation rates on exons versus introns in humans . Nature Communications 11 , 3304 ( 2020 ). Number: 1 Publisher: Nature Publishing Group . OpenUrl CrossRef PubMed ↵ Duret , L. & Galtier , N. Biased gene conversion and the evolution of mammalian genomic landscapes . Annual review of genomics and human genetics 10 , 285 – 311 ( 2009 ). OpenUrl CrossRef PubMed Web of Science ↵ Supek , F. & Lehner , B. Differential dna mismatch repair underlies mutation rate variation across the human genome . Nature 521 , 81 – 84 ( 2015 ). OpenUrl CrossRef PubMed ↵ Dal , G. M. et al. Early postzygotic mutations contribute to de novo variation in a healthy monozygotic twin pair . Journal of medical genetics 51 , 455 – 459 ( 2014 ). OpenUrl Abstract / FREE Full Text ↵ Acuna-Hidalgo , R. et al. Post-zygotic point mutations are an underrecognized source of de novo genomic variation . The American Journal of Human Genetics 97 , 67 – 74 ( 2015 ). OpenUrl CrossRef PubMed ↵ Rahbari , R. et al. Timing, rates and spectra of human germline mutation . Nature genetics 48 , 126 – 133 ( 2016 ). OpenUrl CrossRef PubMed ↵ Dukler , N. , Mughal , M. R. , Ramani , R. , Huang , Y.-F. & Siepel , A. Extreme purifying selection against point mutations in the human genome . Nature communications 13 , 4312 ( 2022 ). OpenUrl CrossRef PubMed ↵ Chen , S. et al. A genomic mutational constraint map using variation in 76,156 human genomes . Nature 625 , 92 – 100 ( 2024 ). OpenUrl CrossRef PubMed ↵ Halldorsson , B. V. et al. The sequences of 150,119 genomes in the uk biobank . Nature 607 , 732 – 740 ( 2022 ). OpenUrl CrossRef PubMed ↵ Murat , P. et al. Dna replication initiation shapes the mutational landscape and expression of the human genome . Science Advances 8 , eadd3686 ( 2022 ). OpenUrl CrossRef PubMed ↵ Li , R. et al. A high-resolution map of non-crossover events reveals impacts of genetic diversity on mammalian meiotic recombination . Nature communications 10 , 3900 ( 2019 ). OpenUrl CrossRef PubMed ↵ Severin , J. et al. ZENBU-Reports: a graphical web-portal builder for interactive visualization and dissemination of genome-scale data . NAR Genomics and Bioinformatics 5 ( 2023 ). ↵ Severin , J. et al. Interactive visualization and analysis of large-scale sequencing datasets using ZENBU . Nature Biotechnology 32 , 217 – 219 ( 2014 ). OpenUrl CrossRef PubMed ↵ Noguchi , S. et al. FANTOM5 CAGE profiles of human and mouse samples . Scientific Data 4 , 170112 ( 2017 ). Number: 1 Publisher: Nature Publishing Group . OpenUrl CrossRef PubMed ↵ Derrien , T. et al. Fast computation and applications of genome mappability . PLoS One 7 , e30377 ( 2012 ). OpenUrl CrossRef PubMed ↵ Hinrichs , A. S. et al. The UCSC genome browser database: update 2006 . Nucleic Acids Res . 34 , D590 – 8 ( 2006 ). OpenUrl CrossRef PubMed Web of Science ↵ Frankish , A. et al. GENCODE reference annotation for the human and mouse genomes . Nucleic Acids Res . 47 , D766 – D773 ( 2019 ). OpenUrl CrossRef PubMed ↵ Cooper , G. M. et al. Distribution and intensity of constraint in mammalian genomic sequence . Genome Research 15 , 901 – 913 ( 2005 ). OpenUrl Abstract / FREE Full Text ↵ Sudlow , C. et al. UK biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age . PLoS Med . 12 , e1001779 ( 2015 ). OpenUrl CrossRef PubMed ↵ Francioli , L. C. et al. Genome-wide patterns and properties of de novo mutations in humans . Nature Genetics 47 , 822 ( 2015 ). OpenUrl CrossRef PubMed ↵ Richter , F. et al. Genomic analyses implicate noncoding de novo variants in congenital heart disease . Nature genetics 52 , 769 – 777 ( 2020 ). OpenUrl CrossRef PubMed ↵ Chihara , L. M. & Hesterberg , T. C. Introduction to confidence intervals . In Mathematical Statistics with Resampling and R, chap . 5 , 103 – 148 ( John Wiley & Sons, Ltd , 2018 ). URL https://onlinelibrary.wiley.com/doi/abs/10.1002/9781119505969.ch5 . https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781119505969.ch5 . OpenUrl ↵ R Core Team . R: A Language and Environment for Statistical Computing . R Foundation for Statistical Computing, Vienna, Austria ( 2022 ). URL https://www.R-project.org/ . ↵ Lawson , C. L. & Hanson , R. J. Linear least squares with linear inequality constraints . In Solving Least Squares Problems , chap. 23, 158 – 173 ( SIAM , 1995 ), 2 edn. URL https://epubs.siam.org/doi/abs/10.1137/1.9781611971217.ch23 . https://epubs.siam.org/doi/pdf/10.1137/1.9781611971217.ch23 . ↵ Mullen , K. M. & van Stokkum , I. H. M. nnls: The Lawson-Hanson Algorithm for Non-Negative Least Squares (NNLS) ( 2023 ). URL https://CRAN.R-project.org/package=nnls . R package version 1.5. ↵ Kent , W. J. , Zweig , A. S. , Barber , G. , Hinrichs , A. S. & Karolchik , D. BigWig and BigBed: enabling browsing of large distributed datasets . Bioinformatics 26 , 2204 – 2207 ( 2010 ). OpenUrl CrossRef PubMed Web of Science ↵ Neph , S. et al. BEDOPS: high-performance genomic feature operations . Bioinformatics 28 , 1919 – 1920 ( 2012 ). OpenUrl CrossRef PubMed Web of Science ↵ Quinlan , A. R. & Hall , I. M. BEDTools: a flexible suite of utilities for comparing genomic features . Bioinformatics 26 , 841 – 842 ( 2010 ). OpenUrl CrossRef PubMed Web of Science ↵ Nassar , L. R. et al. The UCSC genome browser database: 2023 update . Nucleic Acids Res . 51 , D1188 – D1195 ( 2023 ). OpenUrl CrossRef PubMed ↵ Barrett , T. et al. NCBI GEO: archive for functional genomics data sets–update . Nucleic Acids Res . 41 , D991 – 5 ( 2013 ). OpenUrl CrossRef PubMed Web of Science ↵ Consortium , E. P. et al. An integrated encyclopedia of dna elements in the human genome . Nature 489 , 57 ( 2012 ). OpenUrl CrossRef PubMed Web of Science ↵ Komaki , S. et al. iMETHYL: an integrative database of human DNA methylation, gene expression, and genomic variation . Hum. Genome Var . 5 , 18008 ( 2018 ). OpenUrl CrossRef PubMed ↵ Gaffney , D. J. et al. Controls of nucleosome positioning in the human genome . PLoS Genet . 8 , e1003036 ( 2012 ). OpenUrl CrossRef PubMed ↵ Shtumpf , M. , Piroeva , K. V. , Agrawal , S. P. , Jacob , D. R. & Teif , V. B. NucPosDB: a database of nucleosome positioning in vivo and nucleosomics of cell-free DNA . Chromosoma 131 , 19 – 28 ( 2022 ). OpenUrl CrossRef PubMed ↵ Kong , A. et al. A high-resolution recombination map of the human genome . Nat. Genet . 31 , 241 – 247 ( 2002 ). OpenUrl CrossRef PubMed Web of Science ↵ Cortes-Guzman , M. A. reptDir: annotate replication timing direction ( 2024 ). URL https://github.com/MikeACG/reptDir/blob/main/reptDir_0.0.0.9000.pdf . Package github: https://github.com/MikeACG/reptDir . ↵ Morganella , S. et al. The topography of mutational processes in breast cancer genomes . Nature Communications 7 , 11383 ( 2016 ). OpenUrl CrossRef PubMed ↵ Lavinia , G. All cpg sites for hg19. v1 . Accessed: 2021-05-23 . doi: 10.6084/m9.figshare.1415416 . OpenUrl CrossRef ↵ Gardiner-Garden , M. & Frommer , M. CpG islands in vertebrate genomes . J. Mol. Biol . 196 , 261 – 282 ( 1987 ). OpenUrl CrossRef PubMed Web of Science ↵ Qian , S. H. et al. EndoQuad: a comprehensive genome-wide experimentally validated endogenous g-quadruplex database . Nucleic Acids Res . 52 , D72 – D80 ( 2024 ). OpenUrl CrossRef PubMed ↵ Kinsella , R. J. et al. Ensembl BioMarts: a hub for data retrieval across taxonomic space . Database (Oxford) 2011 , bar030 ( 2011 ). OpenUrl CrossRef PubMed ↵ Rauluseviciute , I. et al. JASPAR 2024: 20th anniversary of the open-access database of transcription factor binding profiles . Nucleic Acids Res . 52 , D174 – D182 ( 2024 ). OpenUrl CrossRef PubMed ↵ Xia , X. Position weight matrix, gibbs sampler, and the associated significance tests in motif characterization and prediction . Scientifica (Cairo) 2012 , 917540 ( 2012 ). OpenUrl PubMed ↵ Yang , C. , Bolotin , E. , Jiang , T. , Sladek , F. M. & Martinez , E. Prevalence of the initiator over the TATA box in human and yeast genes and identification of DNA motifs enriched in human TATA-less core promoters . Gene 389 , 52 – 65 ( 2007 ). OpenUrl CrossRef PubMed Web of Science ↵ Kong , A. et al. Rate of de novo mutations and the importance of father’s age to disease risk . Nature 488 , 471 ( 2012 ). OpenUrl CrossRef PubMed Web of Science ↵ GTEx Consortium . The Genotype-Tissue expression (GTEx) project . Nat. Genet . 45 , 580 – 585 ( 2013 ). OpenUrl CrossRef PubMed ↵ PCAWG Transcriptome Core Group et al . Genomic basis for RNA alterations in cancer . Nature 578 , 129 – 136 ( 2020 ). OpenUrl CrossRef PubMed ↵ Wang , T. et al. Identification and characterization of essential genes in the human genome . Science 350 , 1096 – 1101 ( 2015 ). OpenUrl Abstract / FREE Full Text ↵ Weghorn , D. & Sunyaev , S. Bayesian inference of negative and positive selection in human cancers . Nature Genetics 49 , 1785 ( 2017 ). OpenUrl CrossRef PubMed ↵ Futreal , P. A. et al. A census of human cancer genes . Nat. Rev. Cancer 4 , 177 – 183 ( 2004 ). OpenUrl CrossRef PubMed Web of Science ↵ Core , L. J. , Waterfall , J. J. & Lis , J. T. Nascent rna sequencing reveals widespread pausing and divergent initiation at human promoters . Science 322 , 1845 – 1848 ( 2008 ). OpenUrl Abstract / FREE Full Text ↵ Coxe , S. , West , S. G. & Aiken , L. S. The analysis of count data: a gentle introduction to poisson regression and its alternatives . J Pers Assess 91 , 121 – 136 ( 2009 ). OpenUrl CrossRef PubMed Web of Science ↵ Venables , W. N. & Ripley , B. D. Generalized linear models . In Modern Applied Statistics with S , 183 – 210 ( Springer , New York, NY , 2002 ), 4 edn. ↵ Fernandez , G. A. & Vatcheva , K. P. A comparison of statistical methods for modeling count data with an application to hospital length of stay . BMC Med. Res. Methodol . 22 , 211 ( 2022 ). OpenUrl CrossRef PubMed ↵ Little , T. D. Darlington , R. B. & Hayes , A. F. Linear interaction . In Little , T. D. (ed.) Regression Analysis and Linear Models: Concepts, Applications, and Implementation , 377 – 408 ( The Guildford Press , New York, NY , 2016 ). ↵ Brehm , L. & Alday , P. M. Contrast coding choices in a decade of mixed models . Journal of Memory and Language 125 , 104334 ( 2022 ). URL https://www.sciencedirect.com/science/article/pii/S0749596X22000213 . OpenUrl CrossRef ↵ Schielzeth , H. Simple means to improve the interpretability of regression coefficients . Methods in Ecology and Evolution 1 , 103 – 113 ( 2010 ). URL https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/j.2041-210X.2010.00012.x . https://besjournals.onlinelibrary.wiley.com/doi/pdf/10.1111/j.2041-210X.2010.00012.x . OpenUrl CrossRef ↵ Vatcheva , K. P. , Lee , M. , McCormick , J. B. & Rahbar , M. H. Multicollinearity in regression analyses conducted in epidemiologic studies . Epidemiology (Sunnyvale) 6 ( 2016 ). ↵ Fox , J. & Weisberg , S. Regression diagnostics . In An R Companion to Applied Regression , 429 – 434 ( Sage , Thousand Oaks CA , 2019 ), 3 edn. ↵ Fox , J. & Monette , G. Generalized collinearity diagnostics . Journal of the American Statistical Association 87 , 178 – 183 ( 1992 ). URL https://www.tandfonline.com/doi/abs/10.1080/01621459.1992.10475190 . https://www.tandfonline.com/doi/pdf/10.1080/01621459.1992.10475190 . OpenUrl CrossRef Web of Science ↵ Zuur , A. F. , Ieno , E. N. & Elphick , C. S. A protocol for data exploration to avoid common statistical problems . Methods in Ecology and Evolution 1 , 3 – 14 ( 2010 ). URL https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/j.2041-210X.2009.00001.x . https://besjournals.onlinelibrary.wiley.com/doi/pdf/10.1111/j.2041-210X.2009.00001.x . OpenUrl CrossRef ↵ Hartig , F. Dharma: residual diagnostics for hierarchical (multi-level/mixed) regression models . https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMa.html ( 2022 ). Accessed: 2024-09-06 . ↵ Lee , W. , Kim , J. & Lee , D. Revisiting the analysis pipeline for overdispersed poisson and binomial data . J. Appl. Stat . 50 , 1455 – 1476 ( 2023 ). OpenUrl CrossRef PubMed ↵ Hartig , F. DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models ( 2022 ). URL https://CRAN.R-project.org/package=DHARMa . R package version 0.4.6. View the discussion thread. Back to top Previous Next Posted February 05, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Transcription start sites experience a high influx of heritable variants fuelled by early development Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Transcription start sites experience a high influx of heritable variants fuelled by early development Miguel Cortés Guzmán , David Castellano , Clàudia Serrano Colomé , Vladimir Seplyarskiy , Donate Weghorn bioRxiv 2025.02.04.635982; doi: https://doi.org/10.1101/2025.02.04.635982 Share This Article: Copy Citation Tools Transcription start sites experience a high influx of heritable variants fuelled by early development Miguel Cortés Guzmán , David Castellano , Clàudia Serrano Colomé , Vladimir Seplyarskiy , Donate Weghorn bioRxiv 2025.02.04.635982; doi: https://doi.org/10.1101/2025.02.04.635982 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genomics Subject Areas All Articles Animal Behavior and Cognition (7637) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41968) Biophysics (21460) Cancer Biology (18603) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13385) Ecology (19910) Epidemiology (2067) Evolutionary Biology (24327) Genetics (15614) Genomics (22513) Immunology (17741) Microbiology (40423) Molecular Biology (17193) Neuroscience (88646) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4825) Physiology (7647) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)