Characterizing cytosine methylation of polymorphic transposable element insertions using the human pangenome resources

Both PacBio HiFi and Nanopore ONT reads provide accurate DNA methylation estimation

The HPRC published the first human pangenome reference graph representing 47 genetically diverse individuals (Liao et al. 2023). Out of 47, the DNA methylation status of 32 were determined using PacBio HiFi reads (Methods). In addition to the rich genetic information and various annotations, the HPRC resources also enabled direct methylation analysis using long-read sequencing data. Taking advantage of these resources, we performed additional CpG methylation calling using ONT reads (R9.4.1) from five individuals who were part of the first draft of the Human Pangenome Reference (HG00741, HG00621, HG01952, HG01978, and HG03516) and the WGBS of these same five individuals (Methods). After calling CpG methylation, we aligned all reads to the human reference genome GRCh38 (Fig. 1A).

Figure 1.

Comparing CpG methylation calling accuracy between PacBio, ONT, and WGBS. (A) Study design of methylation calling benchmarking. We performed CpG methylation calling for five HPRC samples (HG00621, HG00741, HG01952, HG01978, and HG03516), using PacBio, ONT, and WGBS. We mapped all reads to the human reference genome GRCh38 to compare their differences. (B) Whole-genome methylation percentage correlation between WGBS, HiFi, and ONT methylation of HG00741. The whole-genome methylation percentage distribution is represented as a density plot. The correlation heatmap and metrics (mean average error, 20% concordance, and Pearson's correlation) are plotted on the bottom left and top right, respectively. (C) Methylation percentage distribution of the GM12878 15 ChromHMM states across all five samples. The 15 ChromHMM states are listed on the x-axis, and methylation calling results from different methods are colored differently. (D) WashU Epigenome Browser screenshot of the two WGBS replicates, HiFi, and ONT methylation of HG00741 on the 10 kb window around the TP53 gene promoter. All methylC tracks y-axis percentage was from 0% to 100%. (E) PCA plot of the top two principal components of the whole-genome methylation from the five HPRC samples. For each sample, we performed two replicates of WGBS in addition to HiFi and ONT methylation calling. We used different shapes to represent different methods (circle, triangle, and cross for ONT, HiFi, and WGBS, respectively) and different colors to represent different samples (red, yellow, green, blue, and purple for HG00621, HG00741, HG01952, HG01978, and HG03516, respectively). (F) Number of falsely discovered differentially methylated regions (DMRs) using different threshold cutoffs. Four different thresholds (20% average methylation difference with at least five CpG sites, 20% average methylation difference with at least 10 CpG sites, 30% average methylation difference with at least five CpG sites, and 30% average methylation difference with at least 10 CpG sites) are represented on the x-axis. The number of DMRs are represented on the y-axis. The number of DMRs between different methods under different conditions are plotted using different colors.

Open in new tabDownload PowerPointLink to figure

We first compared DNA methylation calls between HiFi, ONT, and WGBS for all CpG sites in the reference genome GRCh38. For HiFi reads, we adopted the deep learning–based mathematical model (PB-model) developed and suggested by PacBio to estimate methylation percentage of individual CpG sites (https://github.com/PacificBiosciences/pb-CpG-tools). We calculated mean average error (MAE), 20% concordance, and Pearson's correlation between WGBS, HiFi, and ONT for each sample in pairwise comparisons (for HG00741 comparison, see Fig. 1B; for the other four, see Supplemental Fig. S1A–D). We found both HiFi and ONT methylation calls highly correlated with their matching WGBS results (MAE: 0.1062 and 0.0975; 20% concordance: 0.846 and 0.886; Pearson's correlation: 0.901 and 0.925 for HiFi and ONT, respectively, for HG00741) in all five samples (Fig. 1B; Supplemental Fig. S1A–D). We examined the distribution of methylation level across different genomic contexts (defined by the GM12878 ChromHMM states) and found strong concordance among WGBS, HiFi, and ONT across all five samples (Fig. 1C; Ernst and Kellis 2012). We then used the WashU Epigenome Browser to illustrate the comparisons between WGBS, HiFi, and ONT (Fig. 1D).

To gain a global view of the difference and similarity across methods, we performed principal component analysis (PCA) and plotted the first two principal components in a 2D plot for visualization (Fig. 1E). In all five samples, the two matching WGBS replicates consistently clustered together; HiFi and ONT methylation results in the same samples also clustered together. The PCA distance between WGBS and the long-read methylation results seemed to be larger than the distances between HiFi and ONT.

To better understand the scale of differences between the different methods, we defined both differentially methylated loci (DMLs) and differentially methylated regions (DMRs) between HiFi, ONT, and WGBS, as well as between two WGBS replicates of the same sample (Wu et al. 2015). Because the DMRs here were defined using the same samples, they represented the false-positive discoveries and allowed us to estimate an empirical false-discovery rate. Using different thresholds, we found that the numbers of falsely discovered DMRs between HiFi and ONT were comparable with the numbers derived from the two WGBS replicates (Fig. 1F). Thus, the methylation calling difference between HiFi and ONT was quite small, generally comparable with the differences we observed between the two WGBS replicates.

To further understand the methylation differences between long-read and short-read methods, we examined DMLs and DMRs called between the two methods. We found that the DMLs between long and short reads exhibiting a >50% methylation difference were predominantly contributed by CpG sites with abnormally high short-read coverage (>50×) (Supplemental Fig. S1E), suggesting short-read mapping in these regions was more challenging than in others. In contrast, long-read sequencing coverage or GC content did not correlate with falsely discovered DNA methylation differences (Supplemental Fig. S1F,G). We also found that low complexity and simple repeats contributed to a higher fraction of DMRs with a methylation difference >40% (Supplemental Fig. S1H). Because these regions also suffer more from imprecise mapping of short reads, we reasoned that the small discrepancy between long-read- and short-read-based methods in assessing DNA methylation was systematic and likely a result of reference mapping bias (Pollard et al. 2018). A systematic benchmarking analysis with a focus on alignment bias will be required in the future to fully understand this phenomenon.

The methylation patterns of TE insertions

Because the long-read sequencing methylation calling was conducted on raw reads without aligning to the reference genome, it enabled methylation calling of nonreference regions previously inaccessible with Illumina sequencing–based methods (at least not trivially). Thus, we investigated the methylation status of polymorphic TE insertions and their correlation with the methylation levels of their flanking regions. Here we took advantage of the recent release of the draft Human Pangenome Reference made available by the HPRC, which included HiFi reads with methylation calling from 32 HPRC Year1 samples as well as their genetic variations (Methods) (Liao et al. 2023). We further annotated full-length insertions combining the TE insertion annotation pipeline MELT-LRA, PALMER2, and xTEA (Methods) (Gardner et al. 2017; Chu et al. 2021; Zhou et al. 2024). We classified TE insertions into five categories: Alu, L1, SVA, LTR5, and HERVK. In total, we identified 6752 Alu, 410 L1, 334 SVA, 16 solo LTR5, and nine HERVK insertions from the 32 HPRC samples (Supplemental Table S1). For convenience, we termed these TE insertions nonreference TEs.

Using custom pipelines (Methods), we extracted nonreference TEs and the HiFi reads derived from these regions (Fig. 2A). We first assessed if we had sufficient long-read coverage on nonreference TEs. We calculated the relative CpG site coverage of each TE insertion compared with the genome-wide average coverage based on the HiFi data and observed the two expected peaks for each TE class, representing coverages of heterozygous and homozygous insertions (Supplemental Fig. S2A). The average genome-wide coverage was 45.9×; thus, we believe that the coverage on TE insertions was sufficient for the purpose of our analysis. We also extracted HiFi methylation of reference TEs for comparison (Fig. 2A).

Figure 2.

The methylation profile of different transposable elements (TEs) in the human genome. (A) Schematic representation of our pipeline to extract reference and nonreference insertions along with their methylation annotation. We developed a script (bam.extractInsertion.py) to extract truncated DNA sequences from reads while retaining positional methylation tags (subseq modBAM file). We then aligned this modBAM to the corresponding TE consensus and calculated the methylation profile of the TE insertion. We also extracted sequences from the reference genome and the methylation percentage from the bedMethyl file to create another modBAM file of reference TEs and their methylation (bam.creation.py). The methylation would be anchored to the TE consensus after alignment. (B) Methylation profile of L1 insertions. The x-axis represents positions on the L1Hs consensus, and the y-axis represents the methylation percentage with 95% standard error of the mean illustrated as shaded areas. Reference L1pre-Ta, reference L1Ta, and nonreference L1 insertions are labeled with different colors. Only full-length elements (truncation <50 bp on both ends) are included for profiling. (C) Methylation profile of Alu insertions with 95% standard error of mean. The x-axis represents positions on the AluY consensus, and the y-axis represents the methylation percentage. AluY, AluYa5, AluYa8, AluYb8, and AluYb9 on the reference and nonreference Alu insertions are labeled with different colors. (D) Methylation profile of SVA insertions with 95% standard error of mean. The x-axis represents positions on the SVA_F consensus, and the y-axis represents the methylation percentage. Reference SVA and nonreference SVA insertions are labeled with different colors. The VNTR of SVA consensus are labeled on the graph with a gray box. (E) The first 500 bp methylation of L1 insertions whose first 20–50 bp was truncated. The start position is labeled on the left, and the insertion coordinates on GRCh38 are labeled on the right. The “refTE” indicates it was an L1 insertion on the GRCh38 reference. Each row represents one of the 32 samples. The methylation percentage of each CpG site is illustrated by the blue–red (low to high) color scale.

Open in new tabDownload PowerPointLink to figure

For the five HPRC individuals with WGBS, HiFi, and ONT methylation results, we phased both HiFi and ONT reads to maternal and paternal assemblies and mapped them onto maternal or paternal haplotypes, respectively (Methods) (Patterson et al. 2015; Holt et al. 2024). Using the reference-assembly alignment, we created comparative epigenome browser datahubs to interactively examine the methylation of nonreference insertions (Methods).

In addition to characterizing nonreference TE insertions, to provide a contrast, we also analyzed the methylation of related TE subfamilies within GRCh38 reference genome (termed reference TE insertions). Polymorphic TE insertions include all nonreference and some reference TE insertions, and we plotted both of them in the Figure 2. We extracted 52,802 AluY, 2746 AluYa5, 36 AluYa8, 2111 AluYb8, and 301 AluYb9 from GRCh38 using RepeatMasker annotation. Similarly, we extracted 124 SVA_E, 331 SVA_F, and 481 LTR5 as reference TEs. We also downloaded 146 active L1s on GRCh38 from L1 Base 2 (Penzkofer et al. 2017). Active L1s were further classified into L1Ta and L1Pre-Ta using the defining variant at position 5932 (A as L1Ta and G as L1pre-Ta) based on L1base 2 (Boissinot et al. 2000). Intact AluY, L1HS, and SVA each have a region that contains enough CpG sites to be considered a CGI; thus, we have abundant data to characterize the methylation percentage fluctuation of each TE. We first compared DNA methylation profiles of full-length nonreference TE insertions with reference TE insertions. Nonreference TE insertions represent the youngest TE insertions in the human population, whereas reference TEs from the same subfamilies represent slightly older but still human-specific insertions given that nonreference TE insertions are generally of low allele frequency and that reference TEs are mostly of high allele frequency or are fixed in the human population (Logsdon et al. 2025). By anchoring methylation to a consensus-based alignment, we significantly reduced the noise in methylation signals compared with previous studies (Ewing et al. 2020; Gerdes et al. 2023). This allowed us to observe sharp methylation fluctuations between adjacent CpG sites within TEs (Fig. 2B–D). Our findings suggested that certain CpG sites, such as position 5875 in L1HS, were consistently lowly methylated. Overall, the DNA methylation profiles of reference and nonreference TE insertions showed high methylation levels and were similar in all four classes (for Alu, L1, and SVA profiles, see Fig. 2B–D; for the LTR5 profile, see Supplemental Fig. S2B). However, a single CpG site in each of Alu, L1, and SVA (at position 213, 5875, and 1166, respectively) exhibited statistically significant differences between nonreference insertions and reference copies. Also profoundly, the internal CG-rich, variable number of tandem repeats (VNTR) region of SVA was methylated lower in the nonreference insertions compared with their counterparts in the reference copies (Fig. 2D). The overall methylation plots for all but Alu insertions (which were too numerous to display individually) were available in Supplemental Figure S2, C through E.

To our knowledge, the relatively low methylation of VNTR within the SVA elements has never been reported before. Ewing et al. (2020) found that the VNTRs are methylated higher than their flanking regions in SVA, and we confirmed their observation within reference SVAs (Supplemental Fig. S2D). However, the methylation of VNTR of nonreference SVAs in our LCL samples was noticeably lower compared with both their flanking regions and their counterparts in the reference genome (Fig. 2D). Aligning full-length SVAs, often >2 kb long, to the 1.3 kb consensus may create some alignment issue, especially in the VNTR region. However, we do not believe this could explain our observation: The long-read methylation calling was done during basecalling process, independent of alignment, and the presence of many lowly methylated CpGs was visible across many nonreference SVAs (Supplemental Fig. S2F). We also found no correlation between length and methylation level of VNTR regions (Supplemental Fig. S2G).

The transcription factor YY1 is an important regulator of L1 activity in mammals (Saha et al. 2024). It was reported that L1 insertions missing the YY1 binding site, located at positions from 12 to 20 in the L1HS, could escape 5′ methylation in a tissue-specific manner (Sanchez-Luque et al. 2019). It was also proposed that YY1 binding affects the methylation of L1 in a cell type–specific manner (Lanciano et al. 2024). To investigate whether the YY1 binding site could affect polymorphic L1 methylation in LCL cells, we examined the 5′-end methylation of 10 L1 insertions whose first 20–100 bp were truncated (Fig. 2E).

The YY1 motif is located at position 12–20 in L1 consensus, and none of the 10 L1s here contained the YY1 binding site. Despite that, we found five loci, missing 20 bp to 28 bp from the 5′ end, were hypermethylated at the L1 promoter. In contrast, three out of the five loci with >32 bp truncated at the 5′ end were hypomethylated. Two loci, missing 39 and 73 bp, respectively, remained hypermethylated. Our data thus suggest that missing YY1 binding site is not sufficient for the demethylation of L1 promoter, at least not in LCL. In addition we found individual DNA methylation differences of the same L1 insertions. The reference L1 insertion with the first 31 bp at the 5′ end truncated at Chr 13: 29,641,706–29,647,706, denoted as Chr13delta31 by Sanchez-Luque et al. (2019), displayed variation among our 32 samples: Its promoter was hypomethylated in most of the samples but hypermethylated in HG02717 (Supplemental Fig. S2H). Because of the redundant nature of TE insertions, we could not apply other functional genomic analysis, including Hi-C, to explain the methylation differences in these regions.

Low methylation correlation between TE insertions and their flanking regions

We next investigated the relationship between DNA methylation of the insertions and their flanking regions. We compared the average methylation levels of nonreference insertions, categorized based on their TE annotation, with the methylation levels of their 300 bp flanking regions on both the 5′ end and 3′ end (Fig. 3A). We defined all insertions not annotated as TEs as non-TE insertions. TE insertions, particularly Alu and SVA insertions, exhibited high methylation levels regardless of the methylation levels of their flanking regions (Fig. 3B). On average, we observed weak correlations between TEs and their flanking regions (linear regression slope < 0.2, R²∼ 0.1). We did not calculate the correlation of HERVK insertions because of the insufficient number of nonreference HERVK insertions. This pattern was in sharp contrast to non-TE insertions, whose methylation levels were highly correlated with their nearby flanking regions (slope = 0.88, R² = 0.61) (Fig. 3B). The pattern for non-TE insertions did not change much when we separated them based on whether they contained tandem repeats, on whether they were in CGIs, or on their CpG density (Supplemental Fig. S3A). Reference TE insertions displayed a similar pattern to that of nonreference TEs (Supplemental Fig. S3B). Among reference TE insertions, LTR5 exhibited the highest linear correlation slope of approximately 0.5 (Supplemental Fig. S3B).

Figure 3.

The methylation correlation between human-specific TE insertions and their flanking regions. (A) Schematic representation of average methylation calculation of all insertions and associated 300 bp flanking regions. (B) Methylation correlation between insertions and 300 bp flanking regions. The average methylation of insertions is represented on the y-axis; the methylation of flanking regions, on the x-axis. Alu, LTR, L1, SVA, and non-TE insertions are represented by red, green, blue, purple, and gray dots, respectively, whereas non-TE insertions are labeled gray, serving as a background. Linear regression lines with standard errors are also plotted and labeled with the same color scheme. (C) Methylation correlation between insertions within CpG islands and 300 bp flanking regions. The average methylation of insertions are represented on the y-axis; the methylation of flanking regions, on the x-axis. Alu and non-TE insertions are represented in the same way as they are in B. (D) The WashU Epigenome browser views of an Alu insertion between hypomethylated and hypermethylated regions. The insertion was absent from HG03516 paternal haplotype but present in all the other nine haplotypes. A genome-align track is used to illustrate the insertion site. The Alu insertion methylation is boxed, showing a gradual change of methylation from left to right. Both HiFi and ONT methylC tracks display a methylation percentage from 0% to 100%. (E) Schematic representation of insertion-flanking methylation regression as a function of insertion-flanking distance. For each distance from 100 bp to 10 kb, we performed linear regression between insertion and flanking average methylation and estimated the linear regression slope and standard errors. (F) Linear regression slope between insertion and flanking methylation as a function of the insertion-flanking distance from 10 kb upstream to 10 kb downstream. Standard errors were represented by error bars. Alu, LTR, L1, SVA, and non-TE insertions are plotted in separated plots with different colors.

Open in new tabDownload PowerPointLink to figure

CGI is an important genomic feature, and its methylation status plays a significant role in gene regulation. When we restricted the above analysis to non-TE insertions within annotated human CGIs (as defined in the Introduction), we observed a stronger correlation between the insertions and their flanking methylation (slope = 0.97, R² = 0.90) (Fig. 3C). This correlation extended even to Alu insertions within these CGIs. In contrast to the general hypermethylation of Alu insertions genome-wide, the methylation of Alu insertions within them correlated highly with their flanking methylation (slope = 1.16, R² = 0.82) (Fig. 3C).

Despite the majority of Alus being hypermethylated, a small percentage were hypomethylated. For instance, at Chr 11: 59,565,869, a polymorphic Alu insertion within a hypomethylated CGI was hypomethylated. We plotted the methylation of this Alu along with its 300 bp flanking regions across all 32 samples (Supplemental Fig. S3C). The insertion site was hypomethylated in all samples, and the Alu insertion itself appeared to have adopted the hypomethylation status in all five samples in which it was present. This region can be further examined in sample HG03516, whose Alu insertion was heterozygous (Supplemental Fig. S3D). The regional hypomethylation, as well as the presence of the Alu insertion, was supported by HiFi, ONT, and WGBS analysis (Supplemental Fig. S3D,E).

Several Alu elements exhibited a methylation slope that connected hypomethylated and hypermethylated regions. We illustrated four such cases in Supplemental Figure S3F. At Chr 1: 162,791,517, Chr 2: 218,333,959, and Chr 9: 71,308,789, the 300 bp upstream of the insertion was hypermethylated, whereas the downstream region was hypomethylated. In contrast, at Chr 4: 163,587,802, the upstream region was hypomethylated, whereas the downstream region was hypermethylated. In each of these four cases, the Alu insertion displayed a gradual transition between hypo- and hypermethylation, with some variation observed among samples. We displayed the antisense Alu insertion at Chr 1: 162,791,517 in our comparative browser by phasing reads to paternal and maternal alleles (Fig. 3D). The ONT and HiFi data showed high agreement with each other on the methylation level of this nonreference insertion. Downstream from the Alu was hypomethylated (left side), whereas upstream of it (right side) was hypermethylated in all samples. In addition to that, there were variations among individuals and between haplotypes at this locus. Although most instances of this insertion were highly methylated throughout the Alu element, the methylation level in the HG03516 maternal allele exhibited a transition from low to high in the middle of the element. In contrast, this insertion was hypomethylated in the HG01952 paternal allele. No other sequence feature seemed to be associated with the methylation pattern difference across individual alleles around Chr 1: 162,791,517.

Our analysis thus far has revealed a strong correlation in DNA methylation levels between non-TE insertions and their 300 bp flanking regions and revealed a much weaker correlation when the insertions were TEs. We next investigated the insertion-flanking methylation correlation beyond 300 bp as a function of the distance between the insertion and the nearby regions. For each nonreference insertion, we calculated DNA methylation levels in its 10 kb flanking regions in both directions at a resolution of 100 bp nonoverlapping windows (Fig. 3E). We grouped the average methylation of all 100 bp windows from different insertions by their distance from the insertion sites and performed linear regression between the insertion and flanking regions. Because the slopes of the linear regression can reflect the methylation correlation between these regions, we displayed the slope of the linear regression as a function of the distance between the flanking windows and the insertion (Fig. 3F). We used non-TE insertions as a control in the analysis because they can be used to represent all factors affecting nearby methylation correlation including linkage disequilibrium decay. This analysis revealed that the methylation correlation slope between non-TE insertions and their flanking windows quickly faded away as the distance between the flanking window and insertion site increased, with an estimated elbow point being ∼2.5 kb away from the insertion. The slope decayed to and stayed at around 0.1 when the distance reached ≥5 kb. TEs once again displayed a distinct pattern. There was no observable correlation between L1, SVA, and HERV insertions and their flanking regions. We did however observe a weak correlation between Alu insertions and their flanking regions, but it also quickly faded away (Fig. 3F). The pattern of Alu was mostly driven by Alu insertions outside of annotated human CGIs (Supplemental Fig. S3G). This observation could be confounded if TE insertions were longer than non-TE insertions. To rule out this possibility, we investigated how the size of non-TE insertion could affect the insertion-flanking methylation correlation. If the methylation correlation is merely a function of the insertion size, we would expect ∼300 bp and ∼1000 bp non-TE insertions to behave similarly to Alu and LTR insertions, for example. We categorized non-TE insertions by their sizes into five bins and then calculated the regression slope in each bin (Supplemental Fig. S3H). Despite a slight decrease of the peak methylation correlation as a function of increased non-TE insertion size, the correlation was much higher than that of TE insertions of comparable length (e.g., the slope distribution of 300–1000 bp non-TE insertions was much higher than that of Alu and LTR insertions). This result suggests that the methylation of non-TE insertions, but not TE insertions, is highly correlated with that of their flanking regions.

Most TE insertions do not change the methylation level of their flanking regions

Some nonreference insertions could be deletion events in the human population (Lee et al. 2020b). However, TE insertions are essentially homoplasy-free (Batzer and Deininger 2002) and represent recent polymorphic insertions in the human population. Therefore, they are excellent candidates to investigate how new insertions might have influenced the methylation of their inserted regions.

Previous studies have suggested that repressive chromatin marks including DNA methylation deposited on newly inserted TEs can spread beyond TE/genome boundary (Turker 2002; Sentmanat and Elgin 2012; Zhang et al. 2012). On the other hand, hypomethylated L1 could drag down the methylation at its 5′ end (Grandi et al. 2015). It was also found that both hypermethylated and hypomethylated L1 insertions can impact nearby methylation in different cell lines (Lanciano et al. 2024).

The HPRC resources provided a unique opportunity to investigate the impact of TE insertions on nearby regions. For each polymorphic nonreference insertion relative to GRCh38, we separated HiFi reads from the 32 HPRC Year1 samples into two groups: alleles with the insertion and alleles without the insertion (Fig. 4A). We then calculated the average methylation level of the inserted regions (insertion), the 300 bp flanking regions of the insertion from both ends (insertion flanking), and the same flanking regions of the allele without insertion (empty flanking) (Fig. 4A).

Figure 4.

Most of the newly inserted TEs did not alter the methylation pattern of their surrounding genomic context. (A) Schematic representation of a polymorphic insertion with aligned long reads. The allele without the insertion site (top) and the allele with the insertion site (bottom) were represented as reference allele and insertion allele. The average methylation level of the insertion, insertion flanking regions, and empty flanking regions were calculated from reads from all 32 samples. We then calculated flanking methylation difference by subtracting the methylation level of empty flanking region from that of insertion flanking region. (B) Methylation level of empty flanking (x-axis), insertion flanking (y-axis), and non-TE inserted regions (blue–red color scale as shown in C). Insertions were separated into four facets by whether they are tandem repeats and whether they are CpG islands. Linear regression R²-values were shown on each panel. (C) Similar to B, but the methylation level of different TE insertions was presented and separated by the TE insertion types. R² were displayed in the top left corner on each panel. (D) Dot plot illustrated the relationship between methylation level of insertion (y-axis) and methylation difference between empty and insertion alleles of the flanking regions (x-axis). The average methylation difference of all non-TE insertions was labeled at 0.25%. Red dotted lines indicated the methylation difference threshold of ±25%, and 99% of non-TE insertions fall within this range. The insertion methylation density was displayed on the right side of the figure. (E) TE insertion flanking methylation differences plotted in the sample way as in D. Different TE insertions were distinguished by different colors. The percentage of all insertions within 25% differences, and the average flanking methylation differences of different TEs were labeled. The insertion methylation density was shown on the right of the figure. The average methylation of each TE class was labeled with the same color scheme. (F) Alu and L1 insertion flanking methylation differences are not higher when the TEs were hypermethylated. We separated all Alu and L1 insertions to above average and below average based on the average TE methylation. The flanking methylation differences are plotted as a violin plot. The flanking methylation differences distribution is similar in all groups, and each of them has a wide range. (G) Methylation level of 300 bp upstream of the empty allele (x-axis), 300 bp upstream of the insertion allele (y-axis), and the first 500 bp inserted methylation (blue–red color scale) of L1 and SVA insertions displayed similar to B and C.

Open in new tabDownload PowerPointLink to figure

To illustrate the relationship between these three regions, we plotted the methylation levels of insertion flanking and empty flanking on the x–y-axis using a dot plot. We colored the dots based on the methylation level of insertions. We separated non-TE insertions by whether they were annotated as tandem repeat or CGIs and classified TE insertions by TE classes. In general, the methylation levels of insertion, insertion flanking, and empty flanking correlated with each other for non-TE insertions (Fig. 4B). However, for TE insertions, the flanking methylation levels of insertion and empty sites showed a strong correlation independent of the methylation of the insertion (Fig. 4C).

To better understand the subtle changes in flanking methylation levels between insertion and empty sites, we defined the flanking methylation difference for each insertion by subtracting the flanking methylation of the empty site from the flanking methylation of the insertion site (Fig. 4A). We reasoned that an allele without a TE insertion represents the ancestral allele (Batzer and Deininger 2002). If the insertion increases flanking methylation, the difference would be positive; if it decreases methylation, the difference would be negative. Conversely, if the insertion has no impact on flanking region methylation, the difference would be zero.

We first plotted the correlation between flanking methylation variations and the methylation of non-TE insertions (Fig. 4D). In this panel, the x-axis represents the differential flanking methylation, calculated as the methylation level of the allele containing the insertion minus the methylation of the empty site (the allele without the insertion). The y-axis represents the methylation level of the insertion itself. The overall methylation differences between orthologous flanking regions with or without insertions were predominantly centered around zero (population average was 0.25%): 99% of the data points (28,344/28,621) fell within the −25% to +25% range. Notably, these variations were independent of the methylation level of the insertion itself. However, there indeed existed multiple instances in which hypermethylated insertions were associated with increased flanking methylation, and hypomethylated insertions were linked to decreased methylation (Fig. 4D). Collectively, these findings strongly indicate that, although there are noted exceptions, non-TE insertions generally do not significantly impact flanking region methylation.

We then investigated the same relationship in the context of TE insertions (Fig. 4E). Compared with non-TE insertions, the range of methylation differences was slightly wider (96%, 5075/5263, within −25% and 25%), indicating significantly increased variations in flanking region methylation associated with TE insertions (F-test to compare two variances: P-value < 2.2 × 10⁻¹⁶). On average, flanking methylation differences centered around zero for LTR and SVA insertions. However, the average flanking methylation differences for Alu and L1 insertions were significantly larger than zero (3.2% and 3.9%, respectively) (Fig. 4E). Furthermore, there were more Alu and L1 insertions with methylation increase >25% than decrease >25% (112 vs. 62, and four vs. one, respectively) (Fig. 4E). When we separated Alu and L1 insertions by above and below average Alu/L1 methylation, we did not find much difference between them (Fig. 4F). The range of methylation differences we observed was also quite large. For example, the mean methylation level of inserted Alus was 82.7%, and Alu insertions with a methylation level greater than this mean value were associated with a flanking region methylation difference ranging from −42.74% to +42.66% (Fig. 4F), with 21 of them associated with an increased flanking region methylation by >30%. Similarly, the mean methylation level of inserted L1s was 69.2%, and L1 insertions with methylation level greater than this mean value were associated with a flanking region methylation difference ranging from −30.1% to +27.8%, but none of them was associated with an increased flanking region methylation by >30% (Fig. 4F). These examples only weakly support the hypothesis that DNA methylation deposited on newly inserted TEs could spread and increase local DNA methylation level. Contrary to this conventional view, our genome-wide assessment suggests that the spreading is limited, and this more drastic influence of TEs on nearby genomic methylation may be an exception rather than the rule.

It has been reported that L1 promoter methylation can influence the nearby upstream methylation (Lanciano et al. 2024). To test this hypothesis using our data, we created a dotplot similar to Figure 4C, but using the 5′-end 500 bp of L1/SVA insertions along with the 300 bp upstream region genomic sequences (Fig. 4G). Notably, despite the existence of a few exceptions, the upstream methylation of empty sites consistently matched those with insertions, indicating that the vast majority of L1/SVA insertions did not affect their flanking upstream methylation.

To further validate and generalize our findings, we aggregated all polymorphic insertions and plotted the average methylation trendline based on their relative position in relationship to the insertion sites from the reference genome (Supplemental Fig. S4A). We found a subtle increase of methylation flanking the polymorphic L1/Alu insertions, with the highest increase right at the insertion site (3.1% for L1 and 3.2% for Alu). This flanking methylation difference between insertion and empty site decreased as the distance to insertion position increased, and the methylation difference diminished at ∼1 kb away from the insertion site (Supplemental Fig. S4A). The exact same pattern was recapitulated by both the ONT methylation data (phased ONT reads) and the WGBS data (WGBS reads of homozygous insertions to the individual's own genome assemblies) from the five HPRC individuals (Methods) (Supplemental Fig. S4B,C). Taken together, these results suggest that polymorphic Alu and L1 insertions are associated with slightly elevated DNA methylation levels at their insertion sites with highly localized and small effect size.

The exceptions: TE insertions associated with flanking region methylation change

Despite most TE insertions not affecting the methylation of their insertion sites, we also found numerous TE insertions associated with flanking methylation changes (Supplemental Table S3). For instance, the hypermethylated promoter of L1 insertion at Chr 3: 85,527,420 (average promoter methylation level 88%) was associated with an increase in methylation of the CpG site upstream of the insertion. We plotted the methylation of each CpG site of the 5′-end 1000 bp of the L1 along with both flanking regions across all 32 samples in Supplemental Figure S5A and summarized them in Figure 5A. In samples without this insertion, the upstream CpG site was mostly hypomethylated (average methylation 24.3%), whereas in samples with the insertion, it was generally intermediately methylated (average methylation 42.9%) (Fig. 5A). Conversely, the L1 insertion with a hypomethylated promoter (average methylation level 46.3%) at Chr 2: 108,834,214 was associated with a reduction in methylation of its upstream region (from 94.2% to 74.7%) (Fig. 5B; Supplemental Fig. S5B).

Figure 5.

Several examples of TE insertion have been observed to alter the methylation levels of their neighboring regions. (A) The violin plot illustrated the methylation levels of the L1 promoter at Chr 3: 85,527,420. It compared the upstream CpG methylation of the insertion allele and the empty allele. Each dot represented the methylation percentage of a specific CpG site in a particular sample. A t-test was conducted to compare the matching upstream methylation levels, and the P-value is indicated on the figure. (B) The violin plot of L1 promoter methylation at Chr 2: 108,834,214 is similar to A. The upstream CpG methylation of the insertion allele and the upstream methylation of the empty allele were shown. Each dot represented the methylation percentage of a specific CpG site in a particular sample. A t-test was performed and labeled accordingly. (C) The violin plot of Alu methylation at Chr 1: 40,059,013 is similar to A. The downstream CpG methylation of insertion allele and the downstream methylation of empty allele were shown. Each dot represented the methylation percentage of a specific CpG site in a particular sample. A t-test was performed and labeled accordingly. (D) The WashU Epigenome Browser view of a hypermethylated Alu insertion within a hypomethylated valley that was present in HG00621 but absent in HG00741 characterized by WGBS, HiFi, and ONT methylation. All the methylC tracks displayed methylation percentages ranging from 0% to 100%; open chromatin peaks derived from ATAC-seq were shown as blue boxes. (E) Zoom-in view of the methylation of the same Alu element from D. The boxed three CpG sites within 200 bp of the insertion site underwent a change in methylation status, transitioning from hypomethylation to hypermethylation in association with the Alu insertion.

Open in new tabDownload PowerPointLink to figure

The impact of TE insertions on nearby genomic sequence can also be appreciated by examining heterozygous insertion in the same individual, for example, the L1 insertion at Chr 3: 85,527,420 in both HG00621 and HG00741 (present only in the paternal haplotype). Both long-read methods identified the CpG site 46 bp upstream of the L1 insertion was hypermethylated (average 81%) but remained hypomethylated in the remaining eight haplotypes without the insertion (average 21%) (Supplemental Fig. S5H). In contrast, this specific polymorphic L1 was reported by Lanciano et al. (2024) to increase upstream CpG methylation in the MCF-7 cell line but was hypomethylated and had reduced proximal methylation in the 2102Ep cell line. Taken together, these results suggest that this specific polymorphic L1 insertion adopts cell type–specific DNA methylation, which may spread to its immediate region upstream.

We confirmed the impact of multiple L1 loci, including the insertion at Chr 3: 85,527,420 we described above. However, there were also instances in which the relationship between methylation differences and L1 insertions was less evident. For instance, the methylation of the empty allele at CpG site Chr 14: 24,523,788, upstream of the L1 insertion at Chr 14: 24,523,704, exhibited significant variability across the HPRC samples. Meanwhile, the same CpG site was intermediately methylated among samples with the insertion allele, and the difference was not statistically significant (P = 0.42) (Supplemental Fig. S5C,D). At CpG site Chr 5: 137,679,380, upstream of the L1 insertion at Chr 5: 137,679,083, whereas the methylation difference between the insertion and empty alleles was significant (P = 0.0056), the high methylation variation of the empty allele made it difficult to conclude that the insertion allele methylation was influenced by the L1 insertion (Supplemental Fig. S5E,F).

Our study identified thousands of Alu insertions, and most of them were hypermethylated. In some cases, their methylation appeared to be influenced by their genomic context (Fig. 3). Can Alu insertions influence the methylation of their nearby regions? We found that although the large majority of Alu insertions seemed to have no impact on nearby DNA methylation, exceptions did exist. For example, polymorphic Alu insertion at Chr 1: 40,059,013 was hypermethylated (average methylation 94%) and so were three downstream CpG sites within 200 bp (average methylation 58.1%). In samples without the Alu insertion, these three CpG sites were hypomethylated (average methylation 12.2%) (Fig. 5C; Supplemental Fig. S5G). A closeup view revealed that this Alu insertion was homozygous in HG00621 but absent in HG00741 (Fig. 5D,E). This Alu was inserted within a hypomethylation valley (Chr 1: 40,058,800–40,060,800, methylation level 4%). In HG00621, the Alu insertion was hypermethylated (methylation level 85%), and methylation level of the hypomethylation valley was elevated to 21% (Supplemental Fig. S5H; Li et al. 2022; Zhuo et al. 2023). The three CpG sites 200 bp 3′ to the Alu insertion site were lowly methylated in HG00741 (18%, 5%, 4%) but highly methylated in HG00621 (89%, 88%, 93%). On the other hand, multiple CpG sites >300 bp downstream from the insertion site seemed to have remained their hypomethylated state (three CpGs, 3% in HG00741, 5% in HG00621) (Supplemental Fig. S5H). We also mapped Assay for transposase-accessible chromatin using sequencing (ATAC-seq) data to the same genomes and found that the reduction of the hypomethylation valley in HG00621 following the Alu insertion coincided with the reduction of open chromatin ATAC-seq peak (Fig. 5D,E). The impact of this Alu insertion to open chromatin suggested its likely influence on the local cis-regulatory landscape.

Genome assembly and scaffolding

The genome assemblies used in this study were produced for the draft human pangenome, combining PacBio HiFi long-read and Illumina short-read sequencing for de novo assembly with parental short-read data used to for phasing of these assemblies. All assemblies are available at GitHub (https://github.com/human-pangenomics/HPP_Year1_Assemblies). Phased contig assembly was performed using hifiasm (Cheng et al. 2021). Postassembly, the RagTag toolset was employed for reference-guided scaffolding, organizing contigs into structured, phased chromosomal assemblies (Alonge et al. 2022). Validation was performed using Bionano optical mapping in order to detect and amend potential assembly inaccuracies (Lam et al. 2012). For each of the 10 assemblies derived from the five diploid samples, we created a whole-genome pairwise alignment BAM file with GRCh38 using minimap2 (Li 2018). All raw and processed sequencing data from the HPRC are available via the the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/) under accession number PRJNA701308.

Whole-genome bisulfite sequencing

The cytosine methylation was estimated by mapping to reference genome GRCh38 or the individual assemblies using the standard WGBS pipeline Bismark as described by Krueger and Andrews (2011) and Lee et al. (2020a). The raw reads are available via the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) with accession GSE261315.

ATAC-seq

ATAC-seq were conducted following standard protocol with paired-end sequencing on Illumina NextSeq 500 (Buenrostro et al. 2015). The sequencing results were analyzed using ATAC-seq integrated analysis pipeline (AIAP) (Liu et al. 2021a).

PacBio HiFi methylation

PacBio single-molecule HiFi circular consensus sequencing (CCS) reads with kinetics tags can be generated using pbccs (https://github.com/nlhepler/pbccs). We then predicted cytosine methylation probability per molecule in all CpG context by applying primrose to PacBio HiFi reads (https://github.com/PacificBiosciences/jasmine). The CpG methylation probability was saved as the ML and MM tags in the BAM file following standard BAM format convention (http://samtools.github.io/hts-specs/SAMtags.pdf). The HiFi methylation BAM files of the 32 samples are accessible at https://s3-us-west-2.amazonaws.com/human-pangenomics/index.html?prefix=submissions/560047a6-6d16-4b0c-aac9-7d0c83e2188e--HIFI-METHYLATION-READS/. Unaligned HiFi BAM files were aligned to human reference genome GRCh38 using PacBio pbmm2 (https://github.com/PacificBiosciences/pbmm2), and CpG methylation percentage on the reference genome were calculated using PacBio pb-CpG-tools (https://github.com/PacificBiosciences/pb-CpG-tools). We ran the methylation percentage calculation using default model “PB-model.”

Nanopore ONT methylation

We downloaded raw sequencing signal of ONT Nanopore (FAST5 files) of the same five samples and performed methylation calling using the basecaller guppy with the methylation model remora (https://github.com/nanoporetech/remora). CpG methylation probability of ONT reads was estimated using Nanopore basecaller Guppy v6 with remora methylation model (v9.4.1; https://community.nanoporetech.com). The CpG methylation probability were saved as the standard ML and MM tags in the BAM format (http://samtools.github.io/hts-specs/SAMtags.pdf). ONT unaligned BAM files were aligned to human reference genome GRCh38 using minimap2 (https://github.com/lh3/minimap2), and the methylation percentage was calculated using modbam2bed (https://github.com/epi2me-labs/modbam2bed).

Methylation comparison

With genome-wide methylation percentage of all three methods calculated, we made a pairwise correlation heatmap and PCA plot and calculated the MAE, 20% concordance, and Pearson's correlation using R (Wickham 2009; Emerson et al. 2013). We also applied DSS to define pairwise DMRs (Feng et al. 2014; Wu et al. 2015).

Long-read phasing

We used HiPhase for HiFi reads phasing (Holt et al. 2024). To determine in each phase block which haplotype is paternal and which is maternal, we downloaded HiFi reads used for constructing paternal and maternal assemblies from hifiasm output gfa files (HG00621, HG00741, HG01952, HG01978, HG03516). We labeled the haplotype paternal if the number of reads labeled with it from the same phase block used for paternal assembly is more than 10 times higher than that used for maternal assembly; otherwise, it is maternal. Once we separated HiFi reads to paternal and maternal with whatshap (Patterson et al. 2015), we aligned them to the individual paternal and maternal assemblies, accordingly using pbmm2 (https://github.com/PacificBiosciences/pbmm2). At last, we calculated the genome-wide methylation percentage on individual assemblies using pb-CpG-tools (https://github.com/PacificBiosciences/pb-CpG-tools).

To phase ONT reads according to the individual assemblies, we created a diploid VCF file for each sample using dipcall (Li et al. 2018; https://github.com/lh3/dipcall) and haptagged all reads using haplotype information from the VCF file. We then separated paternal and maternal reads using whatshap (Patterson et al. 2015) and aligned them to paternal and maternal assembly using minimap2 (Li 2018, 2021), respectively. At last, we converted the BAM file to methylbed files with modbam2bed (https://github.com/epi2me-labs/modbam2bed).

Individual genome visualization

For each pairwise genome alignment between GRCh38 and individual genome, we created a genome-align track file using the pairwise genome alignment BAM file (Zhuo et al. 2023). We built individual genome assembly of both maternal and paternal chromosomes of the five samples (10 assemblies in total) and aligned each of them to the GRCh38 to create genome-align tracks. We also phased HiFi reads using HiPhase using alignment to GRCh38 and assigned haplotypes to paternal/maternal for each phasing block based on a hifiasm assembly construction graph (Cheng et al. 2021; Holt et al. 2024). We then aligned maternal and paternal HiFi reads to the maternal/paternal assembly separately for each of the five samples. Applying the genome-align track we anchored gene, repeat, and CGI annotations from GRCh38 onto the individual genomes. We colored each CpG site by their methylation prediction (red: methylated; blue: unmethylated) per HiFi read using the modbed tracks and displayed piled up methylation level using the methylC track. Both modbed and methylC tracks were mapped to the individual genome. Combining genome-align, methylC, and modbed track files, we created three datahubs on the WashU Epigenome browser representing HiFi, ONT, and WGBS data (https://epigenomegateway.wustl.edu/browser2022/?genome=hg38&hub=https://wangcluster.wustl.edu/~xzhuo/hifi_methylation/hifi.all.json, https://epigenomegateway.wustl.edu/browser2022/?genome=hg38&hub=https://wangcluster.wustl.edu/~xzhuo/hifi_methylation/ont.all.json, https://wangcluster.wustl.edu/~xzhuo/hifi_methylation/wgbs_individual/wgbs.all.json) (Zhou et al. 2014; Li et al. 2022, 2023).

After aligning WGBS reads to the individual assembly as described above, we also plotted the WGBS methylation of nonreference Alu insertion at Chr 11: 59,565,869 using Methylartist (Ewing et al. 2020; Cheetham et al. 2022).

Structural variation and TE annotation

We obtained the structural variations of each sample from HPRC consortium (Liao et al. 2023; https://s3-us-west-2.amazonaws.com/human-pangenomics/index.html?prefix=submissions/B581EBA7-8BDE-4C7C-9DEA-78B99A051155--Yale_HPP_Year1_Variant_Calls/). We merged all the called SVs from the 32 samples using BCFtools merge (Danecek et al. 2021). To further annotate TE insertions among all structural variations, we ran MELT-LRA and PALMER2 on HiFi assemblies and xTEA on the same PacBio HiFi library mapped to GRCh38 (Gardner et al. 2017; Chu et al. 2021; Zhou et al. 2024). To generate a high-quality full-length TE insertion set, we filtered out Alu insertions <300 bp, SVA < 1000 bp, L1 < 5500 bp, and LTR5 < 900 bp. For Alu, L1, and SVA insertions, we included MEIs called by at least two of three methods. MELT-LRA did not call HERVK insertion, and we manually verified 25 out of 34 HERVK insertions identified by either xTEA or PALMER2. We summarized MEI caller output in Supplemental Table S1. We then intersected the coordinates of TE insertions with SV and annotated which SVs are TE insertions (Supplemental Table S2).

We developed a script to extract nonreference insertions. Using a BED file of annotated insertion sites, it can extract inserted and soft-clipped sequences from raw reads and generate a new unmapped modBAM file of inserted sequences with base modification tags (MM and ML) (bam.extractInsertion.py). We can then align the unmapped modBAM file to their respected TE consensus (L1HS, AluY, SVA_F, LTR5_Hs, and HERVK-int from Dfam database) to create CpG methylation pileup results (Wheeler et al. 2012; Hubley et al. 2016).

To calculate average methylation of certain regions, we extracted the methylation prediction value of each CpG site of each read from BAM files with methylation tags with a Python script (Python bam.Mmtag.regions.oo.py -b <bam> -r <insertion.bed> -o <ins.methylC.reldist.txt> -f <bed>) (see Code availability).

To evaluate the relationship between CpG methylation and indels, we aggregated methylation prediction of all 32 samples and calculated the average CpG methylation of each CGI (Quinlan and Hall 2010).

Reference TE methylation

We extracted 52,802 AluY, 2746 AluYa5, 36 AluYa8, 2111 AluYb8, and 301 AluYb9 from GRCh38 using RepeatMasker annotation. Similarly, we extracted 124 SVA_E, 331 SVA_F, 481 LTR5 as reference TEs. We also downloaded 146 active L1 on GRCh38 from L1 Base2 (Penzkofer et al. 2017). We then take their DNA sequences from the GRCh38 reference and their corresponding methylation from piled up bedMethyl files and generated new modBAM files to bind methylation percentage to the sequence as the modified base (bam.creation.py). At last, we mapped the new BAM file to TE consensus using minimap2 to align the CpG methylation percentage to the TE consensus (L1HS, AluY, SVA_F, LTR5_Hs, and HERVK-int) (Wheeler et al. 2012; Hubley et al. 2016).

Reference or nonreference insertions are relative to the reference genome chosen for alignment, and nonreference insertions can be aligned to the genome if we chose the individual assembly as the “reference.” Therefore, the sample process can be applied to extract “non-ref” insertions on the GRCh38 after aligning long reads to the individual assembly.

Insertion and flanking methylation calculation

We selected nonreference insertions from annotated SVs and used a custom Python script (see Code availability) to calculate the average methylation of the insertion and flanking regions (python read.methylation.mean.py -i <input> -o <output> -s <start> -l <length> -c).

CGI methylation and indel frequency calculation

The unmasked CGI annotation file was downloaded from UCSC (https://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/cpgIslandExtUnmasked.txt.gz). For each CpG site on GRCh38, we calculated the median methylation percentage of the 32 HiFi methylation. We then calculated the average methylation percentage of each CGI using the median methylation.

After separating CGIs to hypermethylated and hypomethylated ones using 50% methylation as the cutoff, we calculated the indel frequency of different sizes within CGIs. To create the background indel frequency in the genome, we shuffled CGIs in the genome and repeated the calculation.

We also repeated the same calculation using the same methylation data and indel callset from The 1000 Genome Project and HGSVC project (The 1000 Genomes Project Consortium 2015; Ebert et al. 2021).

Code availability

Custom scripts are available at GitHub (https://github.com/xzhuo/modbamUtil) and as Supplemental Code.