Development and validation of a Bayesian multiread allocation algorithm

We initially developed our SmartMap algorithm and software for application in ChIP-seq using a set of internally calibrated ChIP-seq (ICeChIP-seq) datasets. These datasets were previously generated by our lab and, with one exception, were previously published as components of past studies [29,30]. We chose to use ICeChIP-seq datasets because the included internal standards allow for computation of antibody specificity and for normalization to calculate the histone modification density (HMD), or the absolute proportion of nucleosomes at a given genomic locus bearing the targeted histone modification. These additional factors which we can compute using ICeChIP-seq datasets afford us additional points of quantitative comparison to assess differences between uniread and SmartMap analyses. However, this tool is not designed solely (or even primarily) for use with ICeChIP-seq datasets; the SmartMap algorithm does not make special use of the internal standards. Rather, this software is designed to be usable for NGS workflows more broadly.

The workflows for uniread analyses (typical of ChIP-seq) and our SmartMap analysis are shown in Fig 1A. For both analyses, the immunoprecipitation (IP) and MNase-seq Input sequences are aligned to the appropriate reference genome and are filtered to select for properly mapped reads in a proper pair. At that point, the two methods diverge. In the uniread analysis, which represents our published analysis pipeline for ICeChIP-seq data [29–31], any reads that don’t align uniquely are discarded, and the remainder are used to compute genome-wide read depth in the IP and Input, fold-change, and (if internal standards are present) HMD.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Summary of the SmartMap analysis workflow and algorithm.

(A) Flowchart outlining the workflow for traditional ChIP-seq (or ICeChIP-seq) analysis [29–31] utilizing only unireads (left, green) vs. the workflow for SmartMap analysis utilizing multireads with an iterative Bayesian reweighting algorithm (right, blue). (B) Schematic showing the Bayesian reweighting algorithm utilized in the SmartMap analysis. Each mapping associated with a read is assigned a weight such that the weight is greater for those mappings associated with loci of greater map weight density. For more detailed description of the algorithm, see Methods.

https://doi.org/10.1371/journal.pcbi.1008926.g001

In the SmartMap analysis, however, rather than discarding ambiguously mapped reads, we instead feed our alignments into our iterative Bayesian reweighting algorithm, outlined in Fig 1B. Our algorithm, like other alignment weighting algorithms [2,35,36,39], is motivated by the assumption that regions with more alignments are more likely to be the true source of an multiread than those with fewer alignments. In addition, like BM-Map [36], SmartMap utilizes both paired-end sequencing information and alignment quality in making these assessments. Accordingly, our tool first assigns each alignment a weight proportional to its alignment quality, computed from the alignment software output. We then iteratively reassign weights to each alignment of each read; alignments with higher alignment quality and more overlapping alignments are assigned higher weights, and those alignments with lower quality and fewer overlapping alignments are assigned lower weights (Fig 1B). After the specified number of reweighting cycles, the resulting weights are used to compute the read depth for the IP and the Input genome-wide, which can then be used to compute fold-change or, if applicable, HMD in a similar manner as the uniread analysis. For computational efficiency, we use binary-indexed (Fenwick) trees to store genomic coordinates and associated alignment weights, much like the previously described CSEM [35]. Our implementation of these binary-indexed trees is modified to enable use of paired-end sequencing reads and, if needed, operate in a strand-specific manner.

To test this method, we created a set of simulated 50bp paired-end sequencing reads from a defined set of randomly selected genomic loci (the “true origin” loci) and used the simulated dataset to conduct uniread and SmartMap analyses (Fig 2A). The read simulation tool produces reads with “sequencing” error and also includes coverage at off-target loci to better represent the noise and off-target capture inherent in a biological experiment. Notably, the simulation enabled us to obtain the true distribution of reads (the Gold Standard), allowing us to compute the error associated with each analysis method (Fig 2A). This is particularly important because we wish to avoid over-refining the multiread weights with our inferential analysis; accordingly, the Gold Standard dataset allows us to evaluate the accuracy of reweighting.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Characteristics of validation dataset.

(A) Schematic outlining the workflow to validate and optimize SmartMap. A set of six million randomly selected 200bp loci were used to simulate paired end reads. The true read depth distribution was then compared to both uniread and SmartMap analyses, with each analysis conducted in both “scored” and “unscored” modes, per Methods. (B, C) Number of (B) alignments or (C) reads vs. number of alignments per read for the validation datasets. (D) Mean absolute error of read depth at true origin loci in SmartMap scored mode vs. number of reweighting iterations (E) Genome browser view showing the read depth in the (top) uniread, (center) SmartMap (0 iterations), and (bottom) SmartMap (1 iteration) datasets of an example locus.

https://doi.org/10.1371/journal.pcbi.1008926.g002

We were particularly interested in the ability of SmartMap to recover read depth at regions of differing mappability. To investigate this relationship, we used the UMAP50 score as a measure of read mappability. The UMAP50 score for a given genomic coordinate is computed as the proportion of the 50mers covering the genomic coordinate of interest that are unique in the set of all 50mers from the genome [25]. For example, if the sequences of two of the fifty 50mers containing the genomic coordinate of interest were non-unique across the genome of study, then the UMAP50 score would be 48/50, or 0.96. As such, a genomic coordinate with a UMAP50 score closer to 1 is uniquely identified by a greater proportion of the 50mers spanning it than is a coordinate with a lower UMAP50 score, and a higher UMAP50 score can thus be interpreted as a more easily mappable region. Many of the true origin loci had low mappability scores (S1A Fig), with the distribution of mappability scores being similar to that of the human genome at large (S1B Fig), making this dataset useful for validating the SmartMap algorithm.

The first step of our analysis was to align the simulated 50bp paired-end reads to the genome. We used Bowtie2 [40] with a maximum of 51 alignments reported per read and charted the distributions of the number of alignments per read (Fig 2B and 2C). Notably, we observed that there were many reads that did not uniquely align to the genome; approximately 17.1% of the simulated reads mapped to more than one locus (Fig 2B and 2C and Table 1).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Alignment statistics for the datasets used in this study.

https://doi.org/10.1371/journal.pcbi.1008926.t001

Our first goal was to determine the optimal number of iterations to use for our SmartMap analyses. To test this, we computed the mean absolute error of SmartMap read depth at the true origin loci with varying numbers of reweighting cycles, as compared to the Gold Standard read depth. Surprisingly, we found that the lowest error occurred after only one reweighting cycle (Fig 2D), with genome browser views showing refinement of peak structure (Fig 2E), which is particularly important given the importance that has been placed on peak breadth [41]. This stands in stark contrast with previous works, which have used up to 200 iterations of reweighting [35]. Our analysis here, however, shows that may be suboptimal, suggesting that applying Bayesian alignment reweighting more than once may over-refine the data.

