Genomic Characterization of Large Heterochromatic Gaps in the Human Genome Assembly
[['Long arrays of near-identical tandem repeats , termed satellite DNAs , compose the predominant sequence feature within constitutive heterochromatin in complex genomes [1] .', 'Genomic regions enriched in satellite DNA are often found in or near centromeres [1] , [2] , and , in some species , they can account for greater than half of the entire genetic material [3] .', 'Yet despite the abundance and potential biological functions of satellite-rich regions , reference assemblies typically exclude them due to the inability of standard assembly algorithms to uniquely align and assemble near-identical tandem repeats .', 'Efforts to characterize satellite DNA arrays across diverse taxa have proven useful for surveying highly represented copies of tandem repeats , typically reporting a single consensus reference or a small number of representative satellite sequences within a given repeat family [4]–[6] .', 'Only in a few instances have satellite DNAs been studied using a comprehensive inventory of sequences that include all instances of a particular repeat family within the genome , a necessary step for understanding the sequence composition and evolutionary processes that govern sequences within satellite-rich regions of complex genomes [7] .', 'Additionally , attempts to detect and characterize novel tandem repeats are often dependent on long read lengths and on the absence of interspersed stretches of exact sequence identity within each tandem repeat [8] .', 'This presents a considerable problem for characterizing satellite DNA families that are defined by an irregular repeat unit length composed of complex arrangements of short repeats .', 'As a result , this type of complex satellite DNA remains largely uncharacterized , even within the well-characterized human genome .', 'Further , unlike those methods employed to characterize satellite families with well-defined tandem repeat lengths ( e . g . [5] , [6] , [9] ) , limited sequence-based tools exist to explore the nature of short , irregular repeat sequences across diverse genomic datasets .', 'Human Satellites 2 and 3 ( here , referred to as HSat2 , 3 ) define large , heterochromatic blocks adjacent to many centromeric regions in the human genome [10]–[14] .', 'Due to their complex arrangement of divergent 5-bp repeats and the inconsistent lengths of their tandem repeat units , these sequence families remain poorly understood in the human genome , with only a limited sequence characterization of each family collected from a small number of experimentally derived-clones ( e . g . [14]–[16] ) and from surveys of small stretches of HSat2 , 3 sequences included in the human reference assembly [17] , [18] .', 'HSat2 , 3 likely originated from an ancestral pentameric simple sequence ( CATTC ) [14]; however , notwithstanding this commonality , hybridization-based studies [11]–[14] , [16] , [19] , [20] have shown that these two satellite families largely occur in spatially distinct arrays on different chromosomes .', 'Further , characterization of clones derived from different HSat2 and HSat3 sites showed that each array typically has a distinct array length , repeat unit length , and sequence composition [15] , [16] , [21]–[25] .', 'While some HSat2 , 3 arrays are nearly identical to arrays on other chromosomes ( e . g . , on the acrocentric chromosomes [15] , [16] , [21] , [22] , [24] , [25] ) , others can be readily distinguished by hybridization-based methods [10] , [16] , [19] .', 'For example , previous studies determined that the large , multi-megabase array of HSat2 adjacent to the centromere on chromosome 1 is largely defined by a chromosome-specific 1 . 77 kb repeat unit [16] .', 'Similarly , a male-specific 3 . 6 kb HSat3 repeat unit comprises the DYZ1 satellite array on the long arm of the Y chromosome [26]–[28] .', 'These limited experimental results suggest the possibility of additional definable subfamilies of HSat2 , 3 within the human genome , although no study to date has created a comprehensive inventory of these subfamilies using whole-genome sequence datasets .', 'The length and sequence composition of satellite DNA arrays are known to be variable among individuals , providing a potentially rich source of information to study human variation .', 'The rapid expansion and contraction of satellite arrays by mechanisms such as unequal crossover or sequence conversion [29] , [30] is expected to alter the frequency and occurrence of sequence variants between individuals .', 'Consistent with this , previous cytogenetic studies have demonstrated that the large , HSat2 , 3-rich heterochromatic bands on chromosomes 1 , 9 , 16 , and Y can vary significantly in size among different individuals [31]–[34] .', 'Further , because the DYZ1 array can occupy more than half the length of the Y chromosome [26] , variations in the size of this array can dramatically alter the total size of the Y chromosome , and limited studies have suggested that Y-chromosome size may vary across human populations [34] , [35] .', 'Although HSat2 , 3-rich regions of the genome are defined primarily by long , uninterrupted arrays of satellite sequences , several small-scale studies have provided evidence for ‘islands’ of non-satellite sequences that , due largely to their association with satellite-rich regions , are not included in chromosomal reference assemblies [36] , [37] .', 'Therefore genome-wide studies of satellite-rich domains may uncover both repetitive and non-repetitive DNA currently missing from reference assemblies .', 'The absence of satellite-rich regions from reference assemblies prevents detailed study of their potential function ( s ) .', 'In at least some organisms , it is becoming increasingly clear that sequences within these regions can contribute to important cellular functions , such as meiotic chromosome pairing and maintenance of centromere function [38] , [39] .', 'Moreover , epigenetic and transcriptional derepression of HSat2 , 3 sequences is a common feature of cancer cells and stressed cells , coupling the chromatin state of these sequences with genome stability [40]–[45] .', 'Satellite DNAs are also overrepresented in extrachromosomal circular DNAs [46] and in cytoplasmic membrane-associated DNA , a poorly understood class of DNA fragments that have been shown to exit the nucleus and associate with the plasma membrane , where they can be transcribed [47] .', 'Without any sort of reference sequence for satellite-rich regions , however , it is difficult to study their expression and potential regulation at finer scales .', 'Additionally , it is unclear if variation between individuals – in both the satellite and non-satellite components of these regions – has any influence on the function and epigenetic behavior of these regions .', 'A more complete genomic census of satellite DNAs would also aid studies of repeat-associated chromosomal rearrangements , as HSat2 , 3 sequences are commonly enriched at the breakpoints of Robertsonian translocations [48]–[51] .', 'To address these current limitations , we present here an algorithm to predict HSat2 , 3 satellite subfamilies and their chromosomal localizations within the human genome , principled on sequence features typical of satellite DNA arrays .', 'For each satellite subfamily , we generate a database of whole-genome shotgun ( WGS ) reads capturing the full spread of satellite and non-satellite sequences associated with that subfamily , useful as a nonlinear reference sequence for further studies of these regions .', 'We then apply this reference database to short-read WGS sequence data from several hundred individuals to examine HSat3 array size polymorphisms within the human population .', 'This work presents an initial genomic definition of classical HSat2 , 3 sequences , which underlie the largest heterochromatin-assigned gaps in the human reference assembly , and it offers a framework to study similar satellite-rich regions characteristic of most other complex genomes .'], ['To build a reference library of all the HSat2 , 3 sequences present in the human genome , we first applied RepeatMasker [52] to all WGS reads from the genome of a single male donor sequenced at ∼7 . 8x coverage with paired Sanger sequencing reads [53] , yielding a comprehensive database of 396 , 228 HSat2 , 3-containing reads ( Dataset S1 ) .', 'We verified the comprehensiveness of this database using an independent method that tests for enrichment of the defining CATTC 5-mer on each read ( see Methods ) .', 'Assuming uniform genomic coverage , this database corresponds to roughly 1 . 5% ( 48 Mb ) of the human genome .', 'To identify shared sequence homology among the HSat2 , 3-containing reads without direct alignment , we first represented each read as a k-mer frequency vector as follows:', '1 ) slide a k bp window along each strand-oriented read one base at a time ,', '2 ) count the number of occurrences of each distinct k-mer , and', '3 ) divide by the length of the read .', 'We chose a k of 5 to reduce memory complexity and to reflect the pentameric nature of the ancestral HSat2 , 3 repeat sequence ( CATTC ) .', 'Additionally , we demonstrate that a modest increase ( where k\u200a=\u200a4–7 ) provided limited , if any , improvement in our final subfamily classifications ( Figure S1 ) .', 'As there are 1 , 024 ( 45 ) possible 5-mers in DNA , we effectively mapped each read onto a vector in a feature space with 1 , 024 dimensions ( Figure 1 ) .', 'In this feature space , reads with similar 5-mer frequencies will tend to be closer together , allowing us to estimate pairwise sequence distance simply as the Euclidean distance between two read vectors ( as in [7] ) .', 'Furthermore , reads mapped onto this feature space can be easily projected in lower dimensions using techniques such as principal component analysis ( PCA ) , allowing one to visualize any predominant clusters .', 'Next , we performed pairwise comparisons of k-mer frequency vectors to identify HSat2 , 3 reads that are near-identical in sequence composition and thus likely to originate from the same satellite subfamily .', 'Further , as satellite arrays consist primarily of uninterrupted tandem repeats that span regions commonly multiple megabases in size ( i . e . greater than the insert length between paired reads ) , one would expect reads originating from a given satellite array to have mate pairs that predominantly contain similar repeat sequences from the same array .', 'That is , read clusters that define each HSat2 , 3 satellite subfamily are expected to have high sequence similarity as well as high ‘self-mate-pair’ frequencies .', 'To generate unsupervised read clusters using these criteria , we developed a modified k-nearest-neighbor spectral clustering algorithm ( see Methods ) .', 'To identify clusters based on read similarity , we determined the 50 nearest neighbors for each read vector in the 5-mer frequency feature space , resulting in a weighted undirected graph with vertices corresponding to read vectors , edges corresponding to nearest-neighbor relationships in either direction , and edge weights inversely proportional to sequence similarity quantified as the respective normalized Euclidean distance .', 'Importantly , we then overlaid maximally weighted edges corresponding to mate-pair relationships between reads , yielding a graph representing both sequence and mate-pair relationships within the read dataset .', 'Next , we performed normalized spectral clustering to cut this graph into two large groups while severing the least total edge weight possible ( Figure', '1 ) [54] .', 'After this binary split , we then calculated the self-mate-pair frequency of each resulting cluster , which is the proportion of all mate pair edges from each cluster that have both read vectors in the same cluster .', 'We accepted binary divisions resulting in self-mate-pair frequencies >0 . 8 ( shown by paired read simulation to identify satellite arrays greater than 100 kb , as described in Text S1 ) .', 'The resulting daughter clusters were evaluated in the same manner , extending recursive binary divisions until a daughter cluster failed to satisfy the self-mate-pair frequency threshold .', 'After the initial binary division , we subdivided the read database into HSat2 and HSat3 satellite clusters ( as shown in the uppermost panel of Figure 2 ) , with self-mate-pair frequencies of 0 . 997 and 0 . 998 ( HSat2 and HSat3 , respectively ) , consistent with previous observations that these sequences , although related by an ancestral pentameric repeat sequence , can be distinguished by sequence and spatial organization in the human genome [13] .', 'After recursive binary divisions , this method yielded 14 distinct HSat2 , 3 subfamilies ( three corresponding to HSat2 and 11 to HSat3 ) , and we computed coverage-based estimates of the abundance of each subfamily within the genome ( Figure 2 , Table 1 ) .', 'Comparisons with previously experimentally verified HSat2 , 3 clone sequences [15] , [16] , [21]–[25] , [55]–[57] showed high concordance with our predicted HSat2 , 3 subfamily datasets , indicating that the reported subfamilies comprise genomically distinct subsets , not simply different parts of the same satellite arrays ( Table S1 ) .', 'To further evaluate the read clusters corresponding to these subfamilies , we computed the self-mate-pair frequencies of all pairwise combinations of the final fourteen clusters ( Table S2 ) .', 'Similar to self-mate-pair frequencies in the initial HSat2 and HSat3 clusters , the majority ( 10/14 ) of the subfamily predictions are defined by self-mate-pair frequencies >90% when compared with all other clusters , providing strong support for these cluster definitions .', 'The remaining four subfamilies ( HSat2A1 , 3A2 , 3B1 , and 3B2 ) show somewhat lower self-mate-pair frequencies ( 73%–83% ) , likely owing to the fact that these subfamilies share subregions with high homology to other subfamilies ( Figure S2 and Table S3 ) .', 'The resulting HSat2 , 3 read clusters each represent a deep genomic inventory of the sequences that comprise a physically distinct HSat2 , 3 subfamily , without use of a consensus sequence .', 'These satellite reference databases are expected to contain essentially all satellite sequence variants as well as representative samples of any non-satellite sequences associated with each subfamily .', 'To study the sequence features specific to each HSat2 , 3 subfamily , we analyzed the nature of individual cluster topologies and performed direct sequence comparisons within and between each corresponding satellite reference dataset .', 'Five subfamilies ( HSat2A1 , 2A2 , 2B , 3A5 , and 3A6 ) are observed to have ‘ring-like’ topologies visible in their PCA projections .', 'We hypothesized that reads within these ring-like projections may collectively represent a predominant repeat unit whose length is longer than a single WGS read and whose subregions are sufficiently diverged to be distinguished in the feature space .', 'As depicted in Figure 3 , if a highly represented tandem repeat unit , shown subdivided into parts a , b , c , is prevalent , then groups of reads are expected to overlap a-b ( red to orange ) , b-c ( orange to yellow ) , and then c-a ( yellow to red ) , resulting in a circular transition from the end of one repeat unit to the start of the next , found immediately upstream , conceptually similar to the notion of a multimeric , higher-order repeat , defined previously for other satellite DNAs [58] .', 'To test this hypothesis , we simulated a read database from a concatenated array of a published 1 . 77 kb HSat2 repeat unit sequence corresponding to HSat2A2 [15] , [59] .', 'By plotting our simulated data relative to our full sequence database for subfamily HSat2A2 ( Figure 3 ) , we show that the general topology corresponds to the presence and prevalence of a tandemly repeating unit .', 'In a similar manner , we demonstrate that the previously published repeat unit for DYZ1 [28] was useful in describing the topology of HSat3A6 ( Figure S3a ) .', 'Notably , many subfamily cluster predictions do not have a similar ring-like topology , perhaps due to short or inconsistent repeat units , as demonstrated for a shorter repeat unit in the HSat3A4 subfamily ( 1366 bp , pTRS-47 , Figure S3b ) .', 'Sequence relatedness among individual subfamilies is observed within the hierarchical tree describing their subsequent subgraph divisions ( as shown in Figure 2 ) .', 'For example , HSat2A1 and HSat2A2 share a parent subgraph and are closely related on the sequence level , suggesting recent evolutionary sequence divergence .', 'Notably , we show that HSat3 can be divided into two large and very well-separated families ( designated as HSat3A , B ) , each of which is composed of multiple subfamilies .', 'To further estimate the level of similarity between subfamilies , we computed the fraction of overlapping 24-mers between each pair of subfamily read databases .', 'In line with previous observations , HSat2 and HSat3 subfamilies are generally highly distinguishable at this sequence level ( i . e . <1% of overlap ) , with the most prevalent overlap ( up to 57% ) observed between closely related subfamilies ( Table S3 ) .', 'Inversely , this multi-way comparison across all HSat2 , 3 clusters is useful in predicting those 24-mers that are highly enriched in each subfamily relative to all others ( Table 1; Table S4; Dataset S2 ) , making them useful for classifying short sequencing reads among subfamilies and for designing subfamily-specific oligonucleotides for experimental applications .', 'To estimate the distribution of each HSat2 , 3 subfamily across the genome , we utilized available HSat2 , 3 sequences present in the current chromosome reference assemblies ( GRCh37 ) and in Genbank , as well as whole chromosome shotgun ( WCS ) sequence libraries ( representing data from 19 flow-sorted chromosomes [60] , [61] , Tables S5 and S6 ) .', "To compute the coverage of each HSat2 , 3 subfamily in each chromosome's WCS dataset , we first computed 5-mer frequency vectors for 186 , 961 HSat2 , 3-containing WCS reads , as described earlier , and assigned each one to the same subfamily as its nearest read within the initial reference database .", 'We accounted for experimental contamination from non-target chromosomes using a Bayesian model with an uninformative prior ( Table S7 as discussed in Methods ) .', 'Consistent with the chromosomal assignments of published HSat2 , 3-containing clones , we observe chromosome 1 enrichment for HSat2A2 [15] , chromosome Y enrichment for HSat3A6 [32] , and acrocentric chromosome enrichment for HSat3A4 [55] , [57] ( Figure 4; Table 1 ) .', 'Further , we identified HSat3B5 as the predominant sequence feature of the large , pericentromeric heterochromatic band on chromosome 9 , which was previously known to contain HSat3 but lacked any complete clone sequence [15] .', 'We also describe a previously uncharacterized subfamily specific to the Y chromosome , HSat3A3 .', 'In line with published high-stringency , hybridization-based mapping experiments , we find that many of our predicted HSat2 , 3 subfamilies are distributed across more than one chromosome [13] , suggesting the possibility of inter-chromosomal sequence homogenization .', 'In addition to WCS assignments , we studied those HSat2 , 3 sequences currently present in the human reference assembly ( GRCh37/hg19 ) [62] , including those found directly adjacent to heterochromatin gaps ( 983 kb total ) , as well as a much smaller number interspersed within chromosome arm assemblies ( 31 kb total ) .', 'We assigned these assembly sequences to subfamilies by searching for perfect matches to our subfamily-specific HSat2 , 3 24-mer libraries ( Table S8 ) .', 'Notably , we found that individual subfamilies are observed in small amounts in several locations in the human genome assembly , with the majority observed directly adjacent to ( within 500 kb of ) the centromeric or heterochromatin-assigned gap [18] , [63] .', 'While most HSat2 , 3 annotations on the assembly are concordant with our WCS-based chromosomal localization predictions , our WCS coverage data lack power to rule out the existence of small ( <50 kb ) interspersed arrays of HSat2 , 3 on some chromosomes .', 'Some genomic regions enriched with HSat2 , 3 have been shown previously to contain islands of non-satellite sequences missing from the reference assembly [36] , [37] .', 'To detect and characterize the presence of such sequences more systematically , we first obtained all unmapped scaffolds from the HuRef assembly that contain at least one uniquely mapping 24-mer ( 18 , 780 scaffolds , totaling 43 Mb , including 22 Mb of non-RepeatMasked sequence; see Methods ) .', 'To determine the chromosomal origin of each scaffold , we used a Bayesian model taking as input sequence coverage from flow-sorted chromosome ( WCS ) datasets and an uninformative prior ( see Methods , Dataset S3 ) .', 'Using this WCS coverage model , we were able to predict the localization of 9 , 950 HuRef scaffolds to at least one chromosome , totaling 13 . 7 Mb of non-repetitive sequence .', 'We also localized 555 unmapped scaffolds found in the “decoy sequences” utilized by the 1000 Genomes Project Consortium [64]–[66] .', 'We compared our localizations with available admixture mapping data from two recent studies and observed our chromosomal assignments to be highly concordant ( 88 . 89%; 313/352 ) ( Dataset S3 ) .', 'Additionally , we evaluated unmapped scaffolds with reported chromosomal localizations in GenBank ( including GRCh37 and GRCh38 ) and found roughly ∼80% ( 247/310 ) of our chromosomal predictions to be concordant .', 'Furthermore , 33/35 unmapped scaffolds with available fluorescence in situ hybridization ( FISH ) mapping data had concordant chromosomal assignments [53] , [67] ( Table S10 ) .', 'We observed that the majority of discordant scaffold placements could be attributed to acrocentric chromosome localization ( where chromosomes 13 , 14 , 21 , and 22 are known to share considerable sequence identity ) and/or to low-confidence chromosome localization ( e . g . scaffolds that had lower-confidence admixture mapping assignments or lacked whole-chromosome shotgun sequencing data ) .', 'We identified 17 large scaffolds ( totaling 1 . 3 Mb of non-repetitive sequence ) with three or more mate-pair connections to the same HSat2 , 3 subfamily , providing evidence for physical proximity to regions enriched in HSat2 , 3 .', 'We evaluated the chromosomal assignments of eleven scaffolds with available mapping results from admixture studies [66] , [68] and found them all to be concordant in their chromosomal localizations .', 'Additionally , we report four out of the remaining six scaffolds to have chromosome assignments concordant with available FISH mapping data [53] ( Table S9 ) .', 'As an example , a single HuRef unmapped scaffold containing roughly 60 kb of segmentally duplicated sequence paralogous to multiple chromosomes was shown to define an inter-satellite sequence region on chromosome 1 by WCS enrichment , PCR-based mapping ( Figure S4 ) , and paired-read mapping to HSat2 and adjacent alpha satellite ( D1Z5 ) ( Figure 5 ) .', 'Similar to previous studies [36] , [37] , [53] , we observe that the majority ( 86% ) of non-satellite sequences that are associated with HSat2 , 3 arrays represent low-copy , segmentally duplicated sequences , including eight annotated gene predictions ( Table S9 ) .', 'The introduction of inter-satellite regions is expected to improve the completeness of ongoing genomic studies aimed to understand the evolution and biological functions of these sequence paralogs .', 'Satellite arrays have been shown previously to vary in size within the human population [31]–[34] , [69]–[71] .', 'Therefore , to illustrate the utility of our HSat2 , 3 reference database for studying satellite array size variation using short-read sequencing data , we investigated the total size of the Y-chromosome heterochromatic locus classically known as DYZ1 ( corresponding to subfamily HSat3A6 ) using genomic sequence data from 396 male individuals from the 1000 Genomes Project [64] , [65] .', 'Predicted DYZ1 array sizes ( see Methods , Table 2 , Dataset S4 ) varied over an order of magnitude ( 7-98 Mb , with mean 24 Mb ) , consistent with previous experimental observations of Y-chromosome length variability [27] , [70] , [71] .', 'Because HSat3A6 is on the non-recombining portion of the Y chromosome , it is inherited as a single haplotype linked with previously characterized sequence polymorphisms that distinguish Y haplogroups ( reviewed in [72] ) .', 'Thus , one predicts that different haplogroups would contain different HSat3A6 array length distributions given their independent evolutionary histories since splitting .', 'Further , the variability within each haplogroup allows one to observe how rapidly satellite array lengths evolve on the time scales since Y haplogroup divergence .', 'As illustrated in Figure 6 , variable HSat3A6 array length distributions are observed both within and between haplogroups , consistent with a model of rapid satellite array length evolution .', 'To test whether array size distributions differ significantly within any haplogroup compared with all other haplogroups , we performed a two-sided , two-sample Wilcoxon rank-sum test for each pairwise combination of haplogroups , using a Holm correction for multiple testing .', 'Six pairwise combinations show highly significant location shifts ( p<0 . 001 ) .', 'For example , array sizes in the Asian-enriched O haplogroup are stochastically greater than array sizes in the African-enriched E haplogroup ( p\u200a=\u200a3 . 05×10-7 ) , with a median difference of 7 . 6 Mb .', 'Notably , the largest array sizes observed occurred within Finnish males ( in haplogroups I , N ) , but other samples with large arrays appear in haplogroups enriched in Asian and African populations ( E , O , R ) .'], ['Achieving more complete reference genome assemblies will require better characterization and representation of heterochromatic regions , which are typically composed of long arrays of satellite repeats along with embedded non-satellite sequences .', 'Here , we present a computational framework for studying both the satellite and non-satellite components of heterochromatic genomic regions , optimized for satellite families that are composed of a complex arrangement of simple repeats , a common feature of satellite families across diverse taxa .', 'Within this study , we have focused in particular on Human Satellites 2 and 3 , which are composed of diverged arrangements of 5-bp repeats and constitute roughly 1 . 5% of the human genome , typically occupying heterochromatic blocks adjacent to centromeric regions .', 'We used an alignment-free approach to detect and cluster highly similar HSat2 , 3 sequences among genome-wide collections of WGS paired reads , yielding 14 subfamilies distinct in sequence composition and physical location in the genome .', 'Further , by utilizing sequence libraries of isolated human chromosomes , we predict the genomic distribution of these subfamilies , presenting chromosome-assigned satellite reference databases capable of describing the sequence variation within multi-megabase HSat2 , 3 arrays .', 'Our computational framework is broadly applicable for clustering satellite DNA paired-read libraries , which can be initially compiled from WGS data according to their enrichment of simple repeat patterns ( in our case , ‘CATTC’ ) .', 'Further application of our strategy to less characterized , or novel satellite DNA families is possible by performing an initial survey of short repeat enrichment in a given genome using methods presented here and/or by employing existing software designed for short repeat discovery [4] , [8] to generate an initial read database for subfamily classification .', 'As an extension of these reference genomic databases , we developed a library of HSat2 , 3 subfamily-specific 24-mers useful for estimating the size of each subfamily across short-read , whole-genome datasets of individual genomes from diverse human populations .', 'Frequent non-homologous sequence exchange is expected to promote rapid expansion and contraction of repeat variants within each satellite DNA array .', 'Focusing on a single , haploid HSat3-rich region on the long arm of the Y chromosome , we provide satellite array length estimates , identifying a substantial difference in array length across 396 male individuals from distinct Y haplogroups .', 'Moreover , we demonstrate that this variation is observed within closely related members of the same haplogroup , illustrating the potential for recent sequence expansion and contraction within these sites .', 'This analysis and the methods described herein provide an opportunity to characterize not only satellite DNA but also non-satellite DNA sequences found in close proximity to each HSat2 , 3 subfamily .', 'In doing so , we have localized roughly 1 . 3 Mb of previously unmapped HSat2 , 3-associated sequences to particular chromosomes , providing a new potential source of assayable polymorphic markers to tag the heterochromatic gaps of the human genome for studies of recombination , population genetics , and disease association .', 'These scaffold sequences , along with our satellite reference database and our library of subfamily-specific 24-mers , also provide potential mapping targets to study the epigenetics and expression of these poorly characterized regions of the human genome .', 'This work thus serves as an initial reference database for HSat2 , 3-enriched regions within a single diploid genome , providing a foundation for ongoing genomic studies aimed at understanding their potential biological function and cellular regulation .'], ['We downloaded all ∼32M unassembled WGS sequencing reads generated by the HuRef genome project [53] and we processed them using RepeatMasker ( [52]; parameters: slow , cross-match ) .', 'Candidate HSat2 , 3-containing reads were identified as having RepeatMasker annotations matching “HSATII , ” “ ( CATTC ) n , ” or “ ( GAATG ) n” .', 'These represent the canonical reference sequences for HSat2 and HSat3 ( note: the latter 2 repeat types are simply reverse complements representing HSat3 ) .', 'A total of 455 , 513 reads had at least one match to these repeat types , with the lowest percent identity match of any annotated segment on any read being 61 . 4% .', 'The mean number of matched bases per read was 701 , compared to the average HuRef read length of 877 base pairs , consistent with the expectation that reads derived from long , contiguous domains of HSat2 , 3 will tend to contain HSat2 , 3 along their entire length .', 'Of these reads that have mate pairs ( ∼65% ) , 93 . 3% have mate pairs also containing HSat2 , 3 , consistent with the expectation that mate pairs derived from long contiguous domains of HSat2 , 3 will tend to both contain HSat2 , 3 .', 'To evaluate the comprehensiveness of this dataset , we employed an independent method of identifying HSat2 , 3-containing reads , in which we calculated the enrichment above background of the ancestral HSat2 , 3 5-mer , CATTC , on each read .', 'First , we counted the number of occurrences of CATTC ( or GAATG , whichever was more frequent ) on each read .', 'To model the background expectation of the number of CATTC occurrences on each read , we used a Poisson distribution with parameter 2 . 37 ( the observed mean number of CATTC on each read ) .', 'Using this background model , we selected a threshold expected to yield 1 false positive in the entire HuRef read dataset of ∼32 million reads , classifying any read with at least 14 occurrences of CATTC/GAATG as HSat2 , 3 .', 'We found the resulting dataset to be 97% concordant with the RepeatMasker dataset ( excluding only 1 . 8% of the RepeatMasker dataset and adding only 5244 new reads ) .', 'Given this high concordance , we utilized the RepeatMasker dataset for downstream analyses .', 'To examine further whether there are diverged HSat2 , 3 subfamilies not detected by the canonical RepeatMasker consensus sequences or by CATTC enrichment , the 455 , 513 candidate HSat2 , 3-containing reads were scanned for matches to any other tandem repeats present on at least 100 reads .', 'Tandem repeats that represent 1-away variants of the canonical CATTC repeat ( e . g . CACTC , CATTT , TGGA , etc . ) were then identified among all HuRef reads .', 'This yield
['Introduction', 'Results', 'Discussion', 'Methods']
['The largest gaps in the human genome assembly correspond to multi-megabase heterochromatic regions composed primarily of two related families of tandem repeats , Human Satellites 2 and 3 ( HSat2 , 3 ) .', 'The abundance of repetitive DNA in these regions challenges standard mapping and assembly algorithms , and as a result , the sequence composition and potential biological functions of these regions remain largely unexplored .', 'Furthermore , existing genomic tools designed to predict consensus-based descriptions of repeat families cannot be readily applied to complex satellite repeats such as HSat2 , 3 , which lack a consistent repeat unit reference sequence .', 'Here we present an alignment-free method to characterize complex satellites using whole-genome shotgun read datasets .', 'Utilizing this approach , we classify HSat2 , 3 sequences into fourteen subfamilies and predict their chromosomal distributions , resulting in a comprehensive satellite reference database to further enable genomic studies of heterochromatic regions .', 'We also identify 1 . 3 Mb of non-repetitive sequence interspersed with HSat2 , 3 across 17 unmapped assembly scaffolds , including eight annotated gene predictions .', 'Finally , we apply our satellite reference database to high-throughput sequence data from 396 males to estimate array size variation of the predominant HSat3 array on the Y chromosome , confirming that satellite array sizes can vary between individuals over an order of magnitude ( 7 to 98 Mb ) and further demonstrating that array sizes are distributed differently within distinct Y haplogroups .', 'In summary , we present a novel framework for generating initial reference databases for unassembled genomic regions enriched with complex satellite DNA , and we further demonstrate the utility of these reference databases for studying patterns of sequence variation within human populations .']
['At least 5–10% of the human genome remains unassembled , unmapped , and poorly characterized .', 'The reference assembly annotates these missing regions as multi-megabase heterochromatic gaps , found primarily near centromeres and on the short arms of the acrocentric chromosomes .', 'This missing fraction of the genome consists predominantly of long arrays of near-identical tandem repeats called satellite DNA .', 'Due to the repetitive nature of satellite DNA , sequence assembly algorithms cannot uniquely align overlapping sequence reads , and thus satellite-rich domains have been omitted from the reference assembly and from most genome-wide studies of variation and function .', 'Existing methods for analyzing some satellite DNAs cannot be easily extended to a large portion of satellites whose repeat structures are complex and largely uncharacterized , such as Human Satellites 2 and 3 ( HSat2 , 3 ) .', 'Here we characterize HSat2 , 3 using a novel approach that does not depend on having a well-defined repeat structure .', 'By classifying genome-wide HSat2 , 3 sequences into subfamilies and localizing them to chromosomes , we have generated an initial HSat2 , 3 genomic reference , which serves as a critical foundation for future studies of variation and function in these regions .', 'This approach should be generally applicable to other classes of satellite DNA , in both the human genome and other complex genomes .']
['algorithms', 'computer and information sciences', 'mathematics', 'genetics', 'applied mathematics', 'biology and life sciences', 'population biology', 'computing methods', 'genomics', 'computerized simulations', 'physical sciences', 'computational biology']