ARTICLE https://doi.org/10.1038/s42003-023-05004-9 OPEN Y chromosome sequence and epigenomic reconstruction across human populations Paula Esteller-Cucala 1,7 ✉, Marc Palmada-Flores1,7, Lukas F. K. Kuderna1, Claudia Fontsere1, Aitor Serres-Armero1, Marc Dabad 2, María Torralvo 1, Armida Faella1, Luis Ferrández-Peral1, Laia Llovera1, Oscar Fornas 3,4, Eva Julià 3, Erika Ramírez3, Irene González3, Jochen Hecht3, Esther Lizano 1,5, David Juan1 & Tomàs Marquès-Bonet 1,2,4,5,6 ✉ Recent advances in long-read sequencing technologies have allowed the generation and 1234567890():,; curation of more complete genome assemblies, enabling the analysis of traditionally neglected chromosomes, such as the human Y chromosome (chrY). Native DNA was sequenced on a MinION Oxford Nanopore Technologies sequencing device to generate genome assemblies for seven major chrY human haplogroups. We analyzed and compared the chrY enrichment of sequencing data obtained using two different selective sequencing approaches: adaptive sampling and flow cytometry chromosome sorting. We show that adaptive sampling can produce data to create assemblies comparable to chromosome sorting while being a less expensive and time-consuming technique. We also assessed haplogroup- specific structural variants, which would be otherwise difficult to study using short-read sequencing data only. Finally, we took advantage of this technology to detect and profile epigenetic modifications among the considered haplogroups. Altogether, we provide a fra- mework to study complex genomic regions with a simple, fast, and affordable methodology that could be applied to larger population genomics datasets. 1 Institut de Biologia Evolutiva (CSIC-Universitat Pompeu Fabra), Doctor Aiguader 88, Barcelona, Spain. 2 CNAG-CRG, Centre for Genomic Regulation (CRG), Barcelona Institute of Science and Technology (BIST), Baldiri i Reixac 4, Barcelona, Spain. 3 Centre for Genomic Regulation (CRG), Barcelona Institute for Science and Technology (BIST), Doctor Aiguader 88, Barcelona, Spain. 4 Universitat Pompeu Fabra (UPF), Doctor Aiguader 88, Barcelona, Spain. 5 Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, Cerdanyola del Vallès, Spain. 6 Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluís Companys 23, Barcelona, Spain. 7These authors contributed equally: Paula Esteller-Cucala, Marc Palmada- Flores. ✉email: paula.esteller@upf.edu; tomas.marques@upf.edu COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio 1 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 H uman sex chromosomes have been traditionally excluded et al.17, whose haplogroup (A0) represents one of the deepest- from genome-wide studies1,2. This exclusion is particu- rooting known haplogroups. Additionally, we also generated larly pronounced for the Y chromosome, the study of Illumina short-read data for the same flow-sorted chromosomes which could be key in understanding differences in disease sus- (Supplementary Data 1). ceptibility between men and women3–5. However, the Y chro- The ONT data available for the cell lines ranged from 6.4 to mosome is now considered important not only for male-specific 10.3 Gb, of which 7–33% mapped to the Y chromosome in the traits but also for the study and characterization of common reference. Moreover, we also generated 6.4–35 Gb of Illumina complex diseases4. Sex-limited chromosomes, defined as those data for all the chromosome sorting extractions. This is a notably unique to a heterogametic genome6, are usually harder to high amount of data, especially considering that the Y chromo- assemble since they are haploid and thus have half the sequencing some sequence represents less than 1% of the known sequence in depth when sequenced together with other autosomal and GRCh38. homogametic chromosomes. Moreover, their repetitive nature, The Y chromosome enrichment specificity was assessed by filled with ampliconic regions and heterochromatin, poses an aligning the basecalled data to the human reference genome additional challenge for assemblers7. assembly GRCh38 and calculating the normalized coverage on The first Y chromosome assemblies were generated by means each chromosome accounting for the gaps of the reference of bacterial artificial chromosomes (BACs), which are labor- genome and the ploidy of each chromosome (“Methods”; Fig. 1c intensive and time-consuming approaches8–10. Indeed, the Y and Supplementary Figs. 1 and 2). The Y chromosome-specific chromosome sequences in the GRCh38 assembly11–13 are a enrichment factor of the six samples showed high variability, as it composite of BAC clones11 from a male that belongs to the R1b ranged from 15- to 50-fold, whereas the A0 haplogroup was over haplogroup14 and pseudoautosomal (PAR) regions from the 100-fold enriched (Supplementary Data 2). As noted in Kuderna X-chromosome. et al.17, we found that chromosome 22 partially co-sorts with To facilitate the assembly process and also to avoid the use of chromosome Y, showing enrichment values slightly higher such costly techniques, one can decrease the potential inter- than 1. chromosomal assembly overlaps by specifically enriching the chromosome of interest. This can be done by physically isolating Adaptive sampling as a strategy to enrich specific chromo- the chromosome using flow cytometry (chromosome somes. A limiting factor in chromosome sorting is the need to sorting)6,15–17. Alternatively, other selective sequencing methods, culture hundreds of millions of cells in order to enrich the chro- such as adaptive sampling on Oxford Nanopore Technologies mosome of interest effectively16,17. To overcome this limitation, we (ONT) devices18, can potentially be used. explored the potential of adaptive sampling to specifically enrich Chromosome sorting allows the chromosome of interest to be the Y chromosome. This approach was done for one of the cell lines sequenced on different platforms after its physical isolation by (haplogroup H) for which chromosome sorting data had also been flow cytometry. This separation is possible because different generated. We used the nucleotide sequences of the Y chromosome chromosomes have specific fluorescence intensity19. On the other (chrY) and the contig chrY_KI270740v1_random (chrY_random) hand, adaptive sampling allows for the sequencing of specific as provided in the GRCh38 assembly as the target sequences to DNA regions by sequence enrichment or depletion of off-target enrich. To obtain comparable coverages using the two methodol- reads during sequencing without the need for previous chromo- ogies (~18x), we ran two ONT MinION flowcells with adaptive some enrichment20,21. Any DNA molecule that does not corre- sampling. In both experiments, we showed that the Y chromosome spond to the genomic region of interest will be ejected from the was preferentially enriched to the other chromosomes (Supple- pore, thereby preventing any further sequencing. To obtain a de mentary Fig. 3). Although the Y chromosome enrichment factor novo assembly, it is also important to avoid whole-genome value with chromosome sorting doubles the one in adaptive sam- amplification (WGA), as this process can introduce chimeras, pling for this cell line (Fig. 1c; Supplementary Data 2), adaptive bias the assembly process22, and prevent the detection of epige- sampling proves to be a cheaper and less time-consuming strategy. netic modifications. Long-read whole-genome sequencing enables the assessment of previously unsolved repeats and thus allows the generation of Y chromosome assembly across haplogroups and enrichment more contiguous assemblies. Currently, ONT can achieve the techniques. We obtained Y chromosome data and assemblies longest read lengths compared to any other existing sequencing corresponding to seven different Y chromosome haplogroups technology23–25. Moreover, ONT allows the detection of DNA using chromosome sorting on distinct cell lines (“Methods”; (and RNA) modifications based on the different current signals of Supplementary Figs. 4–6; Supplementary Data 3). The coverage the nanopores26,27. Taken together, this technology is able to used by the assembler to generate each assembly ranged from 13 resolve gaps, allowing for the true completion of chromosomes or to 50x, with a mean assembly coverage of over 28x. The resulting even genomes28–30. Here we assess the performance of two assemblies spanned from 18.95 to 22.08 Mb in length, being enrichment methods to sequence and assemble the Y chromo- 16–28% shorter than the length of chromosome Y in GRCh38. somes from seven major human haplogroups. Moreover, we We also observed that assemblies with higher continuity (contig provide insights into their structural variation and epigenomic N50) tend to have higher values of read length N50 and mean landscape, showing that enrichment techniques coupled with read lengths (Supplementary Data 3). ONT can be used to study variation between population datasets. Our assemblies had a similar amount of contigs compared to a previously published African Y chromosome assembly (hap- logroup A0)17, which is lower than the number of contigs of the Results GRCh38. The N50 across our assemblies ranged from Data production. Complete Y chromosomes from six different 1.40–2.67 Mb and are thus within the same order of magnitude human haplogroups were isolated as previously described17 as the Y chromosome in GRCh38 (6.91 Mb). These results (Fig. 1a, b). In brief, chromosomes were obtained from lym- suggest that creating de novo assemblies primarily based on long phoblastoid cell lines (LCLs) used in the 1000 Genomes Project31 reads mapping to a reference chromosomal assembly might lead (1kgp) and sequenced on the ONT MinION. We also made use of to shorter assemblies, less fragmented but with lower continuity the Y chromosome sorted ONT data generated by Kuderna values (such as lower contig N50). However, these results also 2 COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 ARTICLE Fig. 1 Study design, enrichment factor, and assemblies. a Summary of the samples, methodologies, and analyses used in the study. b Phylogenetic tree of the human Y chromosomes used in the study. Split times are taken from Jobling and Tyler-Smith79. kya, kilo years ago. c Enrichment factor values of the H haplogroup from data generated using chromosome sorting and adaptive sampling. The chrY shows higher enrichment with chromosome sorting than adaptive sampling for the haplogroup compared. The dashed vertical line equal to 1 denotes no chromosomal enrichment. For the autosomal chromosomes, the mean enrichment value is displayed, and error bars represent the standard deviation (n = 22 autosomal chromosomes). d Dot-plots of the manually scaffolded Y chromosomes compared to the resolved MSY region of GRCh38. The large-scale deletion in the J haplogroup is most likely due to its low coverage. Source data are provided in Supplementary Data 8. COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio 3 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 show that by using ONT and Illumina platforms, it is possible to reads show an overall better performance than methods based on generate assemblies almost as continuous as the current GRCh38 short-read data37. Taking advantage of our data, we assessed the reference, for which much more effort and resources were landscape of structural variants in the Y chromosome in the seven devoted11–13,32. haplogroups. For that, we used two approaches: one based on Furthermore, we manually scaffolded each haplogroup to long-read mapping (Sniffles38,39) and another based on assembly create a single scaffold based on genome-to-genome alignments comparison (Assemblytics40). to the GRCh38 Y chromosome (Fig. 1d). Comparing our First, we identified different structural variants based on how assemblies to the GRCh38 in the male-specific region of the Y the reads align to a reference genome using Sniffles. We used chrY chromosome (MSY) shows how most of the MSY sequence and the chrY_random sequence from the GRCh38 as the classes were assembled for the most part. As previously reference. After merging the indel calls (“Methods”), we identified reported17, the ampliconic region was the most fragmented and 803 unique variants (801 indels), including 166 structural variants the least complete, most likely due to collapsed repeats. Of note, (at least 50 bp in size, Supplementary Data 5). The number of for the J haplogroup, we observed a big gap in the region variants ranged from 103 to 536 events per haplogroup. comprising 6.7–9.3 Mb of the Y chromosome. This region Moreover, Sniffles detected one translocation in the H haplogroup includes an X-degenerate region and most of its adjacent and one duplication event shared between five haplogroups (all ampliconic region. When inspecting the reads mapping to this but A0 and A1a, which are basal relative to the others). The region, we found few reads present, thus possibly explaining why detected duplication is located in the position chrY:56,673,215, at we could not accurately assemble this region. the end of a gap. This indicates that the reference is missing a Compared to the previously assembled African chrY that also region of around 98,295 bp, similar to the sequence close to the made use of chromosome sorting data, we were able to generate a gap. Most of the events were indels of 10–50 bp in size (Fig. 2a). longer and more contiguous assembly starting from the same raw Out of the 801 variants found, there were 320 insertions and 481 fast5 reads (Supplementary Data 3). This demonstrates the value deletions (including three insertions and nine deletions from of combining up-to-date basecalling and assembly tools, which chrY_random). We then analyzed the repeat content of these are constantly evolving for long-read data33–35. variants and found that all indels from chrY_random consisted of Apart from chromosome sorting (CS), data from the H LTR12B-like elements (from LTR retrotransposons type ERV1) haplogroup (GM21113 cell line) was also generated using adaptive (Supplementary Data 6). In the Y chromosome, we found 116 sampling (AS). In order to generate and compare the assemblies insertions and 153 deletions with repetitive elements (RE), and in between the two enrichment methods, and given the unequal the majority of these cases, the variants were almost entirely amount of data generated between them (83.5 Mb difference), we spanned by repeats (Supplementary Data 6). restricted the comparison to assemblies generated using the same Next, we manually investigated the longest events detected number of bases (“Methods”). The resulting assemblies showed using Sniffles and confirmed the longest deletion of 6314 bp and similar values in metrics such as genome span (CS: 21.8 Mb, AS: longest insertion of 6023 bp. Both were found in haplogroup A0 22.0 Mb), contig N50 Mb (CS: 2.7 Mb, AS: 2.6 Mb), and L50 (both and belonged to different X-degenerate regions (Supplementary L50 = 3 scaffolds). Moreover, the AS-based assembly led to a Figs. 9 and 10). The longest deletion had only 38% repeat content slightly more fragmented assembly (44 sequences) compared to the but included different types of repeats, such as ALU-like elements CS-based one (31 sequences) (Table 1 and Supplementary Fig. 7). (AluSc and AluSx), L1-like elements (L1PA13), and Simple Altogether, we have generated assemblies for seven Y Repeats (CCTTn). The repeat content of the longest insertion chromosome haplogroups with similar contiguity to previously (98.7% of the sequence) only comprised one single LINE L1-like published assemblies. Moreover, we also show that adaptive element, L1HS, which is one of the few retrotransposition- sampling can be used for generating assemblies that are competent human-specific L1s41,42 (Supplementary Data 6). comparable to those generated by chromosome sorting. Recently, Even though the most common RE in indels are simple repeats a complete T2T human Y chromosome assembly was for both insertions and deletions, proportionally more insertions published36. This new chromosome assembly showed a higher were caused by smaller REs (for example, ALUs or Satellite-like), proportion of repetitive elements (RE, Supplementary Data 4), and more deletions were detected in regions with long terminal most of which were satellite sequences. Namely 85% (53 Mb) of repeats (LTRs), mostly ERV1 LTR12B-like elements (Supple- this new T2T-Y was annotated as repetitive36, considerably more mentary Data 6). Interestingly, regions enriched with repeat than the 31% (17.5 Mb) reported for GRCh38-Y and also the 68% arrays harboring LTR12B motifs have been previously reported as (13–15 Mb) we found in the assemblies we generated in this study CNV hotspots43. (Supplementary Data 4). Moreover, this assembly has resolved After merging and regenotyping the panel of indels, we the ampliconic regions (Supplementary Fig. 8). observed that the cell line belonging to haplogroup J was the one having more undetermined genotypes (i.e., positions with no genotype information). This correlates with a lower sequencing The landscape of structural variants across the human Y depth for this sample. We also observed that 14% of the variants chromosome phylogeny. As expected by the nature of these data, genotyped in all haplogroups shared the same genotype, which methods to detect structural variants which make use of long was different from the reference (Fig. 2b). Haplogroups A0 and Table 1 Y chromosome assembly statistics in different enrichment methods. Selective sequencing method to generate the long-read Assembly span (bp) Scaffold N50 (bp) Scaffold L50 Number of sequences data Adaptive sampling (AS) 21,955,745 2,612,207 3 44 Chromosome sorting (CS) 21,794,102 2,666,112 3 31 AS+CS 22,007,578 2,640,901 3 42 Assembly metrics of Y chromosome assemblies for the GM21113 cell line (haplogroup H) using adaptive sampling and chromosome sorting, and an assembly using all the data available for GM21113. 4 COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 ARTICLE Fig. 2 Profiling of structural variants. a Number of the structural variant events, insertions, and deletions called by Sniffles or Assemblytics for the different haplogroups. Variants are grouped into three categories depending on their length: from 10 up to 50 bp, from 50 up to 500 bp, and equal to or over 500 bp. b Overlap on the alternative calls between haplogroups. As expected by their evolutionary distance, haplogroups A0 and A1a show higher haplogroup- specific variants. Only variants with genotype calls for all haplogroups have been included (n = 726 variants). c Correlations between genotype calls using Sniffles (ONT-based) or graphtyper (Illumina-based) when calling the same set of structural variants. Phi coefficients range from 1 to −1, where 1 indicates complete association. Only correlation values that are statistically significant (p-value < 0.05) after Bonferroni multiple testing correction are shown. Source data are provided in Supplementary Data 8. A1a harbor the most haplogroup-specific variants, concordant to We observe 194 to 406 insertions per haplogroup with their genetic distance to the reference. We also manually assessed Assemblytics40 compared to the 45 to 255 insertions detected previously reported events for the A0 haplogroup17 and with Sniffles38,39, and 206 to 424 deletions against 57 to 288, confirmed that they were restricted to this cell line (Supplemen- respectively. However, similar amounts of structural variant tary Figs. 11 and 12). This indicates that structural variants found indels are detected by Sniffles in all haplogroups (between 52 and in only one haplogroup might not be representative of widespread 110) but for the J haplogroup (21) compared to Assemblytics structural variants of a chromosome but rather delimited to one (between 47 and 80, J haplogroup having 55 structural variants). specific population or group of individuals. These results, together with the fact that the J haplogroup is the Second, we identified structural variants based on the one with less data generated, suggest that mapping-based comparison of the obtained chrY assemblies to the reference structural variation detection methods may not be able to detect chrY GRCh38 using Assemblytics v1.2.140 (Supplementary as many structural variants compared to assembly-to-assembly Data 7). This allowed for the detection of 557 to 1019 putative comparison-based methods when having limited sequencing variants, for which 202 to 258 were at least 50 bp in size (Fig. 2a). depth. In that situation, generating a de novo assembly and We also found between 1 and 4 structural variants bigger than using Assemblytics can lead to the identification of larger indels. 50,000 bp for the cell lines studied, a type of variant that the Moreover, Assemblytics categorizes structural variants into mapping-based method may not detect because it would require multiple types including indels and other specific variants, such constant coverage along a long region of the reference. as tandem and repeat expansions or contractions, while Sniffles COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio 5 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 was only able to capture one duplication event and one chromosomes50 (Supplementary Fig. 18). X-degenerate regions translocation, besides indels. showed the characteristic bimodal distribution of frequency We further assess if the indels found using long reads could be values with a median close to 0.7, whereas all other regions similarly genotyped using short-read data. For that, we genotyped showed much less defined distributions. We also inspected the the variants confidently called by the ONT mapping-based behavior of methylation according to the epigenetic annotation of approach using Illumina data generated for the same cell line the CpG of each of the sequence classes. For that, we divided the extractions. With the Illumina data, we were able to replicate over CpGs into four mutually exclusive categories (Fig. 3b and one-quarter of the indels found in the nanopore data (214 out of Supplementary Figs. 19 and 20): those in CpG islands (CGI), the 801 indels). A significant positive association was seen CpG shores, CpG shelves, and other inter-CGI regions (open sea). between the predicted genotypes using ONT to those observed CGI in the X-degenerate and X-transposed regions were using Illumina data for each cell line (Fig. 2c). The observed phi predominantly unmethylated, while all the other regions were coefficients (correlation values for binary variables) range mostly methylated. Open sea regions showed intermediate between 0.35 and 0.54, not close to the highest correlation value methylation levels for all sequence classes but the X-degenerate, of one44,45. This is expected, given that most of the variants in our whose median 5mC frequency reached 0.75. As expected by the panel cannot be called by the Illumina data. The intra-haplogroup dynamic nature of the human methylome, CpG shores and correlation is generally higher than that inter-haplogroup. shelves showed intermediate values transitioning from CGI and However, when testing inter-haplogroup associations for geno- open sea regions (Supplementary Figs. 21–23)51–53. types derived from the same methodology, strong associations DNA methylation is associated with gene expression54, and in were also detected (Supplementary Fig. 13). To explore which of addition to the epigenetic annotation, we also inspected the 5mC the indels could be observed in a panel of human variation, we frequency patterns across different gene annotations (Supple- genotyped the same indels in the male samples present in the mentary Figs. 24 and 25). Most annotated genes are present in the 1000 Genomes Project (“Methods”). As expected, intra- X-degenerate and ampliconic sequence classes (Supplementary haplogroup associations were typically positive and significant, Fig. 26), and consistent with the different expression profiles of generally having stronger associations than inter-haplogroup the genes in LCLs (retrieved from GTEx55) in each of these two comparisons, while having a limited correlation due to the fact of sequence classes, we observed clear distinct methylation patterns comparing Illumina vs ONT calls. (Supplementary Fig. 13). The in their TSS, UTRs, and intragenic CpGs (Fig. 3c). Not use of additional cell line replicates for each haplogroup might surprisingly, we found 5’UTRs to be the most constrained gene have yielded higher associations, as it would have helped discard feature across samples, which would directly link its methylation instances of variation that arose during immortalization. status to gene expression (Supplementary Fig. 27). Moreover, we found a direct relationship between upstream CGI methylation status with gene expression (Supplementary Fig. 26). Finally, we CpG methylation across the phylogeny. ONT sequencing relies explored those cases in which differential methylation could have on the identification of different current signals when the DNA an effect on gene expression. We encountered a region with high passes through the pore, so it is possible to go beyond the iden- methylation dispersion fully spanning a protein-coding gene tification of the four canonical nucleotide bases and detect other (Fig. 3d and Supplementary Fig. 28). In that location, haplogroup modifications in the DNA. As such, we studied the 5mC land- A1a was found to be undermethylated compared to the other scape in the seven Y chromosome haplogroups. We used nano- haplogroups, and although this difference was only modest, it polish v1.1246 to call the methylation status of 5-methylcytosines could potentially modulate the expression levels of the gene (5mC) at CpG positions from the nanopore current signal. located in this region. This gene is NLGN4Y, which is a long gene Assessing the Y chromosome methylome using long reads is that spans over 300 kb and is expressed in the brain and other beneficial for exploring regions that are traditionally inaccessible tissues, including LCLs (τNLGN4Y = 0.714). Interestingly, this gene using short-read techniques, such as the PAR, X-transposed has been proposed as a candidate for autism spectrum regions, and even the ampliconic regions. disorder56,57. As expected, we found CGIs located upstream of For that, we performed quantile normalization on the this gene to be unmethylated (CGI_1 and CGI_2), which would methylation values across samples with a minimum coverage of be consistent with the expression of this gene in LCLs. On the 4x (Supplementary Fig. 14). We observed consistent methylation contrary, a CGI potentially regulating an overlapping non-coding patterns along Y chromosomes across samples, indicating a gene in the opposite strand and with no expression in LCLs was strong overall correlation on the methylation status (Fig. 3a and shown to be fully methylated in all cell lines (CGI_3). Supplementary Figs. 15 and 16). These results alone already Experimental validation of the expression of these genes could demonstrate the potential of ONT sequencing in studying DNA ultimately confirm whether such differences in methylation truly methylation in challenging genomic regions. However, 5mC affect gene expression levels or are merely coincidental. frequency values could not recapitulate the expected phylogeny, Altogether, we show that ONT can be used to study 5mC either chromosome-wise or segregating by sequence class or across different cell lines, as it consistently recapitulates the main epigenetic annotation (Supplementary Fig. 17). Given that methylation patterns observed across mammalian genomes. methylation levels might vary within the population, age, Moreover, we also show that it can be helpful for the study of environmental exposures, and cell culture conditions47–49, and traditionally challenging genomic regions, particularly those the absence of replicates for each of the haplogroups considered, present in the Y chromosome, although future work including this observation could be due to differences in any of these replicates will clarify the actual sensitivity of this approach to variables. However, given the uncertainties about the cell lines detect differentially methylated regions. and age of the individuals from which they were generated, we are unable to discern the 5mC variation, which accounts for the different haplogroups from that which could be caused by other Discussion factors. As expected by the nature of the sequence classes11, the Here, we present a panel of ONT data for seven cell lines that X-degenerate region, which harbors single-copy genes and mostly represent the major human Y chromosome haplogroups. We ubiquitous expression, showed 5mC frequency values which have generated assemblies for each of them and studied their resembled most of those normally seen in mammalian autosomal diversity, focusing on structural variation and methylation. To 6 COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 ARTICLE Fig. 3 Methylation landscape across the Y chromosome phylogeny. a Frequency of 5mC in the seven cell lines along the resolved MSY of the GRCh38. The methylation levels are calculated as the median 5mC frequency value in 250 kb sliding windows for each cell line. The sequence classes, the genes annotated, and the standard deviation (SD) of the methylation levels across cell lines are also shown. The standard deviation of the 5mC frequency is represented in a white-to-black scale, in which a darker color denotes a higher standard deviation value. b Median methylation value per cell line segregated by CpG annotation and sequence classes. CpG annotations are mutually exclusive regions that comprise: CpG islands (CGI), CpG shores (up to 2 kb away from the end of the CGI), CpG shelves (up to 2 kb away from the end of the CpG shores), and inter-CGI or open sea regions (where all remaining CpG are allocated). c 5mC frequencies on different gene features in X-degenerate and ampliconic sequence classes. Gene annotation features shown are TSS (region of 200 bp upstream of the transcription start site), both UTRs, and intragenic regions (which combine all exonic and intronic regions without considering the first gene exon). Extended version with the number of CpG sites in Supplementary Figs. 24 and 25. d Methylation frequencies in 3 CpG islands (CGI) surrounding the NLGN4Y and NLGN4Y-AS1 genes. Empty circles show the mean 5mC frequency per CGI, whereas smaller colored points indicate the individual value in each cell line. Source data are provided in Supplementary Data 8. generate this resource, chromosome sorting data were employed comparable assemblies, while they require different time, cost, and compared to adaptive sampling data, an enrichment tech- and expertise. In terms of enrichment factor values, chromosome nique that is compatible with ONT sequencing data. After gen- sorting shows co-enrichment with chromosome 22. This is erating and comparing the assemblies of the two enrichment mainly due to the fact that both chromosomes have similar sizes. techniques, we showed that both methods could lead to However, this is not seen in adaptive sampling. In fact, samples COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio 7 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 enriched with adaptive sampling show the lowest standard detect longer genomics variants. We consider that for low data deviation of the enrichment factor on autosomal chromosomes. samples, the creation of a de novo Y chromosome assembly may Nevertheless, given the homology of the sex chromosomes, and allow the detection of structural variants that cannot be recog- the fact that adaptive sampling is performed by providing the nized with a mapping-based method, considering the low cov- genomic sequence of chromosome Y, chromosome X shows a erage of reads mapping in those regions. Of note, we are aware higher enrichment factor compared to the other chromosomes that limited by only one replicate per haplogroup, we are not able and the samples enriched by chromosome sorting. Altogether we to discriminate sample-specific variation than that appearing show that adaptive sampling is a viable alternative strategy for the during the immortalization of these cells60. Moreover, the recent enrichment of specific genomic regions. We also emphasize the T2T Y chromosome36 will allow for a more reliable and com- importance of using high molecular weight DNA or long DNA prehensive assessment of the structural variants found in fragments, which are especially convenient for the enrichment of this study. small chromosomes with adaptive sampling. As such, at longer Besides the potential to generate high-accuracy assemblies and DNA fragment sizes, the time the sequencer will be scanning for resolve complex genomic regions like structural variants, ONT on-target regions (i.e., those that belong to the Y chromosome) also allows for studying the epigenome. We have assessed the will be reduced. Therefore, we realize that having started from methylation status of cytosines in a CpG context for our panel of higher DNA fragment sizes for the haplogroup H sample would cell lines. Despite the fact that the epigenome of the Y chromo- have led to higher enrichment efficiencies in the adaptive sam- some has not been deeply studied, we were able to consistently pling enrichment method. replicate the methylation patterns that have been described in One major limitation of our work is the conservative filtering other human autosomal chromosomes61,62. Not surprisingly, we have used to generate the assemblies. Our approach uses the Y with the methylation values obtained, we were unable to recapi- chromosome of the current genome of reference GRCh38 as a tulate our samples’ expected phylogeny. Two main factors can be backbone. All data obtained using adaptive sampling relies attributed to this: the lack of replicates for each haplogroup and heavily on the GRCh38 reference and may include a few reads of also within-population variability48,49,63,64, which, in our case, other chromosomes that start with a similar sequence. On the could also be confounded by epigenetic drift65. Still, methylation other side, chromosome sorting produces data on unresolved differences at the population level are expected to be small in chromosomal regions but includes some undesired full chromo- magnitude48. In this line, we were able to detect small differences some data. As such, restricting our assembly to only those reads in methylation in regions that could influence the regulation of that map to the reference leads to the loss of a fraction of Y specific genes. This is the case of gene NLGN4Y, which we found chromosome potentially informative reads during the filtering to fully overlap with a region with consistently lower methylation process. Conversely, this approach minimizes the retention of in the cell line belonging to haplogroup A1a. However, gene non-chromosome Y data and limits the resulting assembly to the expression data from these cell lines could ultimately reveal Y chromosome only. Compared to a previous assembly created whether such subtle variations in methylation translate into dif- with the same data for the cell line that belongs to haplogroup A0, ferences in gene expression. In the same line, the addition of our approach yielded a more contiguous assembly. As such, it biological replicates for different haplogroups could help uncover shows the potential that re-processing the same raw data with whether these differences are either cell-line or haplogroup- novel approaches might have in the future, especially in the specific. context of the big data era58,59. Nevertheless, we are using lymphoblastoid cell lines (LCLs), Due to its large fraction of heterochromatin, around half of the which are artificially transformed cells, so caution must be taken sequence in the current Y chromosome assembly is unresolved. when extrapolating these findings. But the extent to which the This limitation, together with the fact that we are using a partial generalization of our results could be biased is even more con- reference genome to generate assemblies of a specific chromo- sequential when reporting those findings that are sample-specific. some, hampers the possibility of reconstructing the totality of this As such, an increase in the number of replicates would help to chromosome. In the future, telomere-to-telomere Y chromosome discern which of our findings are artifacts from those which have sequencing would undoubtedly avoid reference biases we a true biological meaning. encountered in this study28. Taken together, here we provide a framework to study complex Due to its large fraction of heterochromatin, around half of the genomic regions. We applied this simple, fast, and affordable sequence in GRCh38 remains unresolved. This limitation, toge- technology to study diverse human population groups. Moreover, ther with the fact that we are using this partial reference genome this approach can be applied to the generation of long-read data to generate assemblies of the latter, hampers the possibility of of other regions or chromosomes of interest. Long-read tech- reconstructing the totality of this chromosome. The recently nologies are achieving longer and more accurate reads, which, published telomere-to-telomere Y chromosome assembly (T2T- together with higher throughputs, will enable more reliable and Y) would undoubtedly avoid reference biases we encountered in comprehensive comparative genomics studies, particularly for this study28. Our assemblies align to this T2T Y chromosome chromosomes that contain high repeat content. As such, it could assembly in more than 89.6% of their sequence while covering be used for the characterization of virtually any species, although more than 68% of the repetitive content. This is slightly more it would be especially advantageous for those rich in complex than what is found in GRCh38 (Supplementary Data 4). Such genomic features. higher repeat content recovery showcases that reference-based approaches like the one followed in this study can include repe- titive regions not present in the reference used. Methods We also took advantage of the long-read data generated to Flow chromosome sorting followed by ONT or Illumina sequencing. Chro- explore the landscape of structural variants in each cell line. For mosome preparation was performed as previously described16,17 in six lympho- blastoid cell lines (purchased from Coriell, see Supplementary Data 1 for specific that, we used two different methods for structural variant calling: details on the cell lines). The libraries to obtain the Illumina paired-end data were one based on long-read mapping and another based on assembly constructed using a SureSelect V6-Post Library Kit. Raw data generated for the comparison. The former allows for two rounds of genotyping, haplogroup A0 (HG02982 cell line) was retrieved from Kuderna et al.17. The data and so the final candidates are potentially more curated. The generated in each MinION run was basecalled using Guppy v5.0.1534 with the super accuracy model dna_r94.1_450bps_sup. latter is based on genome-to-genome comparisons, so it is able to 8 COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 ARTICLE Adaptive sampling for the enrichment of a specific chromosome. We extracted Structural variant detection with long reads. Structural variation was called DNA from cultured cells of haplogroup H (GM21113 cell line) using the Qiagen using Sniffles v2.0.238,39 with a minimum number of reads that support an SV of -s MagAttract HMW Kit. DNA libraries for ONT sequencing were obtained using the 10, fed with the bam files for which we calculated and added MD tags using Ligation Sequencing Kit (SQK-LSK110) and sequenced in two ONT MinION flow- SAMTOOLS v1.967, with the program samtools calmd adding options -uAr -Q. cells (FLO-MIN106 R9.4.1) using a MinION Mk1C with MinKNOW v21.02- We summarized the number of SVs per type and filtered out the SVs considered beta4~xenial. We aimed for the specific enrichment of the chrY and the chrY_ran- ‘IMPRECISE’ by Sniffles. dom by adding their nucleotide sequence as provided in the GRCh38 assembly. This We merged the insertions and deletion separately with a maximum permitted method bioinformatically labels the reads that are being sequenced for enrichment or distance of 100 bp (so that indels located 100 bp upward or downward will be depletion. After a DNA strand enters the pore, the sequencer only needs one second considered a single event) found independently in all the cell lines using (around 420 bases) to decide whether to continue sequencing the DNA if it matches SURVIVOR v1.0.775. We removed any genotype with quality under 25 (MQ) and the region of interest or to eject it if it does not. Each of the strands that enter a pore the events that were homozygous for the reference genotype in all samples. will be labeled as unblock and no decision when they are rejected by the pore or they We used Assemblytics v1.2.140 to find structural variants in the different are so short that their status remains inconclusive, respectively. They will be labeled as assemblies generated by comparing them to the reference GRCh38. We looked for stop receiving when they are on target, thus further sequenced. Only reads labeled as structural variants with sizes between 10 and 100,000 and the unique sequence stop receiving were used in this project. length required to call a variant of 1000. The enrichment obtained with adaptive sampling highly depends on the fragment length of the library. Longer DNA library lengths are preferred, as the adaptive sampling enrichment algorithm takes a fixed amount of time to recognize Indel composition. We used RepeatMasker v4.1.2-p176 with options -s -species whether to enrich a DNA strand. Because of this, in order to target a specific region human over the different assemblies and indel reference panel to annotate repe- of the genome that is particularly small (the Y chromosome represents ~1% of the titive sequences. genome), it will always be better to have few long DNA fragments rather than many short DNA fragments, as the time spent by the sequencer scanning for on- target regions will be reduced. Structural variant genotyping with short reads. We genotyped, using the indels obtained using Sniffles and SURVIVOR as a reference, the structural variants based on Illumina data with the program graphtyper v2.7.5 with the option “genotype_sv” Assessing the performance of two different enrichment methodologies. The and only kept the indels with a quality >0. We genotyped them with the Illumina coverage and enrichment factor for each chromosome were calculated as follows: data generated in this study and with the Illumina data of the 1kgp31 available for the Y chromosome. Mapped bp in chrN Coverage of chrN ¼ ð1Þ Size chrN ðbp; without NÞ Correlation between structural variant detection using long or short reads. To assess the reproducibility of the structural variant calls obtained with ONT data in short-read data. We took the indels genotyped in the Illumina data and compared Mapped bp in chrN them between platforms. For each structural variant, the ONT genotypes were Enrichment factor of chrN ¼ Total mapped bp ð2Þ assumed to be true positives, and all genotype calls were binarised into presence (1) Size chrN ´ n ðbp; without NÞ or absence (0). For each structural variant, and given that Y chromosomes are Diploid genome size ðbp; without NÞ hemizygous, homozygous and heterozygous alternative calls were considered Since more than 50% of the Y chromosome in GRCh38 is composed of long present, and homozygous reference genotypes absent. stretches of unknown sequence (that in the assembly is seen as N), it is important As we wanted to study the correlation of binary variables, we made use of phi to exclude these regions from the coverage and enrichment calculations. Because of coefficients, also known as Matthews correlation coefficient or MCC. Phi that, Eq. 1 and Eq. 2 only consider chromosome sizes without Ns. Moreover, for coefficients should be interpreted similarly to a Pearson correlation coefficient. For calculating the enrichment and in order to account for the real target space of each the 1kgp data, we considered a structural variant to be present if it was at a chromosome, the size of each of them is multiplied by its ploidy. frequency higher than 0.2. Assembly generation. Basecalled passed reads (Q > 10) were mapped to the Studying methylation using ONT. The methylation status was called using human GRCh38 genome assembly using minimap2 v2.17-r94166 with the option nanopolish v0.13.246, which assigns a log-likelihood ratio to each individual CpG -x map-ont. The resulting bam was indexed using SAMTOOLS v1.1266,67, and site. To avoid adding noise to the methylation results, we only used reads with the the reads mapping either to chrY or chrY_random (chrY-specific reads) highest mapping quality as provided by minimap2 (mapQ = 60) and filtered out all were retrieved. others. We used the default log-likelihood threshold of 2 as implemented in We ran Flye v2.933 using the chrY specific reads with the option --nano-hq as nanopolish v0.12 onward. As suggested by the developers, we called methylation suggested by the developers while using data basecalled using Guppy v534 onward with the option --min-separation 5 to help calling CpG dense regions. The with the super accuracy model. We added the option --scaffold to enable methylation frequency was calculated for each site as the number of mapped reads scaffolding based on the assembly graph and included two internal rounds of predicted as methylated divided by the number of total mapped reads. polishing with the argument -i 2. We filtered out the few instances in which alternative alleles were present in a As we used uncorrected long reads to obtain the draft assemblies, we polished the genomic position with cytosine in the reference sequence. We performed quantile initial assemblies by first using ONT reads. We started with two rounds of Racon normalization on the methylation values across samples with a minimum coverage v1.3.168, using minimap2 v2.9-r72066 with the option -x ont to obtain the mapping file of 4x using the R package preprocessCore v1.56.077. CpG and gene annotations were and adding to Racon the argument -u to keep any unpolished sequences. To further obtained using the R package annotatr v1.24.078. Minor modifications were made improve the assembly, we then ran medaka v1.4.169 using the medaka_consensus to these annotations for different analyses. All these modifications have been program with default settings and the model -m r941_prom_sup_g507. specifically described when used in the text. For the overlapping regions in the gene Additionally, to polish the assemblies with Illumina data, we used HyPo annotations, the priority set was the following: promoters, UTRs (5’, 3’), first exon, v1.0.370, mapping the Illumina reads to the polished assembly. For mapping non-first exons, all introns, and upstream region. short reads to the existing assembly, we used minimap2 v2.9-r72066 with the option -x sr. Reporting summary. Further information on research design is available in the Nature Once we polished the assemblies, we purged them using purge_dups v1.2.566,71, Portfolio Reporting Summary linked to this article. with default parameters and the -2 option. This was done to remove any haplotig present in the assemblies. We obtained the mapping files using minimap2 v2.14- r88366 with the option -x ont to map the ONT reads to the polished assembly and with the options -xasm5 -DP to map the split polished assembly to itself. Data availability For the comparison of the assemblies generated from the two selective All sequencing data generated for this study have been deposited at the European sequencing methods (chromosome sorting and adaptive sampling), we Nucleotide Archive (ENA) under the study accession PRJEB58141. Assemblies are downsampled the data of the adaptive sampling experiment. For that, we used deposited at the ENA under the study accession PRJEB59245. Raw sequencing data for Filtlong v0.2.072 with the option --keep_percent 87.8 so as to retrieve 87.8% of the the A0 haplogroup (cell line HG02982) were retrieved from the ENA study accession AS sequencing data. From that point on, the assembly process was the same as the PRJEB28143 and its assembly from the accession ULGL01000000. The source data one explained (Supplementary Fig. 5). underlying Figs. 1–3 are provided in Supplementary Data 8. All other relevant data are available upon request. Genome-to-genome comparisons. To obtain genome-to-genome alignments, we used MuMmer v3.2373 nucmer tool with options –-maxmatch -l 100 -c 100. To Received: 21 December 2022; Accepted: 31 May 2023; manually scaffold the Y chromosome assemblies, we used the dot-plot viewer dot74. We manually reordered and reoriented the scaffolds, joining them based on the alignments to the GRCh38 reference. COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio 9 ARTICLE COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 References 35. Dida, F. & Gangman, Yi. Empirical evaluation of methods for de novo genome 1. Accounting for sex in the genome. Nat. Med. 23, 1243 https://doi.org/10.1038/ assembly. PeerJ Comput. Sci. 7, e636 (2021). nm.4445 (2017). 36. Rhie, A. et al. The complete sequence of a human Y chromosome. Preprint at 2. Wise, A. L., Gyi, L. & Manolio, T. A. eXclusion: toward integrating the X bioRxiv https://doi.org/10.1101/2022.12.01.518724 (2022). chromosome in genome-wide association analyses. Am. J. Hum. Genet. 92, 37. Mahmoud, M. et al. Structural variant calling: the long and the short of it. 643–647 (2013). Genome Biol. 20, 246 (2019). 3. Wilson, M. A. The Y chromosome and its impact on health and disease. Hum. 38. Sedlazeck, F. J. et al. Accurate detection of complex structural variations using Mol. Genet. 30, R296–R300 (2021). single-molecule sequencing. Nat. Methods 15, 461–468 (2018). 4. Anderson, K., Cañadas-Garre, M., Chambers, R., Maxwell, A. P. & McKnight, 39. Smolka, M. et al. Comprehensive structural variant detection: from mosaic to A. J. The challenges of chromosome Y analysis and the implications for population-level. Preprint at bioRxiv https://doi.org/10.1101/2022.04.04. chronic kidney disease. Front. Genet. 10, 781 (2019). 487055 (2022). 5. Molina, E., Clarence, E. M., Ahmady, F., Chew, G. S. & Charchar, F. J. 40. Nattestad, M. & Schatz, M. C. Assemblytics: a web analytics tool for the Coronary artery disease: why we should consider the Y chromosome. Heart detection of variants from an assembly. Bioinformatics 32, 3021–3023 Lung Circ. 25, 791–801 (2016). (2016). 6. Mank, J. E. The W, X, Y and Z of sex-chromosome dosage compensation. 41. Beck, C. R. et al. LINE-1 retrotransposition activity in human genomes. Cell Trends Genet. 25, 226–233 (2009). 141, 1159–1170 (2010). 7. Tomaszkiewicz, M., Medvedev, P. & Makova, K. D. Y and W chromosome 42. Seczynska, M. & Lehner, P. J. The sound of silence: mechanisms and assemblies: approaches and discoveries. Trends Genet. 33, 266–282 (2017). implications of HUSH complex function. Trends Genet. 39, 251–267 (2023). 8. Hughes, J. F. et al. Chimpanzee and human Y chromosomes are remarkably 43. Poznik, G. D. et al. Punctuated bursts in human male demography inferred divergent in structure and gene content. Nature 463, 536–539 (2010). from 1,244 worldwide Y-chromosome sequences. Nat. Genet. 48, 593–599 9. Hughes, J. F. et al. Strict evolutionary conservation followed rapid gene loss on (2016). human and rhesus Y chromosomes. Nature 483, 82–86 (2012). 44. Chicco, D. & Jurman, G. The advantages of the Matthews correlation 10. Soh, Y. Q. S. et al. Sequencing the mouse Y chromosome reveals convergent coefficient (MCC) over F1 score and accuracy in binary classification gene acquisition and amplification on both sex chromosomes. Cell 159, evaluation. BMC Genomics 21, 6 (2020). 800–813 (2014). 45. Yule, G. U. On the methods of measuring association between two attributes. 11. Skaletsky, H. et al. The male-specific region of the human Y chromosome is a J. R. Stat. Soc. 75, 579 (1912). mosaic of discrete sequence classes. Nature 423, 825–837 (2003). 46. Loman, N. J., Quick, J. & Simpson, J. T. A complete bacterial genome 12. Lander, E. S. et al. Initial sequencing and analysis of the human genome. assembled de novo using only nanopore sequencing data. Nat. Methods 12, Nature 409, 860–921 (2001). 733–735 (2015). 13. International Human Genome Sequencing Consortium. Finishing the 47. Varley, K. E. et al. Dynamic DNA methylation across diverse human cell lines euchromatic sequence of the human genome. Nature 431, 931–945 (2004). and tissues. Genome Res. 23, 555–567 (2013). 14. Mendez, F. L., David Poznik, G., Castellano, S. & Bustamante, C. D. The 48. Fraser, H. B., Lam, L. L., Neumann, S. M. & Kobor, M. S. Population- divergence of Neandertal and modern human Y chromosomes. Am. J. Hum. specificity of human DNA methylation. Genome Biol. 13, R8 (2012). Genet. 98, 728–734 (2016). 49. Husquin, L. T. et al. Exploring the genetic basis of human population 15. Tomaszkiewicz, M. et al. A time- and cost-effective strategy to sequence differences in DNA methylation and their causal impact on immune gene mammalian Y Chromosomes: an application to the de novo assembly of regulation. Genome Biol. 19, 222 (2018). gorilla Y. Genome Res 26, 530–540 (2016). 50. Bell, J. T. et al. DNA methylation patterns associate with genetic and gene 16. Kuderna, L. F. K. et al. Flow sorting enrichment and nanopore sequencing of expression variation in HapMap cell lines. Genome Biol. 12, R10 (2011). chromosome 1 from a Chinese individual. Front. Genet. 10, 1315 (2019). 51. Illingworth, R. S. & Bird, A. P. CpG islands—‘a rough guide’. FEBS Lett. 583, 17. Kuderna, L. F. K. et al. Selective single molecule sequencing and assembly of a 1713–1720 (2009). human Y chromosome of African origin. Nat. Commun. 10, 4 (2019). 52. Irizarry, R. A. et al. The human colon cancer methylome shows similar hypo- 18. Martin, S. et al. Nanopore adaptive sampling: a tool for enrichment of low and hypermethylation at conserved tissue-specific CpG island shores. Nat. abundance species in metagenomic samples. Genome Biol. 23, 11 (2022). Genet. 41, 178–186 (2009). 19. Doležel, J. et al. Chromosomes in the flow to simplify genome analysis. Funct. 53. Ziller, M. J. et al. Charting a dynamic DNA methylation landscape of the Integr. Genomics 12, 397–416 (2012). human genome. Nature 500, 477–481 (2013). 20. Payne, A. et al. Readfish enables targeted nanopore sequencing of gigabase- 54. Lowdon, R. F., Jang, H. S. & Wang, T. Evolution of epigenetic regulation in sized genomes. Nat. Biotechnol. 39, 442–450 (2021). vertebrate genomes. Trends Genet. 32, 269–283 (2016). 21. Kovaka, S., Fan, Y., Ni, B., Timp, W. & Schatz, M. C. Targeted nanopore 55. GTEx Consortium. The GTEx Consortium atlas of genetic regulatory effects sequencing by real-time mapping of raw electrical signal with UNCALLED. across human tissues. Science 369, 1318–1330 (2020). Nat. Biotechnol. 39, 431–441 (2021). 56. Ross, J. L., Tartaglia, N., Merry, D. E., Dalva, M. & Zinn, A. R. Behavioral 22. Pinard, R. et al. Assessment of whole genome amplification-induced bias phenotypes in males with XYY and possible role of increased NLGN4Y through high-throughput, massively parallel whole genome sequencing. BMC expression in autism features. Genes Brain Behav. 14, 137–144 (2015). Genomics 7, 216 (2006). 57. Chen, J., Yu, S., Fu, Y. & Li, X. Synaptic proteins and receptors defects in 23. Udaondo, Z. et al. Comparative analysis of PacBio and Oxford Nanopore autism spectrum disorders. Front. Cell. Neurosci. 8, 276 (2014). sequencing technologies for transcriptomic landscape identification of 58. Dall’Alba, G. et al. A survey of biological data in a big data perspective. Big Penaeus monodon. Life 11, 862 (2021). Data 10, 279–297 (2022). 24. Lang, D. et al. Comparison of the two up-to-date sequencing technologies for 59. Kamble, S. S., Gunasekaran, A., Goswami, M. & Manda, J. A systematic genome assembly: HiFi reads of Pacbio Sequel II system and ultralong reads of perspective on the applications of big data analytics in healthcare Oxford Nanopore. Gigascience 9, giaa123. management. Int. J. Healthc. Manag. 12, 226–240 (2019). 25. Tvedte, E. S. et al. Comparison of long-read sequencing technologies in 60. Nickles, D. et al. In depth comparison of an individual’s DNA and its interrogating bacteria and fly genomes. G3 11, jkab083 (2021). lymphoblastoid cell line using whole genome sequencing. BMC Genomics 13, 26. Rand, A. C. et al. Mapping DNA methylation with high-throughput nanopore 477 (2012). sequencing. Nat. Methods 14, 411–413 (2017). 61. McCarthy, N. S. et al. Meta-analysis of human methylation data for evidence 27. Simpson, J. T. et al. Detecting DNA cytosine methylation using nanopore of sex-specific autosomal patterns. BMC Genomics 15, 981 (2014). sequencing. Nat. Methods 14, 407–410 (2017). 62. Johansson, A., Enroth, S. & Gyllensten, U. Continuous aging of the human DNA 28. Nurk, S. et al. The complete sequence of a human genome. Science 376, 44–53 (2022). methylome throughout the human lifespan. PLoS ONE 8, e67378 (2013). 29. Logsdon, G. A. et al. The structure, function and evolution of a complete 63. Palumbo, D., Affinito, O., Monticelli, A. & Cocozza, S. DNA methylation human chromosome 8. Nature 593, 101–107 (2021). variability among individuals is related to CpGs cluster density and 30. Miga, K. H. et al. Telomere-to-telomere assembly of a complete human X evolutionary signatures. BMC Genomics 19, 229 (2018). chromosome. Nature 585, 79–84 (2020). 64. Galanter, J. M. et al. Differential methylation between ethnic sub-groups 31. 1000 Genomes Project Consortium. A global reference for human genetic reflects the effect of genetic ancestry and environmental exposures. Elife 6, variation. Nature 526, 68–74 (2015). e20532 (2017). 32. Tilford, C. A. et al. A physical map of the human Y chromosome. Nature 409, 65. Horvath, S. & Raj, K. DNA methylation-based biomarkers and the epigenetic 943–945 (2001). clock theory of ageing. Nat. Rev. Genet. 19, 371–384 (2018). 33. Kolmogorov, M., Yuan, J., Lin, Y. & Pevzner, P. A. Assembly of long, error- 66. Li, H. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics prone reads using repeat graphs. Nat. Biotechnol. 37, 540–546 (2019). 34, 3094–3100 (2018). 34. Wick, R. R., Judd, L. M. & Holt, K. E. Performance of neural network basecalling 67. Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics tools for Oxford Nanopore sequencing. Genome Biol. 20, 129 (2019). 25, 2078–2079 (2009). 10 COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio COMMUNICATIONS BIOLOGY | https://doi.org/10.1038/s42003-023-05004-9 ARTICLE 68. Vaser, R., Sović, I., Nagarajan, N. & Šikić, M. Fast and accurate de novo E.L. provided analytical support. M.D. and M.T. helped with data curation and analyses. genome assembly from long uncorrected reads. Genome Res. 27, 737–746 M.P.-F. generated assemblies and structural variant calls. P.E.-C. performed methylation (2017). analyses, supervised all other analyses and generated most figures, including all main 69. nanoporetech/medaka: sequence correction provided by ONT Research. figures. D.J. designed and supervised analyses. P.E.-C. and M.P.-F. wrote the manuscript GitHub https://github.com/nanoporetech/medaka (2018). with input from all co-authors. 70. Kundu, R., Casey, J. & Sung, W.-K. HyPo: super fast & accurate polisher for long read genome assemblies. Preprint at bioRxiv https://doi.org/10.1101/ 2019.12.19.882506 (2019). Competing interests L.F.K.K. is currently an employee of Illumina Inc. All other authors declare no competing 71. Guan, D. et al. Identifying and removing haplotypic duplication in primary interests. genome assemblies. Bioinformatics 36, 2896–2898 (2020). 72. Wick, R. rrwick/Filtlong. GitHub https://github.com/rrwick/Filtlong (2018). 73. Kurtz, S. et al. Versatile and open software for comparing large genomes. Additional information Genome Biol. 5, R12 (2004). Supplementary information The online version contains supplementary material 74. MariaNattestad/dot: Dot: an interactive dot plot viewer for comparative available at https://doi.org/10.1038/s42003-023-05004-9. genomics. GitHub https://github.com/marianattestad/dot (2017). 75. Jeffares, D. C. et al. Transient structural variations have strong effects on Correspondence and requests for materials should be addressed to Paula Esteller-Cucala quantitative traits and reproductive isolation in fission yeast. Nat. Commun. 8, or Tomàs Marquès-Bonet. 14061 (2017). 76. Smit, A. F. A., Hubley, R. & Green, P. RepeatMasker Open-4.0. 2013–2015. Peer review information Communications Biology thanks the anonymous reviewers for http://www.repeatmasker.org (2015). their contribution to the peer review of this work. Primary handling editors: Hélène 77. Bolstad, B. preprocessCore: a collection of pre-processing functions. R package Choquet, Luke R. Grinham, and George Inglis. A peer review file is available. version 1.62.1. https://doi.org/10.18129/B9.BIOC.PREPROCESSCORE (Bioconductor, 2017). Reprints and permission information is available at http://www.nature.com/reprints 78. Cavalcante, R. G. & Sartor, M. A. annotatr: genomic regions in context. Bioinformatics 33, 2381–2383 (2017). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 79. Jobling, M. A. & Tyler-Smith, C. Human Y-chromosome variation in the published maps and institutional affiliations. genome-sequencing era. Nat. Rev. Genet. 18, 485–497 (2017). Open Access This article is licensed under a Creative Commons Acknowledgements Attribution 4.0 International License, which permits use, sharing, M.P.-F. has the support of an INPhINIT Retaining Fellowship from “La Caixa” Foun- adaptation, distribution and reproduction in any medium or format, as long as you give dation (ID 100010434) with code LCF/BQ/DR20/11790032. T.M.-B. is supported by appropriate credit to the original author(s) and the source, provide a link to the Creative funding from the European Research Council (ERC) under the European Union’s Commons license, and indicate if changes were made. The images or other third party Horizon 2020 research and innovation programme (grant agreement No. 864203), material in this article are included in the article’s Creative Commons license, unless PID2021-126004NB-100 (MICIIN/FEDER, UE) and Secretaria d’Universitats i Recerca indicated otherwise in a credit line to the material. If material is not included in the and CERCA Programme del Departament d’Economia i Coneixement de la Generalitat article’s Creative Commons license and your intended use is not permitted by statutory de Catalunya (GRC 2021 SGR 00177). regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ Author contributions licenses/by/4.0/. T.M.-B. and L.F.K.K. conceived the study. O.F., E.J., and E.R. cultured cells and per- formed the flow cytometry. A.F. cultured cells. P.E.-C., C.F., and L.L. performed adaptive © The Author(s) 2023 sampling. I.G. and J.H. performed Illumina sequencing. L.F.K.K., A.S.-A., L.F.-P., and COMMUNICATIONS BIOLOGY | (2023)6:623 | https://doi.org/10.1038/s42003-023-05004-9 | www.nature.com/commsbio 11