We wanted to explore whether these increases in mean absolute error were systematic or driven by random “overshoot” of weight at each locus. In the former case, we might expect to see that the true origin loci would either show systematic increases or decreases in read weight with greater numbers of reweighting cycles. In the latter case, we would expect that the changes to each weight might increase or decrease by too much in the initial iteration, which would present as random, relatively unbiased errors.

To distinguish between these two possibilities, we conducted two analyses. First, we computed mean error of weights at the true origin loci (S2A Fig) rather than the mean absolute error (Fig 2D). If there was a systematic erroneous increase or decrease in the average read depth of each locus, then we would observe a corresponding increase or decrease in mean error with more iterations, respectively. However, what we instead observe is that the mean error is relatively stable from iterations 2–8 (S2A Fig), suggesting that the marked increase in mean absolute error with increasing iterations is not primarily caused by systematic erroneous increases or decreases in locus weight depth. Put differently, it does not appear that the true loci are incorrectly and systematically “pulling in” or “pushing out” read depth with each reweighting cycle.

Second, we explored the possibility that the reweighting “overshoots” the weight adjustment for reads at random. If this was the case, then we would expect that the errors would increase relatively randomly, with both positive and negative errors. Indeed, this is what we observe in our analysis of mean error by iteration (S2A Fig). In addition, we would predict slowing the rate of weight adjustment with each cycle would decrease the amount of overshoot and thereby lead to a lesser increase in error. To test this, we introduced a tunable reweighting rate parameter such that the weights could be changed less with each reweighting iteration. When applying a reweighting rate of 0.25 (wherein the weights only change by 25% as much in normal SmartMap analysis), we found that the mean absolute error was markedly more stable after one iteration (S2B Fig). Indeed, after two cycles of standard SmartMap, the mean absolute error exceeds that of the dataset with no reweighting (Fig 2D); by contrast, with eight cycles of SmartMap with a reweighting rate of 0.25, the mean absolute error is considerably below that of the iteration 0 dataset and comparable to the minimum mean absolute error after one iteration (S2B Fig). This suggests that the increase in error with increasing iterations observed with standard SmartMap may be due to “overshoot” of reweighting, which compounds in magnitude with further reweighting. Interestingly, we found that the mean absolute error with one iteration of standard SmartMap analysis was on par with (and even slightly lower than) that of the slow-reweighting dataset (Figs 2D and S2B), suggesting that this potential overshoot error may not be too detrimental after only one iteration of reweighting. By command line switch, these two algorithms are both available in the SmartMap software.

After determining the optimal number of reweighting cycles, we then compared the SmartMap and uniread analyses of our simulated datasets (Fig 2A). To determine the relative impact of using alignment quality for multiread analysis, we ran SmartMap in both scored and unscored modes. All the SmartMap analyses had greater read depth (and were closer to the Gold Standard dataset) at true origin loci than the corresponding uniread analyses (Fig 3A). Interestingly, the increases in read depth were not uniform across the entire set of loci; indeed, approximately 70% of the true origin loci saw no excess read depth, defined as the difference between SmartMap and uniread read depths (Fig 3B). This is similarly observed in the QQ plot comparing uniread and SmartMap analyses; a shoulder is seen at low uniread depth, with the plot converging onto a slope of unity at higher read depths (Figs 3C and S2C), suggesting that the gains in read depth were primarily at regions of low uniread depth. While SmartMap does not fully recover the read depth of the Gold Standard at the low end of the QQ plot, the shoulder is nonetheless much less prominent than with the uniread analysis (S2D Fig), indicating considerably greater depth recovery. Consistent with that observation, the uniread analyses and the SmartMap analyses both performed well at highly mappable regions, with read depths approximately at the level of the Gold Standard (Fig 3D). However, at regions of lower mappability, the SmartMap analyses recovered a markedly greater proportion of the read depth than did the uniread analyses (Fig 3D and 3E). As expected from prior analyses (Fig 2D), the SmartMap analyses with one iteration of reweighting recovered greater read depth than those with no reweighting (Fig 3D and 3E). Importantly, though they performed similarly at regions of lower mappability, the SmartMap scored analyses recovered greater read depth than their unscored counterparts at regions with moderate mappability scores (Fig 3D and 3E).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. SmartMap and uniread analyses of the validation dataset.

Iteration 0 and iteration 1 refer to SmartMap analysis with 0 and 1 iterations of reweighting, respectively. Scored and unscored refer to whether alignment score was considered in analysis, per Methods. Dashed lines are presented for readability of overlapping curves rather than discontinuities in data throughout this figure. (A) Quantile plot of read depth at the true origin loci, with Gold Standard dataset and analysis conducted in (left) scored mode or (right) unscored mode. (B) Quantile plot of excess read depth in SmartMap datasets relative to corresponding uniread dataset at true origin loci in (left) scored mode and (right) unscored mode. (C) QQ plot of read depth at true origin loci in the SmartMap (1 iteration) scored dataset vs. uniread scored dataset. Color scale represents percentile of each point, from 1^st to 99^th percentiles. (D-E) Median (D) read depth or (E) excess read depth vs. mappability score (UMAP50) [25] of the true origin loci. (F-G) Average read depth (F) at true origin loci and (G) outside true origin loci. (H) Mean absolute error of read depth at true origin loci for each dataset, with Gold Standard as the reference point. (I) Mean proportion of alignments intersecting with the true read of origin for each weight after SmartMap with no reweighting (green) and one iteration of reweighting (red) in scored mode. Dashed line represents line with slope of unity. (J) Mean weighted overlap proportion score between alignments intersecting the true read of origin and the true read locus for each weight after SmartMap with no reweighting (green) and one iteration of reweighting (red) in scored mode. Weighted overlap proportion score is meant to represent the proportion of a read’s weight that maps to the correct location due to a particular alignment and is computed as a weighted geometric mean of the proportion of the alignment covered by the true read and the proportion of the true read covered by the alignment.

https://doi.org/10.1371/journal.pcbi.1008926.g003

Genome-wide, SmartMap analyses had lower on-target read depth than the Gold Standard dataset but were still able to recover greater depth at the on-target loci than corresponding uniread analyses (Fig 3F). Similarly, the SmartMap analyses had marginally higher off-target read depth than the Gold Standard and uniread datasets (Fig 3G); however, the increased off-target depth relative to uniread datasets can be explained by the overall lower read depth in the uniread datasets (S3A Fig). Consistent with the notion that improved priors enhance Bayesian predictions, the unscored SmartMap analyses had lower on-target and higher off-target read depth than the corresponding scored analyses (Fig 3F and 3G), and the no-iteration SmartMap analyses had similarly lower on-target and higher off-target read depth than their one-iteration counterparts.

As another metric to evaluate each analysis, we conducted MACS2 peak calling on each dataset and assessed the degree to which they overlap. The SmartMap analyses had similar (albeit slightly higher) base pair coverage with called peaks relative to the Gold Standard dataset and considerably higher coverage on called peaks than the uniread analyses (S3B Fig), consistent with the genome-browser views that suggest a similar pattern of peak boundary sharpening (S1 Fig). As a measure of sensitivity, we computed the proportion of the Gold Standard peaks that were covered by SmartMap or uniread peaks (S3C Fig). Conversely, to measure specificity, we computed the proportion of SmartMap or uniread peaks that were covered by Gold Standard peaks (S3D Fig). As expected, there was considerably lower coverage by the uniread datasets than the SmartMap datasets, and the one-iteration SmartMap analyses had very slightly lower coverage over the Gold Standard peaks than the no-iteration analyses (S3C Fig). However, the one-iteration analyses were better-covered by Gold Standard peaks than were their no-iteration counterparts (S3D Fig). Together, these data suggest that SmartMap analyses with one iteration of reweighting have a marked increase in specificity relative to the no-iteration analyses at the expense of a slight decrease in sensitivity.

We then evaluated the overall mean absolute error of read depth at the true origin loci relative to Gold Standard. The uniread analyses had the highest average mean absolute error, with all SmartMap analyses outperforming all uniread analyses (Fig 3H). The scored SmartMap analyses also all had lower error than did the unscored analyses, and the one-iteration SmartMap analyses slightly outperformed the no-iteration analyses (Fig 3H). The error in all datasets tended to primarily be concentrated at regions of lower mappability (S3E Fig). Interestingly, though SmartMap with one iteration had lower mean absolute error overall, the no-iteration modality had slightly lower mean absolute error at true origin loci of lower read depth (Figs 3H and S3F). The reason for this difference is not clear; across all read depth classes, the one-iteration analyses had slightly less negative mean error, suggesting that there wasn’t a large-scale difference in over- or underweighting after iteration as a function of read depth (S3G Fig). With that said, we feel it is important to contextualize these results; these differences between the no-iteration and one-iteration analyses are small in magnitude and are comparatively dwarfed by the differences between SmartMap and uniread analyses (Figs 3H and S3E–S3G). Accordingly, though there may be small differences between the no-iteration and one-iteration SmartMap analyses, the one-iteration analyses still performed better in aggregate, and both of these scored SmartMap analyses consistently outperformed their unscored or uniread counterparts.

The above analyses all focused on validating SmartMap from the perspective of the total read depth across a set of genomic intervals. However, given that we had a Gold Standard dataset listing the true positions of each read, we also wished to evaluate whether our reweighting method could improve the estimated of the probability that an alignment was properly mapped–and, by proxy, improve the MAPQ score estimate for each alignment. Without reweighting, the probability of correct alignment ranged from 0–0.67 and 1, with no alignments with correct alignment probability between 0.67 and 1. One iteration of SmartMap reweighting expanded the spectrum of possible alignment weights to the full range of 0–1. Without reweighting, the weight of alignments did not correlate well with the proportion of alignment intersecting the true genomic position, with many large deviations seen from linearity (Fig 3I). By contrast, though one iteration of reweighting still showed some deviations from linearity by this analysis, the weight of alignments more closely concorded with the proportion of the alignments intersecting the true read origin (Fig 3I). This suggests that by this measure, SmartMap reweighting improved the estimates of the probability that the alignment intersects with the true genomic position of the corresponding read. Similarly, we compared the weighted proportion of overlap between the true read positions and any intersecting alignments as a function of alignment weight. This is meant to represent the proportion of a read’s weight that is mapped to the correct location due to a given alignment and incorporates both the confidence of the alignment selection (i.e. the weight) and the overlap of the alignment with the true origin of the read. In both the no-reweighting and one iteration analyses, the overlap proportion score was closely linearly related to the alignment weight, though the reweighted analysis showed a slightly smoother curve with fewer marked deviations from linearity (Fig 3J). This is roughly expected, as the overlap proportion score is itself a function of weight; however, this analysis is comforting insofar as it shows that the SmartMap reweighting does not markedly inflate or deflate the expected weight contribution of a given alignment to a proper intersection with the true origin. Similarly, we find that the unweighted overlap proportion of alignments with the true origin of the read is roughly constant near one for both the no-iteration and one-iteration datasets, though again, the one-iteration SmartMap analysis reduces the deviations from this level (S3H Fig). These analyses suggested that in addition to improving measurement of read depths in aggregate, the SmartMap reweighting procedure can also improve the estimates of correct alignment for individual reads and alignments.

The biological ChIP-seq and MNase-seq datasets presented in the remainder of this work used 50bp read lengths or shorter, which is why we used 50bp read lengths for our simulated dataset. However, in recent years, 100bp read lengths have become commonplace, and indeed, the ENCODE datasets we present later in this work employed paired-end 100bp NGS. As such, we examined the degree to which SmartMap can improve recovery of sequencing depth with longer reads by conducting a similar analysis as the above with a similarly simulated dataset employing 100bp paired-end reads. For facile comparison to the other figures and analyses in this work, we have continued to use the UMAP50 score as our mappability score. This choice is in spite of the fact that UMAP50 measures mappability by 50mers rather than 100mers and will thus underestimate mappability by 100bp reads. Because we are using this metric, regions with lower mappability scores will often be more easily mapped than the score would indicate, blunting differences between SmartMap and uniread analyses. As such, our analyses using the UMAP50 score will offer a very conservative view at the impact of SmartMap analysis on read depth recovery and error.

Despite this conservative choice of mappability score, we still see that SmartMap analysis improves sequencing depth recovery nearly as well with 100bp reads as it does with 50bp reads. The simulated dataset with 50bp reads shows a 13.9% increase in analyzable reads due to the high number of multireads (2B and Fig 2C and Table 1); the simulation with 100bp reads shows a 13.6% increase in analyzable reads and a similar proportion of multireads (S4A and S4B Fig and Table 1). Along the same lines, the two simulations increase read depth over similar proportions of the genome (Table 2). Over the true origin loci, much like the 50bp simulation, the 100bp simulated dataset shows an increase in read depth on quantile plots (S4C Fig) under SmartMap analysis, with this increase in read depth primarily occurring at regions of low UMAP50 mappability score (S4D Fig), conservative though this measurement of mappability is. Much like the 50bp simulated datasets, the increases in read depth under SmartMap analysis are primarily seen at regions of low mappability and low uniread read depths; QQ plots comparing the uniread analysis with the SmartMap or Gold Standard show a shoulder at low uniread depths, with the plot converging onto a slope of unity at higher uniread depths (S4E and S4F Fig). It should be noted that, as with the 50bp simulated dataset (S2D Fig), the SmartMap dataset still fails to fully recover read depth as compared to Gold Standard with 100bp reads (S4G Fig). Nonetheless, the SmartMap analysis still shows considerably lower mean absolute error than does the uniread analysis at true origin loci (S4H Fig), with this decrease in error being particularly prominent at regions with lower UMAP50 mappability scores (S4I Fig). In total, these analyses suggest that even for datasets employing 100bp paired-end sequencing reads, multiread analysis still has nearly undiminished importance and that SmartMap can still markedly improve read depth recovery while decreasing overall error.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Analysis of reads across genomic windows.

https://doi.org/10.1371/journal.pcbi.1008926.t002

The above analyses all restricted Bowtie2 to report a maximum of 51 alignments for computational efficiency. Subsequently, only those reads aligning to fewer than 51 alignments were used for SmartMap analysis. However, this practice excluded more than 7 million reads (Table 1), likely including reads that map to the most highly repetitive regions of the genome. Notably, this is a restriction on alignment itself, not SmartMap; there’s no reason that SmartMap would inherently be unable to handle greater numbers of alignment. Nonetheless, to evaluate the impact of this restriction on the SmartMap datasets, we reanalyzed our simulated 50bp read length dataset with a maximum of 101 alignments per read (hereafter, the k101 dataset) and compared it to the previous analysis (the k51 dataset). To our surprise, the two analyses were highly similar despite the near-doubling in the maximum-alignments threshold in the former dataset. The increase in the number of analyzable reads was nearly identical between the two analyses (S5A and S5B Fig and Table 1), with similar increases in depth over genomic windows (Table 2). At the true origin loci, the SmartMap read depths in both the k51 and k101 datasets were very similar at the level of read depth (S5C and S5D Fig). Mean absolute error relative to the Gold Standard was actually very slightly lower in the k51 dataset, though they were quite similar in magnitude compared to the uniread dataset (S5E and S5F Fig). Even specifically examining repetitive elements, read depth was very similar between the k51 and k101 SmartMap analyses at all repeats (S5G Fig), LINEs (S5H Fig), SINEs (S5I Fig), and Alu elements (S5J Fig), closely approximating the Gold Standard read depth in both cases. Accordingly, though there is still a large proportion of reads that mapped to still greater numbers of loci, we find that at the range we have tested, the SmartMap analyses are robust to differences in maximum-alignments reporting thresholds and that there is little practical difference between restricting datasets to a maximum of 51 or 101 alignments per read besides the additional time and storage space needed for the latter.

To be sure, the reweighting used for SmartMap is not without concerns. In particular, one of the potential problems for SmartMap is the existence of high-signal regions, which can show falsely high read depth in NGS experiments due to sequencing or alignment error [42]. If there are regions of falsely high weight, then those regions could be skewed by the SmartMap reweighting algorithm to report even greater weights, thus exacerbating these artifactually high signals. To assess the degree to which these regions represent an issue for SmartMap, we computed the number of genomic windows with more than 60, 70, 80, or 90 reads in our simulated datasets (Table 3). We used these benchmarks as rough thresholds for defining high-signal regions because the Gold Standard dataset had a maximum read depth across a genomic window of approximately 83 reads. Notably, the Gold Standard did not require sequencing or mapping and should thus not be susceptible to these high-signal artifacts. Unfortunately, one iteration of SmartMap reweighting did increase the proportion of high-signal regions considerably; there were fewer than 600 genomic windows with an average depth of more than 70 in the Gold Standard, Iteration 0, and Uniread datasets, compared to more than 10,000 in the SmartMap dataset with one reweighting cycle. It’s important to note that these regions represent a very small proportion of the genome; only 0.066% of the genomic windows had more than 70 reads on average, and even fewer had more than 80 or 90 reads (Table 3), leaving considerably more than 99.9% of the genome as not having abnormally high-signal attributable to SmartMap. In contrast, almost 15% of the genome is hidden from uniread analysis (Table 2). Nonetheless, we feel it is fair to say that the reweighting algorithm used for SmartMap will increase the weights of multiread alignments at high signal regions, which can exacerbate artifactually high read depths.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Analysis of high-depth regions under SmartMap analysis.

https://doi.org/10.1371/journal.pcbi.1008926.t003

Even so, on the whole, these analyses suggest that SmartMap recovers read depth at a large set of loci that would otherwise be missed by the uniread analyses and that of the SmartMap analyses, one iteration of reweighting with use of alignment scores largely outperforms the other modalities. Accordingly, for the remainder of this work, we use SmartMap analysis with one iteration in scored mode as our default SmartMap method.

Utilizing SmartMap on MNase-seq and ChIP-seq datasets

Having validated our method on the simulated dataset, we turned to the biological samples. We deployed a total of 21 datasets derived from three different organisms for our analysis (Table 1). Of these datasets, six were control ICeChIP Inputs, generated by MNase-seq [29,30], 11 were ICeChIP-seq IP datasets, two were ATAC-seq datasets, and two were RNA-seq datasets. After alignment, the samples showed a 13–50% increase in the number of usable reads for SmartMap analysis relative to uniread (S6 and S7 Figs and Table 1).

To evaluate the impact of our algorithm on the ICeChIP-seq datasets, we first conducted SmartMap and uniread analysis on each of the Input datasets and computed the average read depth on 200bp genomic windows. As with the simulated dataset, the SmartMap analyses of the Inputs had increased read depth relative to the uniread datasets (Figs 4A and S8A), with markedly greater depth in the SmartMap analysis at windows of lower mappability (Figs 4B and S8B). Similarly, this excess read depth was not distributed across all reads, but rather, was concentrated onto 17–35% of windows (Figs 4C and S8C and Table 2), primarily at regions of lower mappability (Figs 4D and S8D). The QQ plots of the SmartMap vs. the uniread read depths showed a shoulder at low uniread depth (Figs 4E and S8E), again suggesting that the increase in read depth from the SmartMap analysis is primarily at loci where the uniread analysis performs poorly. This difference in the distributions of read depths further comments on the importance of analyzing multireads.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. SmartMap and uniread analyses of ICeChIP-seq input depth.

All analyses conducted on 200bp genomic windows for the Inputs defined in Table 2. (A) Quantile plot of read depth for SmartMap and uniread analyses. (B) Median read depth vs. mappability score (UMAP50) for SmartMap and uniread analyses. (C) Quantile plot of excess read depth in SmartMap relative to uniread analysis. (D) Median excess read depth vs. mappability score (UMAP50). (E) QQ plot of read depth in SmartMap vs. uniread analysis. Color scale represents percentile of each point, from 1^st to 99^th percentiles. Dashed line represents line with slope of unity. (F) Quantile plots of depth-normalized log ratio of read depths of biological input replicates under SmartMap and uniread analysis. Graph breaks are present on both the upper and lower ends of the graphs. (G) Mean absolute depth-normalized log ratio for the comparisons presented in panel F.

https://doi.org/10.1371/journal.pcbi.1008926.g004

With our Input datasets, we could also examine the reproducibility of the MNase-seq experiments under uniread and SmartMap analyses. There were three biological replicates of Input in mESC E14 cells (AR7 Replicate 1, AR7 Replicate 2, and AR9), and two biological replicates of Input in K562 cells (AR16 and AR17). For all loci with nonzero read depth, we computed the depth normalized log ratios of reads in a pairwise manner for biological replicates, shown as quantile plots in Fig 4F. These plots are highly similar under SmartMap and uniread analyses across all pairwise comparisons (Fig 4F). Accordingly, the average magnitudes of these ratios are similar between the two analyses–and indeed, are slightly lower in the SmartMap datasets (Fig 4G). This suggests that the two modalities show highly similar estimates of reproducibility of data between biological replicates of MNase-seq.

Having examined the Input datasets, we then used the ICeChIP-seq datasets to compute histone modification densities (HMD) across 200bp genomic windows with both uniread and SmartMap analyses. Interestingly, we noted that the mean HMD was quite similar between the SmartMap and uniread datasets across a broad range of mappability scores (Figs 5A and S9A). However, the median HMD of those same datasets were divergent, with the SmartMap analyses having considerably higher median HMD across bins of low mappability than the uniread analyses (Figs 5B and S9B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. ICeChIP-seq histone modification density in SmartMap and uniread analyses.

All analyses conducted on 200bp tiled genomic windows. (A-B) (A) Mean or (B) Median HMD vs. mappability score (UMAP50) for SmartMap and uniread analyses. (C-D) Scatterplots of (C) specificity or (D) log specificity for uniread vs. SmartMap analyses. Specificity is measured as the enrichment of each on- or off-target internal standard nucleosome as a percentage of on-target enrichment.

https://doi.org/10.1371/journal.pcbi.1008926.g005

The difference between mean and median HMD may be attributable to the fact that HMD is a scaled-version of fold-change of IP over Input. We attribute the median HMD divergence to sparser distribution of read depth in the uniread dataset at lower mappability scores (Fig 4B). As such, there are fewer regions with nonzero read depth in both the IP and Input. The result of this mismatch in read distribution is that more regions have an apparent HMD of zero under uniread analysis. That the mean HMDs are similar between the two analyses suggests that the ratios of the total read depths in IP over Input are similar between SmartMap and uniread analyses. Together, these data suggest that the SmartMap analyses preserve the overall HMD across a wide range of mappability scores while also enabling measurement of HMD at a broader range of loci than do uniread analyses.

One of the major benefits of using ICeChIP-seq data is the ability to measure antibody specificity [29–31]. In ICeChIP, internal standards bearing a variety of different histone modifications can be simultaneously spiked into an experiment, and the relative pulldown efficiency of each modification can be quantified as a proportion of the target to measure the off-target binding of the antibody. We wished to determine whether the SmartMap analyses would yield similar specificity estimates as did the uniread analyses. First, we found that the ratio of the reads from the on-target nucleosome in the IP over the Input was highly similar between the uniread and SmartMap analyses (Table 4). Moreover, the scatterplots of specificity (Figs 5C and S10A) and logarithm of specificity (Figs 5D and S10B) under each modality show slopes close to unity and high coefficients of determination. This further shows that the specificity measurements in SmartMap and uniread analyses are highly similar in both an absolute (Figs 5C and S10A) and a relative (Figs 5D and S10B) sense.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Analysis of ICeChIP Calibrant Barcodes.

https://doi.org/10.1371/journal.pcbi.1008926.t004

Extending the utility of SmartMap to ATAC-seq and RNA-seq

We also found that SmartMap could be applied to ATAC-seq data to obtain more global measurements of chromatin accessibility. To demonstrate this, we used two replicates of K562 ATAC-seq data, originally generated by the Snyder Lab as part of the ENCODE Consortium [15]. As with the ICeChIP-seq datasets, we found that SmartMap analysis could utilize 17–20% more reads than uniread analysis (Table 1); this increased read depth was primarily concentrated at 20–30% of the genome (S11A–S11C Fig and Table 2), particularly those loci with low mappability scores (S11D–S11F Fig). SmartMap and uniread analyses also showed similar levels of reproducibility between the two isogenic replicates, though SmartMap showed slightly lower reproducibility between the two datasets than did the uniread analysis (S11G and S11H Fig). These data suggest that SmartMap is also useful for ATAC-seq data and can reveal accessible regions of the genome at poorly mappable loci that would have been missed by uniread analysis alone.

In addition to the MNase-seq, ChIP-seq, and ATAC-seq datasets, we also sought to apply our SmartMap analysis to RNA-seq experiments. Specifically, we analyzed two replicates of K562 bulk RNA-seq data, originally generated by the Gingeras Lab as part of the ENCODE Consortium [15]. Our SmartMap RNA-seq analyses showed that for each replicate, relative to uniread analysis, there was a 20–35% increase in usable reads (Table 1) concentrated into a minority of distinct Refseq genes (S12A–S12C Fig and Table 2). The reproducibility of the two datasets was also similar between the SmartMap and uniread analyses, though as with the ATAC-seq data, the SmartMap analysis showed marginally lower reproducibility between the two RNA-seq experiments than did the uniread analysis (S12D and S12E Fig). With that said, these differences in read depth are relatively minor in magnitude, especially when normalized to differences in read depth in the SmartMap and uniread analyses. Given the other concerns with using this particular multiread allocation algorithm in gapped reads or spliced transcripts (as noted in the Discussion), it is likely that SmartMap is not optimally configured for use in RNA-seq analysis.

SmartMap drives new biological insights about repetitive DNA elements

With this method, we sought to better explore the role of histone modifications at repetitive regions. Traditionally, the epigenetic profile of repetitive elements is viewed in light of the “genome defense” hypothesis, which suggests that regulation of repetitive elements (and particularly transposable elements) serves to silence the elements and thereby prevent transposition [43]. Consequently, much previous work on this topic has primarily pointed towards repetitive elements being enriched with heterochromatin-associated modifications such as H3K9me2 [44], H3K9me3 [9,43,45–47], and H3K27me3 [43,45,48]. In recent years, some studies have described a role for canonically activating histone modifications at a subset of repetitive elements [49–54]. Indeed, this body of work has suggested that some long interspersed nuclear elements (LINEs) can bear marks such as the transcriptionally activating H3K4me3 modification, particularly early in development [51–53]. Similarly, other work has suggested that a class of mammalian-wide interspersed repeats (MIRs) may be transcriptionally active and play a role in enhancer regulation [54]. Much of this work, however, has relied upon uncalibrated ChIP with antibodies of uncertain specificity, both of which can result in data distortion and biologically incorrect conclusions [30]. Further, the ChIP-seq and RNA-seq studies have used a variety of different methods of aligning and filtering for reads to reach their conclusions, none of which used a method similar to our Bayesian SmartMap analysis, which may further affect the interpretations of the experiments. As such, we sought to use our calibrated and highly specific ICeChIP-seq datasets in conjunction with SmartMap to gain new insights into the epigenetic landscape of repetitive elements and to examine the degree to which uniread analysis yields an incomplete view of the data.

To accomplish this, we examined the histone modification landscape at the promoters of all repetitive elements, LINEs, short interspersed nuclear elements (SINEs), and Simple Repeats. K-means clustering analysis on all repetitive elements revealed four classes of promoters, each with a different histone modification profile: Cluster 1, enriched for H3K27me3 and H3K9me3; Cluster 2, enriched for H3K4me1 and H3K4me2; Cluster 3, enriched for H3K4me2 and H3K4me3; and Cluster 4, which is relatively depleted of histone modifications (Fig 6A). These clusters are roughly reminiscent of the functional classifications of the ENCODE hidden Markov model, where Clusters 1, 2, and 3 correspond to silenced promoters, enhancers, and active promoters, respectively [55]. Interestingly, in all but Cluster 4, a greater proportion of nucleosomes is enriched with H3K27me3 than H3K9me3, despite the previous emphasis on the latter in repetitive element silencing [9,43,45–47], emphasizing the importance of calibration in ChIP-seq studies for comparing different modifications [29,30]. Similar histone modification profiles are seen for the LINEs (S13A Fig), SINEs (S13B Fig), and Simple Repeats (S13C Fig). Across all these classes, Cluster 3 had the highest ATAC-seq signal (S13D–S13G Fig), consistent with the presence of histone modifications associated with transcription and accessible chromatin [10,13,56].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Assessment of histone modifications at promoters of repetitive DNA elements.

(A) Mean histone modification densities (HMDs) about promoters for classes of all repetitive elements, as defined by k-means clustering. Corresponding analyses of LINE, SINE, and Simple Repeat elements in S13 Fig. (B) Heatmap of repeat promoters with newly measurable HMD in SmartMap analysis, sorted on first principal component of repetitive elements. (C) Proportion of each cluster comprised by each repeat class or family for all repeats (left), LINE elements (center), and SINE elements (right). All significance tests performed as post-hoc Bonferroni-corrected pairwise 2x2 chi-square tests. (D) Quantile boxplots of average normalized RNA-seq read depth across LINE elements for each LINE cluster. Solid line with marker represents 90^th percentile; dashed line with marker represents 95^th percentile. Significance test shows difference in median by Bonferroni-corrected pairwise Mood’s median tests. Significance markers: *p<0.01, **p<10⁻⁵, ***p<10⁻¹⁰.

https://doi.org/10.1371/journal.pcbi.1008926.g006

Importantly, SmartMap analysis enabled us to more accurately measure HMD and assign clusters than did unriead analysis. Overall, there were 142,392 promoters with nonzero HMD in the SmartMap analysis that displayed no measurable HMD within 200bp of the TSS across all five histone modifications in the uniread dataset; similarly little HMD was detected within 1kb of the same in the uniread dataset (S14 Fig). This increase in HMD was substantial; most such sites had meaningful levels of histone modifications (Fig 6B). A small subset was primarily H3K4me2/me3 predominant; a larger subset had high levels of H3K4me1/me2, and the remainder were primarily characterized by H3K27me3 and H3K9me3 (Fig 6B). These represent promoters that would have been misclassified as histone-modification-depleted under uniread analysis; it is only through proper allocation of multireads that we can measure their HMDs and assign them to the appropriate cluster of repeat elements.

The distribution of repetitive elements across these clusters revealed interesting patterns. The distribution of the repeat classes or families across the clusters are presented in Table 5 for all repeats, Table 6 for LINEs, and Table 7 for SINEs, and summarized in Fig 6C. Notably, amongst SINEs, MIRs were enriched in Cluster 3 (Fig 6C), consistent with previous descriptions of a class of transcriptionally active MIRs [54]. In addition, Cluster 3 was enriched for Simple Repeats across all repeat promoters, consistent with descriptions of Simple Repeats in and around protein coding genes in the literature [57]. Interestingly, Cluster 3 was enriched for the L2 subtype of LINEs, despite previous work primarily focusing on the role of H3K4me3 at L1 elements [52], representing a novel prediction of transcriptional activity of this family. To this end, using SmartMap analysis of the RNA-seq data, we found that the Cluster 3 LINEs had greater transcriptional activity than did the other clusters (Fig 6D), confirming the transcriptional activity suggested by the presence of H3K4me3. Collectively, these data demonstrate the risk in only focusing on unireads–namely, the risk of missing important classes of genomic features–and highlights the role of multiread analysis of both DNA and RNA in driving new biological discovery.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 5. Clustering of Repetitive Elements.

https://doi.org/10.1371/journal.pcbi.1008926.t005

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 6. Clustering of LINE Elements.

https://doi.org/10.1371/journal.pcbi.1008926.t006

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 7. Clustering of SINE Elements.

https://doi.org/10.1371/journal.pcbi.1008926.t007

Sequencing data sources

Computing average value of BEDGRAPH at target loci

Because the BEDTools map software does not compute base-pair-wise averages of BEDGRAPH signals, the following procedure was used to compute read depth at a list of target loci. Overlapping loci were merged using BEDTools merge, and the resultant list of loci were partitioned into 1bp windows using BEDTools makewindows. The BEDGRAPH was then mapped onto the windows using BEDTools map, and the mapped windows were then mapped with the mean function onto the original list of target loci using BEDTools map.

Mappability estimation and binning

The mappability of a list of loci was computed by computing the average value of the UMAP50 bedgraph for the relevant genome at those loci using the method described above. To compute the number of regions per UMAP50 score, the loci were binned by average UMAP50 score in bins of width 0.01. The number of loci at each bin were then computed to determine the approximate distribution of UMAP50 score across the selected loci.

MACS2 peak calling

Peak calling was conducted on the simulated datasets using MACS2 [64] with the bdgpeakcall function with the relevant BEDGRAPH file and default settings.

Alignment and read filtering and processing

FASTQ files for ICeChIP-seq or the simulated dataset were aligned using Bowtie2 [40] due to its common usage in the field and due to its ability to report alignment scores for each mate for each alignment reported as opposed to for just the best alignment. Bowtie2 alignment was run on the paired-end sequencing samples with the following settings: end-to-end alignment, very-fast preset settings, no discordant alignments, no mixed alignments, report up to 51 alignments, insert size minimum 100bp, insert size maximum 250bp. In the case of the analysis with up to 101 alignments (the k101 dataset), the above settings were used with up to 101 alignments reported per read. The genomes used for alignment were as follows: AR7, mm10 with ICeChIP barcodes series 1; AR8, dm3 with ICeChIP barcodes series 1; AR9, mm10 with ICeChIP barcodes series 2; AR16 and AR17, hg38 with ICeChIP barcodes series 3; simulated datasets, hg38.

FASTQ files for RNA-seq were aligned using Hisat2 [65] for the same reasons as the choice to use Bowtie2. Hisat2 alignment was run on the paired-end sequencing samples with the following settings: no discordant alignments, no mixed alignments, report up to 51 alignments. The genome used for alignment was hg38 with the ENCODE ERCC standards.

FASTQ files for ATAC-seq were aligned using Bowtie2 [40] on the paired-end sequencing samples with the following settings: no discordant alignments, no mixed alignments, report up to 51 alignments, insert size maximum 2000bp.

Alignments were then filtered to select for reads that are paired, mapped in a proper pair, and mate on the reverse strand, corresponding to SAM flags of 99, 163, 355, and 419. For non-strand-specific applications, the selected SAM file records were then extracted into a file containing the following fields: chromosome, start position, stop position, read ID, read alignment score (field labeled “AS:i:”), mate alignment score (field labeled “YS:i:”). For strand-specific applications, the selected SAM file records were extracted into a file containing the following fields: chromosome, start position, stop position, read ID, strand, read alignment score (field labeled “AS:i:”), mate alignment score (field labeled “YS:i:”). The reads were then split into separate BED files based on the number of alignments per read. For downstream uniread analyses, only the reads with 1 alignment were used; for downstream SmartMap analyses, reads with 1–50 alignments were used except for the k101 dataset, in which case reads with 1–100 alignments were used.

The file with 51 alignments per read (or that with 101 alignments per read for the k101 analysis) was not used for downstream analyses to prevent confounding by reads with fewer reported than possible reads.

Uniread and SmartMap analysis of genome coverage

For the uniread analysis, our SmartMap software was used with only the file containing reads with only 1 alignment per read. For the SmartMap analysis, our SmartMap software was used with the files containing reads with fewer than 51 alignments per read.

The SmartMap software uses a set of dual Binary Indexed Trees to store map counts and weights across the genome and uses an iterative Bayesian reweighting algorithm to assign weights to each of the different alignments. These steps are outlined below. Unless otherwise specified, all analyses are conducted with 1 iteration in scored mode. For the strand-specific analyses, there is a separate set of dual Binary Indexed Trees for each strand.

Storage of map counts in the genome.

To facilitate efficient summation and updating of map counts and weights across the genome, each chromosome is stored as a pair of Binary Indexed Trees (BIT), also known as Fenwick Trees. The BIT is a data structure that is efficient for computing prefix sums of an ordered dataset from the beginning of the dataset to the given index. Because we used a 1-based coordinate system for the genome, the datasets to which we refer as being represented by a BIT should be assumed to be 1-based datasets unless otherwise specified.

For a dataset of length L, the BIT is represented by L+1 nodes, which are stored in an array. To increment a dataset represented by a BIT T at index i by the value v, the following algorithm is used. Let T[i] represent the ith node of T. Let lsb(i) represent the lowest significant bit in the binary representation of i. Then: (1) (2)

If the new value of i ≤ L+1, Eqs 1 and 2 are iterated until i > L+1. For brevity, we will refer to this operation to increment the BIT T representing the 1-based dataset by v in the value i as BITUpdate(T, i, v).

To compute the prefix sum of the dataset at index i (i.e. the sum of all values with indices [1, i] of a 1-based dataset), the following algorithm is used, using the above definitions of T[i] and lsb(i). Let the prefix sum be represented by sum, where sum = 0 at the beginning of the algorithm. Then: (3) (4)

If the new value of i > 0, Eqs 3–4 are iterated until i ≤ 0. For brevity, we will refer to this operation to obtain the prefix sum of the T at value i as BITSum(T, i).

To understand how we here use BITs to efficiently store values across the genome and efficiently sum the values across loci, consider the following.

Consider a dataset C represented by BITs T₁ and T₂. If the values of C for indices in range [l, r) are incremented by v, then let the resulting dataset be represented by C’. Let ΔPointSum(C, i) = PointSum(C’, i)–PointSum(C, i).

The prefix sum of the resulting dataset C’ at index i, represented by PointSum(C’, i), is changed in one of three ways:

Case 1: i < l. In this case, the increment on range [l, r) will not change PointSum(C’, i). As such: (5)

Case 2: l ≤ i < r. In this case:

However, per Eq 5, ΔPointSum(C, l-1) = 0. As such: (6)

Case 3: i ≥ r. In this case, the increment on range [l, r) will not change the values of C past index r-1. Accordingly:

However, per Eq 6, ΔPointSum(C, r) = v*r–v*(l-1). As such (7)

These three cases will provide the basis for our use of BITs to store and efficiently sum values across the genome. Each chromosome C in the genome is stored as a pair of BITs, to which we shall here refer as T₁ and T₂. Let L represent the length of the chromosome. Accordingly, T₁ and T₂ have L+1 nodes.

To increment the value associated with the base pairs in the range [l, r) by an increment value v, the following procedure is used.

(8)

The value associated with base pair i is then BITSum(T₁, i).

To obtain the prefix sum of the chromosome dataset C at base pair index i, represented by PointSum(C, i), the following equation is used.

(9)

The sum of the values associated with each base pair in the range [l, r) = [l, r-1] on chromosome C, represented by LocusSum(C, l, r), can thus be described by: (10)

This dual-BIT data structure allows for efficient handling of data with respect to time complexity. The BITUpdate and BITSum steps occur with time complexity O(log L), and the updates to ranges (Eq 8) and range summations (Eq 10) use four BITUpdates and four BITSums, respectively. As such, both range updates and range queries occur with time complexity O(log L).

Histone modification density and specificity computation

Because the ICeChIP-seq datasets have internal nucleosome standards bearing the targeted nucleosome modifications with uniquely identifying DNA “barcodes”, we were able to calibrate our ChIP-seq results to yield the histone modification density (HMD), or the proportion of nucleosomes bearing the modification of interest. HMD for each dataset was computed as follows.

The average value of the BEDGRAPH for each of the calibrant barcodes was computed as above, and these values were grouped by the nucleosome modification associated with the barcode and summed, as previously described [29,30] for both the IP and the Input datasets. The ratio of the summed values for the targeted modification in IP over the same in Input was designated as the target enrichment E_t.

The HMD at each genomic locus was then computed as follows, where IP and Input represent the value of the IP and the Input at that genomic locus: (18)

To generate genome-wide HMD BEDGRAPH files, the IP and corresponding Input genome coverage BEDGRAPH files outputted by the SmartMap software were merged with BEDTools unionbedg, and the HMD computation in (18) was used to compute HMD. Any region with an Input value of zero was set to an HMD of zero.

To compute the specificity, for those ICeChIP-seq datasets with calibrants bearing more than one modification with uniquely identifying DNA barcodes (AR9, AR16, and AR17), the enrichment of every species (E_i) was computed analogously to the E_t, and the specificity (as percent of target enrichment) was computed as: (19)

Alignment overlap analysis

To assess for overlap of alignments, we conducted SmartMap on the simulated dataset with the read weight output setting activated. Using bedtools intersect, we then identified alignments that intersected with the true read origin in the Gold Standard dataset. From this, by weight, we were able to compute three metrics. First, we computed the number of alignments by weight that were present in the intersected dataset as a proportion of the total number of alignments by weight. Second, we computed the alignment weighted overlap proportion score, a measure of the proportion of a read’s overall weight that overlaps with a given true origin of the read due to a given alignment. This is computed as the product of the weight of the alignment with the geometric mean of the proportion of overlap between the true read locus and the alignment locus. Finally, we computed the unweighted overlap proportion score, which is computed as the geometric mean of the proportion of overlap between the true read locus and the alignment locus.

Repetitive element analysis

Repetitive elements for hg38 were obtained from the HOMER list of repeats [66]. The promoter was defined as the most upstream portion of the annotated repeat. This dataset was used for analyzing all repeats; for analyzing LINE elements, SINE elements, or Simple Repeats, the corresponding subset of the repeats was used.

The HMD profiles in Figs 6A and S10 were generated by computing the average HMD (from SmartMap analysis of AR16 and AR17) in 50bp windows from -1000bp to +1000bp relative to the promoter, with HMDs above 100% being set to 100% (because an HMD above 100% is definitionally impossible), and corresponding windows were averaged together to yield the average HMD profile for each set of elements.

To conduct clustering, first, the average HMD of the region -100bp to +100bp relative to each promoter in the relevant dataset was computed using the SmartMap analysis of AR16 and AR17, with HMDs above 100% being set to 100%. The data was then transformed to orthonormal basis by principal component analysis in R with scaling and centering. The resultant coordinate matrix used for k-means clustering, starting with 2 clusters and increasing the number of clusters until the decrease in total within-cluster sum of squares became markedly diminished; for each dataset (all repeats, LINE elements, SINE elements, and Simple Repeats), this occurred with 5 or more clusters and, accordingly, 4 clusters were used for each dataset.

RNA-seq analysis was conducted as follows. The average value of the RNA-seq SmartMap BEDGRAPH datasets were found across each LINE element. These values were then normalized to the SmartMap read count for each replicate (as average reads per million reads analyzed) and averaged to yield the average normalized read depth for each LINE element. These were then grouped by cluster and used to generate the quantile boxplots.

Genome browser visualization

Genome browser visualization was conducted using Integrative Genomics Viewer (IGV) [67].

Comparison to other methods

Statistical analyses

Statistical analysis for Fig 6C was conducted first with chi-square analysis on full contingency table and with post-hoc tests on collapse contingency tables as follows. For each of the datasets in Fig 6C, chi-square test for goodness-of-fit was conducted on the corresponding contingency table presented in Tables 4–6. The p-value for each of these tests was p<2.2x10^-16 and, accordingly, post-hoc tests were conducted. The post-hoc tests consisted of collapsing each contingency table into a set of 2x2 contingency tables with the cluster of interest and family/type of interest compared to all other clusters and/or all other families within the contingency table. Chi-square goodness-of-fit tests were then conducted on each of these 2x2 contingency tables, and the p-values were Bonferroni corrected to adjust for the number of tests. These adjusted p-values for each 2x2 contingency table test were used to label the graphs in Fig 6C as follows: *p<0.01, **p<10⁻⁵, ***p<10⁻¹⁰.

Statistical analysis on Fig 6D was conducted to compare median average normalized RNA-seq depth by cluster. Because the difference between cluster 3 and all of the other clusters appeared to be the most biologically meaningful, only pairwise comparisons were conducted between cluster 3 and the other clusters to limit the number of statistical comparisons and, accordingly, the degree of Bonferroni correction needed. Mood’s median tests were solely conducted as pairwise comparisons between cluster 3 and each of the other clusters with Bonferroni correction to p values with n = 3 for Bonferroni correction. The adjusted p-values for each of these comparisons was p<10⁻¹⁰ and was marked appropriately on the graph.