ARTICLE https://doi.org/10.1038/s41467-018-07885-5 OPEN Selective single molecule sequencing and assembly of a human Y chromosome of African origin Lukas F.K. Kuderna 1, Esther Lizano 1, Eva Julià2,3, Jessica Gomez-Garrido4, Aitor Serres-Armero1, Martin Kuhlwilm 1, Regina Antoni Alandes4, Marina Alvarez-Estape1, David Juan1, Heath Simon 4,5, Tyler Alioto 4,5, Marta Gut4,5, Ivo Gut4,5, Mikkel Heide Schierup6,7, Oscar Fornas 3,5 & Tomas Marques-Bonet 1,4,5,8,9 1234567890():,; Mammalian Y chromosomes are often neglected from genomic analysis. Due to their inherent assembly difficulties, high repeat content, and large ampliconic regions, only a handful of species have their Y chromosome properly characterized. To date, just a single human reference quality Y chromosome, of European ancestry, is available due to a lack of accessible methodology. To facilitate the assembly of such complicated genomic territory, we developed a novel strategy to sequence native, unamplified flow sorted DNA on a MinION nanopore sequencing device. Our approach yields a highly continuous assembly of the first human Y chromosome of African origin. It constitutes a significant improvement over comparable previous methods, increasing continuity by more than 800%. Sequencing native DNA also allows to take advantage of the nanopore signal data to detect epigenetic mod- ifications in situ. This approach is in theory generalizable to any species simplifying the assembly of extremely large and repetitive genomes. 1 Institut de Biologia Evolutiva, (CSIC-Universitat Pompeu Fabra), PRBB, Doctor Aiguader 88, Barcelona, Catalonia 08003, Spain. 2 Institut Hospital del Mar d’Investigacions Mèdiques (IMIM), Carrer del Doctor Aiguader 88, PRBB Building, Barcelona 08003, Spain. 3 Centre for Genomic Regulation (CRG), The Barcelona Institute for Science and Technology, Carrer del Doctor Aiguader 88, Barcelona 08003, Spain. 4 CNAG‐CRG, Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Baldiri Reixac 4, Barcelona 08028, Spain. 5 Universitat Pompeu Fabra (UPF), Doctor Aiguader 88, Barcelona 08003, Spain. 6 Bioinformatics Research Center, Aarhus University, C.F. Moellers Alle 8, DK-8000 Aarhus C, Denmark. 7 Department of Bioscience, Aarhus University, Ny Munkegade 116, DK-8000 Aarhus C, Denmark. 8 Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluís Companys 23, Barcelona, Catalonia 08010, Spain. 9 Institut Català de Paleontologia Miquel Crusafont, Universitat Autònoma de Barcelona, Edifici ICTA-ICP, c/ Columnes s/n, Cerdanyola del Vallès, Barcelona 08193, Spain. These authors contributed equally: Lukas F.K. Kuderna, Esther Lizano, Oscar Fornas, Tomas Marques-Bonet. Correspondence and requests for materials should be addressed to L.F.K.K. (email: lukas.kuderna@upf.edu) or to E.L. (email: esther.lizano@upf.edu) or to T.M-B. (email: tomas.marques@upf.edu) NATURE COMMUNICATIONS | (2019)10:4 | https://doi.org/10.1038/s41467-018-07885-5 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-018-07885-5 R ecombinational arrest in the common ancestor of the X and to assemble the Y chromosome purely from flow-sorted material Y chromosomes led to the degeneration and accumulation did so using the gorilla8, a species with a previously unchar- of large amounts of repetitive DNA on the Y chromosome acterized Y chromosome, meaning that potential biases in the along its evolutionary trajectory1. Furthermore, many sequencing assembly cannot be detected without a gold standard reference to efforts have traditionally chosen female samples, as the hemi- compare with, such as human. Integrating single-molecule zygous nature of the sex chromosomes leads to half the effective sequencing has been shown to produce far superior whole- sequencing coverage on both of them in a male, resulting in genome shotgun (WGS) assemblies than sequencing by synthesis inferior genome assemblies2,3. Together, these causes have led to platforms11–14. Furthermore, the MinION sequencing platform an underrepresentation of Y chromosomes in genomic studies from Oxford Nanopore Technologies has recently been used to and proper characterization of the Y chromosome in only a create the most contiguous human WGS assembly to date15 and handful of mammalian species through a time- and labor- to resolve the structure of the human Y-chromosome cen- intensive clone by clone approach4–7. One strategy to reduce the tromere16. To take advantage of these benefits, we developed a complexity of the assembly problem for the Y chromosome is to protocol to sequence native, unamplified flow-sorted DNA on the isolate it by flow cytometry, thus dramatically reducing the MinION sequencing device. potential amount of overlaps of repetitive regions in the context of the whole genome3. Notwithstanding previous efforts, which sought to do this have faced some drawbacks, as the material has Results been heavily amplified post sorting to increase yield8. Whole- Flow sorting and sequencing. We sorted approximately 9,000,000 genome amplification (WGA) introduces biases that are detri- individual Y chromosomes from a lymphoblastoid cell line mental to genome assembly, such as unequal sequence coverage (HG02982) from the 1000 Genomes Project, whose haplogroup and chimera formation, as well as limited fragment length9. (A0) represents one of the deepest known splits in humans17 (see Moreover, these methods lead to the loss of epigenetic mod- Fig. 1a). Given the large volume in which the chromosomes were ifications that can now be directly determined from the signal sorted, and potential issues with residual dyes that are necessary for data from nanopore sequencers10. Additionally, previous efforts the sorting process, we devised a purification protocol to bring the a b 60,000 30 Hoechst−33258, 350−460/50 (A.U.) Count 10,000 40,000 Coverage (fold) 20 100 20,000 1 10 0 0 0 20,000 40,000 60,000 ch 1 ch 2 r3 ch 4 ch 5 ch 6 ch 7 ch 8 ch r9 ch 10 ch 11 ch 12 ch 13 ch 14 ch 15 ch 16 ch 17 ch 18 ch 19 ch 20 ch 21 ch 2 ch X rY r r r r r r r r2 r ch ch r r r r r r r r r r r r Chromomycin−A3 457−550/50 (A.U.) c run_1 run_2 d 3000 20,000 2000 Read N50 (bases) 1000 15,000 0 Count run_3 run_4 10,000 3000 2000 5000 1000 0 0 run_1 run_2 run_3 run_4 All 0 0 0 0 0 00 00 0 00 00 0 0 10 10 ,0 ,0 10 10 0, 0, 10 10 10 10 Read length (bases) Fig. 1 Flow-sorting and sequencing specificity. a Flow-karyogram of a human genome. The different clusters correspond to different chromosomes. The red circle delimits the cluster corresponding to the Y chromosome used for this project. b Enrichment specificity of the sequencing data. Sequences on the Y chromosome are ~ 110-fold enriched compared with WGS sequencing. Chromosome 22 partially co-sorts with Y. All other chromosomes are depleted. c Read length (log10 scale) distribution of the four runs. d N50 values for all four runs and the combined dataset. Colors in panels c and d correspond to the different runs. Source data are provided as a Source Data file 2 NATURE COMMUNICATIONS | (2019)10:4 | https://doi.org/10.1038/s41467-018-07885-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-018-07885-5 ARTICLE a b HG02982_chrY GRCh38_chrY GorgorY_ver1 Ampliconic Heterochromatic Others X−degenerate X−transposed Gaps Seg. Dup. (intra) Seg. Dup. (inter) c HG02982 HG38 HG02982 Repeat class DNA 6e+05 LINE LTR RC Basepairs Retroposon 4e+05 RNA rRNA Satellite scRNA 2e+05 SINE snRNA tRNA GRCh38 0e+00 Unknown 0 20 40 60 0 20 40 60 Divergence from consensus Fig. 2 Chromosome-Y assembly overview and comparisons. a Dot-plot comparing the resolved MSY of GRCh38 with HG02982. The reconstruction is highly continuous along most sequence classes, with ampliconic regions showing a higher degree of fragmentation. Seg. Dup. (intra) refers to intra- chromosomal segmental duplications, Seg. Dup. (inter) refers to inter-chromosomal segmental duplications. Altogether, ~ 50% of the of the Y chromosomes resolved sequence space in GRCh38—> 13 Mb—are annotated as segmental duplications. b Treemap comparing the contiguity of HG02982 chrY with GRCh38 chrY and the gorilla Y chromosome by Tomaszkiewicz et al. The size of each rectangle corresponds to the size of a contig within each of the assemblies. Neighboring rectangles are colored differently as a visual aid. c Repeat landscape of common, interspersed repeats annotated equally in GRCh38 and HG02982. Common repeats—including very recent ones—are well resolved in HG02982. The exception are satellite sequences, and a population of somewhat divergent (~ 20%) LTR elements, which are absent in HG02982 (see supplementary Figures 7-9). Source data are provided as a Source Data file DNA into conditions suitable for sequencing. We ran four Oxford Canu18. We performed a self-correction by aligning the reads Nanopore MinION flowcells to generate 305,528 reads summing used for assembly and called consensus using Nanopolish10, to over 2.3 Gb of data. The yields per flowcell varied considerably correcting a total of 127,809 positions. Finally, the Illumina from 897.6 Mb to 163.8 Mb (see Fig. 1c and Supplementary library served to polish residual errors within the assembly using Table 2). Sequencing yields were on the lower end of the reported pilon19. By this means, we corrected a further 101,723 single- spectrum (75 Mb–5.5 Gb per flowcell in15), but read N50 surpassed nucleotide positions and introduced 105,640 small insertions and most of them, ranging from 16.8 to 23.8 kb15 (see Fig. 1d, sup- 6983 small deletions. We also explored further polishing options plementary Figures 1–2). Additionally, for the same flow-sorted and found that running one additional round of error correction material we ran an Illumina MiSeq lane for 2 × 300 cycles but with racon20 potentially resolves several remaining errors, despite including four rounds of PCR amplification. To check the also introducing additional discordances (see Supplementary enrichment specificity, we aligned the reads to the human refer- Table 4, Supplementary Notes 1,3 and Supplementary Figures 10– ence genome (GRCh38) and calculated the normalized coverage 13). The final assembly is comprised of 35 contigs, with an N50 of on each chromosome. Taking into account the size of the Y 1.46 Mb amounting to 21.5 Mb of total sequence, in contrast to a chromosome and its haploid nature, we find it to be over 110-fold contig N50 of 6.91 Mb of the GRCh38 Y-chromosome assembly. enriched compared with a random sampling from the human Compared with the gorilla Y-chromosome assembly with a contig genome (see Fig. 1b, Supplementary Figures 3–5, Supplementary N50 of 17.95 kb8, our assembly is two orders of magnitude more Note 2 and Supplementary Data 1–2). contiguous (see Fig. 2b). The Y chromosome is comprised of a set of discrete sequence classes4. To check the completeness of our assembly, we assessed Y-chromosome assembly and comparison with GRCh38. We how well each of them is represented. After retaining only single used the Nanopore data to construct a de novo assembly using best placements, we were able to align 21.1 Mb, or 98.4% of its NATURE COMMUNICATIONS | (2019)10:4 | https://doi.org/10.1038/s41467-018-07885-5 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-018-07885-5 Table 1 Assembly statistics overview Seq. class Aln. HG02982 HG02982 ID HG02982 ID SNP Rec. HG02982 Aln. NA24385 Rec. in Len. w/o gaps (b) SNP (%) + InDel(%) (%) (b) NA24385 (%) (b) Ampliconic 6,146,087 99.91 99.67 62.67 5,242,461 53.46 9,807,089 Heterochromatic 543,005 99.66 99.31 32.77 171,045 10.32 1,656,797 Others 295,160 99.47 99.18 385.59 63,973 83.57 76,547 Pseudo-autosomal 2,219,743 99.58 99.13 78.02 117,626 4.13 2,844,939 X-degenerate 8,537,493 99.95 99.81 98.94 8,238,733 95.48 8,628,904 X-transposed 3,374,011 99.94 99.81 99.21 1,474,610 43.36 3,400,750 Summary of sequence class coverage of HG02982 versus GRCh38, as well as the contigs from NA23385 identified as derived from the Y chromosome. The proportion of recovered sequences and % identity are calculated over the resolved sequences in GRCh38, excluding gaps. There are currently 30.8 Mb of unresolved sequence (represented by the ambiguous base N) in the reference Y chromosome of GRCh38, the vast majority of which belongs to heterochromatin on the q arm Aln.: aligned bases to GRCh38, ID.: percent identical bases in GRCh38, Rec.: recovered proportion from GRCh38, Len.: length in GRCh38 total length, with 99.9% of identical bases on average (see Fig. 2a). male-specific region of the Y chromosome (MSY) region and are We recovered the full-length (~ 99% of the annotated length in therefore able to annotate all single-copy genes. Furthermore, we GRCh38) reconstructions of both the X-transposed and the X- are able to retrieve at least one member of all multi-copy gene degenerate regions. Although the X-degenerate region can be families. For four out of nine of these gene families, we are considered a single-copy region due its distant common ancestry additionally able to resolve further copies within our assembly with the X chromosome, the X-transposed region emerged only (see Supplementary Data 4–5). We also note that four genes after the split between humans and chimpanzees21. The largest (ASMTL, IL3R, P2RY, SLC25) from a comparatively short syntenic sequence class on the Y chromosome is comprised of ampliconic block of around 200 kb are partially missing from our assembly regions, which amount to around 30% of the euchromatic portion due to the aforementioned technical challenges in the PAR-1 and sum to 9.93 Mb. These regions contain eight massive, region. Mapping the raw data onto GRCh38 show that this is an segmentally duplicated palindromes, all of which share >99.9% artifact, presumably due to insufficient coverage in this region. identity between their two copies, with the largest one spanning over 2.90 Mb. We find this region to be the most challenging to Structural variants. We produced a stringent call set of structural reconstruct, with fragmented and collapsed sequences, but are variants (SVs) derived from alignments to GRCh38 using nevertheless able to recover 6.14 Mb, or 62.7% of its length in Assemblytics25. We detect 347 SVs at least 50 bp in size (931 GRCh38. Surprisingly, we recover only 78% of the pseudo- variants at least 10 bp in size) of which 82 are at least 500-bp long autosomal regions (PARs). We observed a rather steep drop-off in (see Fig. 2c, Supplementary Figures 14–15, Supplementary coverage coinciding with the PAR-1 boundary on GRCh38. As we Table 3). The cumulative length of these variants sums to 184 kb. are sequencing native, unamplified DNA, the genomic coverage is We observe a 4.8-fold excess number of deletions versus number directly proportional to the number of copies of the underlying of insertions, amounting to a twofold excess of bases in deletions sequenced region22. We compared the mapped coverage of our versus bases in insertions. Although a deletion bias for nanopore- raw data on GRCh38 and find that PAR-1 exhibits only around based assemblies had previously been reported15, we find the 72% of the average coverage of the whole chromosome (19.8-fold strength of this bias to be decreasing in our analysis, probably versus 27.3-fold). We observe the drop-off in coverage to coincide reflecting improvements in base-calling accuracy. To check the sharply with the PAR-1 boundary (see Supplementary Figure 6). presence of large-scale copy number variation in multi-copy Finally, of the remaining sequence classes, we are able to recover genes, we additionally determined the chromosome-wide copy around 32.8% of the resolved heterochromatic regions, and number based on a read depth approach using the Illumina data. multiple instances of the remaining unclassified sequences We find extensive genic copy number variation, with expansions (referred to as other; see Table 1 and Supplementary Data 3). in five of the nine multi-copy genes, when compared with the To contrast our approach to a long-read WGS assembly, we reference individual. Among these, we find expansions in RBMY, assembled the publicly available PacBio dataset from the PRY, BPY2, and DAZ, all members of the AZFc region locus with Ashkenazim son from the Genome in a Bottle Consortium23, implications for male fertility. Although these expansions are to which has a sequencing depth comparable to ours on the sex some degree represented in our assembly, the precise genomic chromosomes (~ 30X). We identified 193 contigs mapping to the architecture remains challenging to reconstruct. Due to the high Y chromosome, with an N50 of 213 kb, covering 15.3 Mb, or degree of similarity between copies, several of them will be col- around 28% less than by our approach. The WGS fails to lapsed in the assembly specially in the AZFc region. Finally, to assemble roughly 56.6% of the X-transposed region and 47% of assess concordance with previous studies, we compared our SV the ampliconic regions (see Table 1). calls with those generated by the 1000 Genomes Project, which contains the same cell line used for this study26. We manually Comparative gene annotation. We performed a comparative confirm the presence of all structural three variants called in annotation to check the completeness of our assembly at the gene HG02982 in the 1000 Genomes Project in our data by checking level. To this end, we projected all Gencode (v. 27, GRCh38) the overlap of calls produced by orthogonal approaches (see annotations on the Y chromosome onto our assembly and Supplementary Figures 16–18). annotated them there. Due to its peculiar evolutionary trajectory, the gene-space on the Y chromosome is degenerated, and any CpG methylation status. Finally, we called the methylation status remaining genes can generally be classified into two categories: on of 5-methylcytosines (5-mC) at CpG positions from the Nano- one hand there are single-copy genes, which are broadly expressed pore signal data using a recently developed model implemented beyond the testis. On the other, there are multi-copy genes within in Nanopolish10. To assess potential biases on the CpG methy- the ampliconic regions, which are mainly involved in spermato- lation status introduced by our workflow, we also produced genesis24. We recover the complete gene set of the genes in the whole-genome bisulfite sequencing data (WGBS) for the same cell 4 NATURE COMMUNICATIONS | (2019)10:4 | https://doi.org/10.1038/s41467-018-07885-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-018-07885-5 ARTICLE line. We calculated the methylation frequency (i.e., the proportion complexity reduction. Given different chromosomes that are of reads supporting 5-mc at a given CpG) for both datasets. For sufficiently distinguishable in terms of size and GC content, positions where both datasets have at least 10-fold coverage (n = immediate applications are either very complex chromosomes, 4654), we observe a good concordance in the methylation fre- such as the human Y, or extremely large genomes with a very quency with a Pearson’s r of 0.816 (see Supplementary Figure 19 high degree of common repeats, which have long challenged and Supplementary Table 5). Remaining differences might be traditional WGS approaches, such as wheat, the loblolly pine, or attributable to differences in sensitivity, variation in the methy- the axolotl27–30. lation state, or alternative modifications such as 5-hydro- xymethylation, which cannot be distinguished from 5-mC by Methods WGBS10. Additionally, detecting the 5-mC status on the Y Chromosome preparation for flow karyotyping. Mitotic chromosomes in sus- chromosome from long reads in our methodology has the pension were prepared as follows (adapted from31 with some modifications): the advantage of allowing to interrogate regions that are not acces- lymphoblastoid cell line HG02982 (purchased from Coriell, cat. no. HG02982) were cultured in RPMI-1640 medium supplemented with 2mM L-glutamine sible to WGBS with short reads, namely the PAR, the X- (Invitrogen, ref. 21875-034), 15% fetal bovine serum and antibiotics (penicillin and transposed region, and to some degree the Ampliconic regions. streptomycin (Invitrogen, ref. 15140-122)) at initial concentration no <150,000 We interrogated the methylation state of CpG 200-bp upstream viable cells per ml. Near confluence, cells were subcultured to 50%. After 24 h, of the transcription start site (TSS) in protein-coding genes falling the cells were blocked in mitosis by adding Colcemid to the culture (10 µg ml−1 demecolcine solution (Gibco, ref. 15210-040)) to a final concentration of 0.1 µg ml−1 within the different sequence classes of the Y chromosome. Genes and incubated for an additional 6–7 h. To swell and stabilize mitotic cells, they were from the PAR, X-degenerate, and X-transposed regions are centrifuged 5 min at 300 × g at room temperature. The pellet was slowly resuspended expressed throughout the body, whereas Ampliconic genes have in 10 ml hypotonic solution (Hypotonic solution: 75 mM KCl, 10 mM MgSO4, testis specific expression24. In agreement with these patterns, we 0.2 mM spermine, 0.5 mM spermidine. pH 8.0), incubated for 10 min at room temperature. After the incubation in the hypotonic solution, the swollen cells were find the genes within PAR, X-degenerate, and X-transposed centrifuged at 300 × g for 5 min. The cell pellet was resuspended in 1.5 ml of ice-cold regions to show low degrees of CpG methylation at TSS. Within polyamine isolation buffer (PAB: 15 mM Tris, 2 mM EDTA, 0.5 mM EGTA, 80 mM the Ampliconic regions, the distribution of methylation fre- KCl, 3 mM dithiothreitol, 0.25% Triton X-100, 0.2 mM spermine, 0.5 mM spermi- quencies of CpGs at TSS shows an overall high degree of dine. pH 8.0) for 20 min to release the chromosomes. To ensure the integrity of the chromosomes, their morphology was checked methylation and is therefore consistent with the expected before staining them. To this end, the pellet was vigorously vortexed for 30 s to downregulation of these genes in lymphoblastoid cells (see Sup- liberate the chromosomes from the mitotic cells. The suspension was filtered plementary Figures 20–21). Nevertheless, single-copy resolution through a 35 µm mesh filter and stored at 4 °C until its sorting. is not possible due to potential mapping ambiguities. Finally, chromosomes were stained with chromomycin-A3 (Sigma, ref. C2659) and Hoechst 33,258 (Invitrogen, ref. H3569) at a final concentration of 40 µg ml−1 and 5 µg ml−1, respectively, in presence of divalent cations (10 mM MgSO4 (Sigma, ref. 60142)). Staining was performed for at least 8 h at 4 °C, to allow the dyes to Discussion equilibrate. Before the sample analysis on a cell sorter, potassium citrate was added Here, we report the first successful sequencing and assembly of to a final concentration of 10 mM (Sigma, ref. 89306) to enhance peak resolution in native, flow-sorted DNA on an Oxford Nanopore sequencing the flow karyotype. device, without previous amplification. We apply our methodol- ogy to assemble the first human Y chromosome of African origin Chromosome sorting. Flow karyotyping for chromosome sorting was performed on BD Influx cell sorter (Becton Dickinson, San Jose, CA), a jet-in-air cell sorter to benchmark our approach. This is arguably the most challen- that was selected for its relatively easy manual daily fine-tuning and high-resolution ging human chromosome to assemble due to its high repeat and capabilities. Of the five available lasers, only the blue (488 nm laser at 200 mW), segmental duplication content, and hence a good test-case to deep-blue (457 nm laser at 300 mW), and ultraviolet (355 nm laser at 100 mW) explore the possibilities and limitations of this approach. With the ones were used for flow karyotyping. The setup and performance were optimized exception of bacterial artificial chromosome-based assemblies, we using standard 8-peaks Rainbow beads (SpheroTM Rainbow Calibration Particles 3.0–3.4 μm, BD Biosciences, ref. 559123), 1-peak UV beads for UV laser alignment are able to reconstruct the Y chromosome to unprecedented (AlignflowTM Flow Cytometry Alignment 2.7 μm, Molecular Probes, ref. A16502), quality in terms of contiguity and sequence class representation. and 1-peak 457 nm for deep-blue laser alignment (FluoresbriteTM Plain YG We show that we not only outperform previous efforts that Microspheres 1.0 μm, Polysciences, Inc. ref. 17154) were, respectively, used for 488- sought to achieve a similar goal of reconstructing Y chromo- blue, 355-UV, and 457-deep-blue optimal laser alignment and instrument fine tuning to obtain the highest resolution of chromosome detection and sorting. somes8, but also accomplish a better reconstruction on all The threshold for chromosome sorting was set triggering in chromomycin-A3 sequence classes than the Y chromosomal sequences derived from fluorescence on 457 nm laser as primary excitation line and set at approximately a long-read WGS assembly. Additionally, our method is orders of 1800 a.u. Then, chromomycin-A3 fluorescence was used as primary fluorescence magnitude cheaper than reconstructions from WGS data too, reference through a light line of 500 LP filter and collected by a 550/50 nm band- pass filter. Hoechst was excited with the UV laser and its fluorescence was collected especially considering that twice the desired Y chromosomal through a light line of 400 LP filter and by 460/50 BP. All parameters were collected target coverage is needed on the autosomes. Given the current in lineal mode and analyzed with the BD FACSTM Software (v. 1.0.0.0.650, Becton developments in sequencing throughput, a single-MinION flow- Dickinson, San Jose, CA). cell should now be sufficient to assemble a whole human Y We chose a 100 μm nozzle because we found it to have the best piezoelectric frequency/electronic-noise ratio. The piezoelectric frequency was adjusted at chromosome. Furthermore, it is becoming clear that the upper 38.7 KHz. The sample flow rate for chromosome sorting was adjusted at up to 6000 read length boundary is only delimited by the integrity of the events s−1. The gating strategy for chromosome sorting was simple because only a DNA, suggesting the possibility that complete Y-chromosome bi-parametrical dot-plot Hoechst versus chromomycin-A3 fluorescence was used assemblies, including full resolution of amplicons, might be (see Fig. 1a). possible in the near future. Notwithstanding, some challenges to obtain ultra-long reads from flow-sorted chromosomes are still to Purification and concentration of flow-sorted Y chromosomes. For each of the be overcome, as sorting sufficient material for this protocol is a two rounds of purification, the fractions corresponding to approximately 4.5 M Y chromosomes (~ 500 ng of DNA per aliquot) were divided into 1 ml aliquots with substantial endeavor. It also is worth noting that our efforts to an estimated chromosome count of 400,000, corresponding to a DNA con- sequence the same input material on Pacific Biosciences Sequel centration of approximately 0.04 ng µl−1. The approximate total volume per round platform have been fruitless, presumably due to interference of of purification was around 22.5 ml. Each tube containing the flow-sorted DNA was residual dyes with the sequencers optical detection system. treated overnight with 10 µl of proteinase K (20 mg ml−1) at 50 °C. After treatment, the buffer was exchanged, and proteinase K, as well as chromomycin-A3 and Despite the technical challenges of flow-sorting single chromo- Hoechst 33,258 removed by dialysis against 1 liter of TE buffer using a Pur-A- somes, the method described here offers the opportunity to take Lyzer™ Maxi Dialysis column with a molecular weight cut-off of 50 kDa (Sigma- advantage of the benefits of long-range data together with local Aldrich). Dialysis was carried out for 48 h exchanging the buffer every 10–16 h. To NATURE COMMUNICATIONS | (2019)10:4 | https://doi.org/10.1038/s41467-018-07885-5 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-018-07885-5 reduce the volume after buffer exchange, DNA was transferred into 1.5 ml tubes duplicates were removed with Picardtools (v. 2.8.2, https://broadinstitute.github.io/ and concentrated by evaporation in a miVac DNA concentrator (Barnstead Gen- picard/). eVac, Ipswich, UK) up to a volume of approximately 5–10 µl. A final purification bwa mem HG02982_canu.selfcorrected.fasta step was performed by pooling the concentrated DNA into two tubes and sub- reads.p1.fastq reads.p2.fastq | samtools sort jecting it to a solid-phase reversible immobilisation (SPRI) bead purification with a -o reads.paired.mappings.bam -T tmp -; 2X ratio (SPRI beads/sample). DNA was eluted in 9 µl of low TE buffer and pooled java -jar picard.jar MarkDuplicates I=reads. into one tube. Concentrations were determined by absorbance at 260 nm with a paired.mappings.bam O=reads.paired.mappings. NanoDrop 2000 (Thermo Scientific) and by fluorometric assay with the Qubit 2.0 markdup.bam M=reads.paired.mappings.markdup. using the Qubit dsDNA HS kit (Invitrogen) (see Supplementary Table 1). bam Polishing was performed with Pilon (v 1.22)19, resulting in 132,336 residual errors being corrected. Sequencing the flow-sorted chromosomes. The purified DNA was prepared for java -Xmx96G -jar pilon-1.22.jar --threads 12 sequencing following the protocol in the Rapid Sequencing kit SQK-RAD002 --genome HG02982_canu.selfcorrected.fasta (ONT, Oxford, UK). Briefly, approximately 200 ng of purified DNA was tagmented --frags reads.paired.mappings.markdup.bam for 1 min at 75 °C with the Fragmentation Mix (ONT, Oxford, UK). The Rapid --output HG02982_canu.selfcorrected.pileon Adapters (ONT, Oxford, UK) were added along with Blunt/TA Ligase Master Mix --outdir pilon_corrections --changes --vcf (NEB, Beverly, MA) and incubated for 30 min at room temperature. The resulting --tracks --fix all library was combined with Running Buffer with Fuel (ONT, Oxford, UK) and To run racon (v 1.3.1, see Supplementary Table 4), we mapped the Illumina Library Loading Beads (ONT, Oxford, UK) and loaded onto a primed R9.4 Spot- reads onto the polished reference with bwa, sorted the alignments with samtools On Flow cell (FLO-MIN106). Sequencing and initial base calling was performed and removed duplicates as described above. The resulting alignments were with a MinION Mk1B MinKNOW v1.7.10 software package running for 48 h. provided to racon as an input: Estimates for DNA quantification were based on chromosomal counts with cor- racon -u -t 12 reads.fastq mappings.sam responding quantification values from Gribble et al.31. The uncertainties in HG02982_chrY_v1.fasta quantification with Qubit 2.0 or NanoDrop are presumed to be due to residual intercalating dyes present within the sample, which interfere with the quantifica- tion platforms detection systems, with competition of additional intercalants Variant calls. For variant calls, the Illumina data were mapped onto the GRCh38 leading to underestimation on the Qubit 2.0, and the additional presence of aro- or the HG02982 assembly, respectively, and processed the same way as detailed matic groups leading to overestimation on the NanoDrop. above. Variants were called using GATKs Haplotype Caller with the following A total estimated amount of 100 ng of Y chromosome was fragmented on a optional flags: “--genotyping-mode DISCOVERY --sample-ploidy 1”. Covaris ultrasonicator with settings targeting fragments of 450 bp. The library was java -jar gatk-package-4.0.0.0-local.jar prepared using NEBNext Ultra II DNA Library Prep Kit (New England BioLabs) HaplotypeCaller -R reference.fa -I mappings. following the manufacturer’s instructions, including four cycles of PCR bam --genotyping-mode DISCOVERY -O variants. amplification. Agilent BioAnalyzer High-Sensitivity DNA Kit was used to determine vcf the size distribution and molarity. The library was sequenced on an Illumina MiSeq using the v3 kit and 600 cycles resulting in 300-bp paired-end reads. Repeat annotations. Repeat annotations were performed using RepeatMasker (v. 4.0.7) with rmblastn v. 2.6.0+ as the engine. To be comparable, the annotations for both the HG02982, as well as the GRCh38 assembly were performed the same way. Assembly, error correction, and polishing. The initial base-calls (MinKNOW We used the RepBase-20170127 as the repeatmasker database, and Homo sapiens 1.7.10 using Albacore 1.1) from the Nanopore data were assembled with Canu (v as the query species. Divergence of the repeat annotations to their consensus was 1.6)18 without previous read separation of reads deriving from different chromo- calculated using the “calcDivergenceFromAlign.pl” utility included in the Repeat- somes and assuming a chromosome size of 52 Mb. The following parameters were Masker package. used: RepeatMasker -e ncbi -pa 12 -s -species human canu -p HG02982 -d HG02982_canu genomeSize=52 -no_is -noisy -dir ./outDir -a -gff -u reference. m overlapper=mhap utgReAlign=true -nanopore- fa raw raw_data/HG02982/all.joint.fastq The 2.3 Gb of input data resulted in 25X of error corrected reads for assembly, assuming a chromosome size of 52 Mb. The data assembled into 35 contigs, which Whole-genome alignments. Whole-genome alignments to GRCh38 were pro- where self-corrected using the Nanopore input reads. To this end, we re-performed duced using last (v. 914) with the following parameters as suggested by the base calling from the fast5 files using Albacore (v 2.1, available from the nanopore developer for highly similar genomes for indexing and alignments: user community) to be used for variant calling with Nanopolish (v. 0.8.4, https:// lastdb -uNEAR -R01 index reference.fa github.com/jts/nanopolish, 11 December 2017). lastal -e25 -v -q3 -j4 index query.fa > read_fast5_basecaller.py -f FLO-MIN106 -k mappings.maf SQK-RAD002 -i input_folder -s outout_folder -t 8 Single best placements of query sequences were retained using the “last-split” -o fastq,fast5 -q 10000000 -n 100000 --disable_ script included in the last alignment package. Alignments were filtered for a pings maximum mismap probability of 10e–5. The alignments were converted to psl We indexed the reads to be used with Nanopolish: format for further processing. nanopolish index -f fast5.fofn reads.joint. fastq Comparison with WGS PacBio data. The PacBio data from the Ashkenazim Son The reads were mapped onto the raw assembly using bwa mem (v. 0.7.120)32 (Coriel ID NA24385) produced by the genome in a bottle consortium was also with the additional flag -x ont2d and the mappings merged and sorted with assembled using Canu (v. 1.6) using default assembly parameters and assuming a samtools (v. 1.5): genome size of 3.2 Gb: bwa mem -x ont2d HG02982_canu.uncorrected. canu -p NA24385 -d NA24385_canu genomeSize=3. fasta reads.joint.fastq | samtools sort -o 2g -pacbio-raw data/fastq/*fastq.gz reads.joint.mappings.bam -T tmp - gridOptionsExecutive=’--mem-per-cpu=16g The mappings were fed to Nanopolish and corrected in chunks of 50 kb using --cpus-per-task=2’ the helper script “nanopolish_makerange.py” included in the Nanopolish package. After genome assembly, we performed a whole-genome alignment to GRChg38 Variants were called using “nanopolish variants –consensus” with the optional flag and retained single best placements as mentioned above. To identify contigs “--min-candidate-frequency 0.1”. belonging to the Y chromosome, we performed the following filtering steps: for nanopolish_makerange.py HG02982_canu. contigs, which have local best placements on a chromosome different than the Y, uncorrected.fasta | xargs -i echo nanopolish we filtered out those whose proportion of mapped bases is higher on a sequence variants --consensus selfcorrected.{}.fa -w {} from the reference assembly different from the Y chromosome. Additionally, we -r reads.joint.fastq -b reads.joint.mappings. filtered out any alignments with a mismap probability higher than 10e–5. By this bam -g HG02982_canu.uncorrected.fasta -t 4 means, we retained 184 contigs mapping 15,308,468 base pairs on the Y --min-candidate-frequency 0.1 | sh chromosome (see Supplementary Data 6) By this means, we corrected 127,801 positions in the initial assembly. The self- corrected assembly was further polished with the Illumina library. To this end, we trimmed the Illumina reads to get rid of any adapters in the sequences using SV calls. SVs were called with assemblytics25. To this end, we produced whole- trimgalore (v 3.7, https://github.com/FelixKrueger/TrimGalore). genome alignments using nucmer from the Mummer package (v. 3.22)33. The trim_galore --fastqc --paired --retain_ resulting delta file was passed to assemblytics, with the required unique anchor unpaired gzip pair1.fastq pair2.fastq length set to 10000 bp. The trimmed reads were mapped with BWA mem (v.0.7.12)32 in paired-end nucmer -maxmatch -l 100 -c 500 GRCh38.chrY.fa mode and the mappings converted to a sorted bam files using samtools sort. PCR HG02982_chrY_v1.fasta -prefix HG02982_vs_HG38 6 NATURE COMMUNICATIONS | (2019)10:4 | https://doi.org/10.1038/s41467-018-07885-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-018-07885-5 ARTICLE Assemblytics HG02982_vs_HG38_.delta nanopolish call-methylation -v --progress -t HG02982_chrY_v1.vs.hg38_10kanchor.50kmax 8 -r joint_reads.fastq -b joint_reads. 10000 bin/Assemblytics/ mappings.bam -g chrY.fa -w"chrY:$start-$stop" > methylation_calls.tsv Finally, we calculated the methylation frequency and log-likelihood ratios of Read depth duplication detection. We estimated absolute copy number with a methylation at each position: depth of coverage approach using the Illumina data22. We masked all common calculate_methylation_frequency.py -i repeats as identified by RepeatMasker (see above) and tandem repeat finder. We methylation_calls.tsv created non-overlapping 36-mers of the raw reads, which were mapped onto the We filtered out any position with <10 reads in either the WGBS or the assembly using GEM (v 2)34 allowing for a divergence of up to 5%. The read depth Nanopore data. Additionally, any position with a log-likelihood ratio of <2.5 in the was calculated in non-overlapping windows of 1 kb of non-repetitive sequence. Nanopore data were also excluded. After correcting for GC content using mrCanavar (v. 0.51), we normalized by the mean read depth. To assign a copy number to each gene, we calculated the median copy number of all windows intersecting a gene. For the hg38 Y chromosome, a set Code availability. The custom script used for the gene annotation has been of custom single-copy regions needed to be provided to the CN caller as calibra- deposited at Figshare at https://doi.org/10.6084/m9.figshare.7359065.v1. tion. These regions were inferred by subtracting the reference WGAC (whole- genome assembly comparison, UCSC track genomic superdups) segmental dupli- cation track from the whole Y chromosome and keeping only stretches of single- Data availability All raw sequencing data for this study have been deposited at the European copy sequence longer than 2 kb. Nucleotide Archive (ENA) under the study accession PRJEB28143. The assembly is deposited at the ENA under the accession ULGL01000000. The WGS assembly for Gene annotation. The annotation of the HG02982 assembly was performed by the NA24385 individual, the repeat masker tracks for the GRCh38 chrY and trying to assign the genes present in the Y-chromosome annotation of GRCh38 HG02982 assemblies, and the methylation calls from the Illumina WGBS gencode version 27. For this purpose, we downloaded the gff3, the transcript and the Nanopore data are deposited at Figshare at https://doi.org/10.6084/m9. sequences and the protein sequences that corresponded to the Y-chromosome figshare.7358480.v1. The source data underlying Figs. 1a-d and 2a-c are provided annotation and performed transcript and protein mappings with GMAP (v. as a Source Data File. A Reporting Summary for this Article is available as a 2017031735) and exonerate (v. 2.2.036), respectively. Additionally, a numeric index Supplementary Information file. All other relevant data are available upon request. was assigned to each gene in the HG38 Y chromosome according to the order in the chromosome. Next, we combined all the data (transcript mappings, protein mappings and gene synteny) with an in-house script (available at https://doi.org/ Received: 29 June 2018 Accepted: 2 December 2018 10.6084/m9.figshare.7359065.v1) to locate each gene in our assembly and assign parts of the assembly to their corresponding region in the Y chromosome of GRCh38. After following the strategy mentioned above for all the genes, we took a closer look to the protein-coding genes, by manually checking some of the map- pings in order to determine possible errors in the sequence caused by the Nanopore References reads that could introduce frameshifts or internal stop codons in the aminoacidic 1. Charlesworth, B. & Charlesworth, D. The degeneration of Y chromosomes. sequence. Philos. Trans. R. Soc. B Biol. Sci. 355, 1563–1572 (2000). 2. Hughes, J. F. & Page, D. C. The biology and evolution of mammalian Y chromosomes. Annu. Rev. Genet. 49, 507–527 (2015). Illumina WGBS sequencing and methylation calls. Two micrograms of genomic 3. Tomaszkiewicz, M., Medvedev, P. & Makova, K. D. Y and W chromosome DNA from a lymphoblastoid cell line (HG02982) were spiked with unmethylated assemblies: approaches and discoveries. Trends Genet. 33, 266–282 (2017). bacteriophage λ DNA (5 ng of λ DNA per microgram of genomic DNA; Promega) 4. Skaletsky, H. et al. The male-specific region of the human Y chromosome is a and with methylated T7 phage DNA (5 ng of T7 DNA per microgram of genomic mosic of discrete sequence classes. Nature 423, 825–837 (2003). DNA). The DNA was sheared to 50–500 bp in size using Covaris LE220 ultra- 5. Hughes, J. F. et al. Chimpanzee and human y chromosomes are remarkably sonicator, and fragments of 150–300 bp were size-selected using AMPure XP beads divergent in structure and gene content. Nature 463, 536–539 (2010). (Agencourt Bioscience). The libraries were constructed using the KAPA Library 6. Hughes, J. F. et al. Strict evolutionary conservation followed rapid gene loss on Preparation Kit with no PCR Library Amplification/Illumina series (Roche-Kapa human and rhesus y chromosomes. Nature 483, 82–87 (2012). Biosystems) together with the NEXTFLEX® Bisulfite-Seq Barcodes (Perkin Elmer). After adaptor ligation, the DNA was treated with sodium bisulfite using the Epi- 7. Soh, Y. Q. S. et al. Sequencing the mouse y chromosome reveals convergent Tect Bisulfite kit (Qiagen) following the manufacturer's instructions. Enrichment gene acquisition and amplification on both sex chromosomes. Cell 159, for adaptor-ligated DNA was carried out through seven PCR cycles using KAPA 800–813 (2014). HiFi HotStart Uracil+ReadyMix PCR 2x Kit (Roche-Kapa Biosystems). Library 8. Tomaszkiewicz, M. et al. A time- and cost-effective strategy to sequence quality was monitored using the Agilent 2100 Bioanalyzer DNA 7500 assay, and mammalian Y chromosomes: an application to the de novo assembly of gorilla the library concentration was estimated using quantitative PCR using the KAPA Y. Genome Res. 26, 530–540 (2016). Library Quantification Kit for Illumina® Platforms, v1.14 (Roche-Kapa 9. Zhang, K. et al. Sequencing genomes from single cells by polymerase cloning. Biosystems). Nat. Biotechnol. 24, 680–686 (2006). Paired-end DNA sequencing (2×101 bp) of the converted libraries was 10. Simpson, J. T. et al. Detecting DNA cytosine methylation using nanopore performed using the HiSeq 2500 (Illumina) following the manufacturer’s protocol sequencing. Nat. Methods 14, 407–410 (2017). with HiSeq Control Software (HCS) 2.2.68. Primary data analysis, image analysis, 11. Berlin, K. et al. Assembling large genomes with single-molecule sequencing base calling, and quality scoring of the run, was processed using the manufacturer’s and locality-sensitive hashing. Nat. Biotechnol. 33, 623–630 (2015). software Real Time Analysis (RTA 1.18.66.3) and followed by generation of 12. Chaisson, M. J. P. et al. Resolving the complexity of the human genome using FASTQ sequence files by CASAVA. single-molecule sequencing. Nature 517, 608–611 (2015). We used the gemBS pipeline37 using the default parameters to perform the 13. Kuderna, L. F. K. et al. A 3-way hybrid approach to generate a new high- analysis. The reference genome used for the alignment was GRCh38. Methylated quality chimpanzee reference genome (Pan tro 3.0). Gigascience 6, gix098 and unmethylated cytosine conversion rates were determined from spiked-in (2017). bacteriophage DNA (fully methylated phage T7 and unmethylated phage lambda). 14. Bickhart, D. M. et al. Single-molecule sequencing and chromatin The under and over conversion rates for the sample were <1 and ~ 5%, respectively. conformation capture enable de novo reference assembly of the domestic goat Only uniquely mapping reads were retained for downstream analysis. The genome. Nat. Genet. 49, 643–650 (2017). comparison with the Nanopore calls was performed for all canonical CpG sites on 15. Jain, M. et al. Nanopore sequencing and assembly of a human genome with the Y chromosome where there was sequencing data available from both ultra-long reads. Nat. Biotechnol. 36, 338–345 (2018). experiments. The comparison took account of the variable precision of the 16. Jain, M. et al. Linear assembly of a human centromere on the Y chromosome. methylation estimates due to variation in sequencing coverage between sites so that Nat. Biotechnol. 36, 321–323 (2018). low-coverage sites did not affect the comparison. 17. Auton, A. et al. A global reference for human genetic variation. Nature 526, 68–74 (2015). Nanopore methylation calls. The methylation status was called using Nanopol- 18. Koren, S. et al. Canu : scalable and accurate long-read assembly via adaptive k- ish10 as suggested by the developers. To this end, we aligned the Nanopore reads to mer weighting and repeat separation. Genome Res. 27, 722–736 (2016). GRCh38 with minimap238 and sorted with samtools (v 1.5). The calls were per- 19. Walker, B. J. et al. Pilon: an integrated tool for comprehensive microbial formed in 200 kb windows. variant detection and genome assembly improvement. PLoS ONE https://doi. minimap2 -a -x map-ont chrY.fa joint_reads. org/10.1371/journal.pone.0112963 (2014). fastq | samtools sort -T tmp -o joint_reads. 20. Vaser, R., Sović, I., Nagarajan, N. & Šikić, M. Fast and accurate de novo mappings.bam genome assembly from long uncorrected reads. Genome Res. 27, 737–746 samtools index joint_reads.mappings.bam (2017). NATURE COMMUNICATIONS | (2019)10:4 | https://doi.org/10.1038/s41467-018-07885-5 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-018-07885-5 21. Page, D. C., Harper, M. E., Love, J. & Botstein, D. Occurrence of a Instituto de Salud Carlos III, from the Generalitat de Catalunya through the Departa- transposition from the X-chromosome long arm to the Y-chromosome short ment de Salut and Departament d’Empresa i Coneixement, and co-financing by the arm during human evolution. Nature 311, 119–123 (1984). Spanish Ministry of Economy, Industry and Competitiveness (MEIC) with funds from 22. Alkan, C. et al. Personalized copy number and segmental duplication maps the European Regional Development Fund (ERDF) corresponding to the 2014–2020 using next-generation sequencing. Nat. Genet. 41, 1061–1067 (2009). Smart Growth Operating Program. L.F.K.K. is supported by an FPI fellowship associated 23. Zook, J. M.et al. Extensive sequencing of seven human genomes to with BFU2014-55090-P (MINECO/FEDER, UE). M.K. is supported by a Deutsche characterize benchmark reference materials. Sci. Data. 3, 160025 (2016). Forschungsgemeinschaft (DFG) fellowship (KU 3467/1-1). T.M.-B. is supported by 24. Hughes, J. F. & Rozen, S. Genomics and genetics of human and primate Y BFU2017-86471-P (MINECO/FEDER, UE), U01 MH106874 grant, Howard Hughes chromosomes. Annu. Rev. Genom. Hum. Genet. 13, 83–108 (2012). International Early Career, Obra Social “La Caixa” and Secretaria d’Universitats i Recerca 25. Nattestad, M. & Schatz, M. C. Assemblytics: a web analytics tool for the del Departament d’Economia i Coneixement de la Generalitat de Catalunya. D.J. is detection of variants from an assembly. Bioinformatics 32, 3021–3023 supported by a Juan de la Cierva fellowship (FJCI-2016-29558) from MICINN. (2016). 26. Poznik, G. D. et al. Punctuated bursts in human male demography inferred Author contributions from 1,244 worldwide Y-chromosome sequences. Nat. Genet. 48, 593–599 T.M.-B. conceived the study, L.F.K.K., J.G.-G., A.S.A., M.K., S.H., D.J., and T.A. per- (2016). formed computational analysis, E.L. developed and performed the purification protocol, 27. Lukaszewski, A. J. et al. A chromosome-based draft sequence of the hexaploid E.J. and O.F. cultured cells and performed the flow cytometry, R.A.A., M.A.-E., M.G., I. bread wheat (Triticum aestivum) genome. Science 345, 6194 (2014). G., and M.H.S. prepared materials and/or performed the sequencing. L.F.K.K., E.L., and 28. Zimin, A. V. et al. The first near-complete assembly of the hexaploid bread T.M.-B. wrote the manuscript with input from all authors. All authors approved the wheat genome, Triticum aestivum. Gigascience 6, gix097 (2017). manuscript. 29. Neale, D. B.et al. Decoding the massive genome of loblolly pine using haploid DNA and novel assembly strategies. Genome Biol. 15, R59 (2014). 30. Nowoshilow, S. et al. The axolotl genome and the evolution of key tissue Additional information formation regulators. Nature 554, 50–55 (2018). Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- 31. Gribble, S. M., Ng, B. L., Prigmore, E., Fitzgerald, T. & Carter, N. P. Array 018-07885-5. painting: a protocol for the rapid analysis of aberrant chromosomes using DNA microarrays. Nat. Protoc. 4, 1722–1736 (2009). Competing interests: The authors declare no competing interests. 32. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows- Wheeler transform. Bioinformatics 25, 1754–1760 (2009). Reprints and permission information is available online at http://npg.nature.com/ 33. Kurtz, S. et al. Versatile and open software for comparing large genomes. reprintsandpermissions/ Genome. Biol. 5, R12 (2004). 34. Marco-Sola, S., Sammeth, M., Guigó, R. & Ribeca, P. The GEM mapper: fast, Journal peer review information: Nature Communications thanks the anonymous accurate and versatile alignment by filtration. Nat. Methods 9, 1185–1188 reviewers for their contribution to the peer review of this work. Peer reviewer reports are (2012). available. 35. Wu, T. D. & Watanabe, C. K. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics 21, 1859–1875 Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in (2005). published maps and institutional affiliations. 36. Slater, G. S. C. & Birney, E. Automated generation of heuristics for biological sequence comparison. BMC Bioinforma. 6, 31 (2005). 37. Merkel, A. et al. gemBS: high throughput processing for DNA methylation Open Access This article is licensed under a Creative Commons data from bisulfite sequencing. Bioinformatics https://doi.org/10.1093/ Attribution 4.0 International License, which permits use, sharing, bioinformatics/bty690 (2018). adaptation, distribution and reproduction in any medium or format, as long as you give 38. Li, H. Genome analysis Minimap2 : pairwise alignment for nucleotide appropriate credit to the original author(s) and the source, provide a link to the Creative sequences. Bioinformatics 34, 3094–3100 (2018). Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the Acknowledgements article’s Creative Commons license and your intended use is not permitted by statutory This study was supported by the Spanish Ministry of Economy and Competitiveness with regulation or exceeds the permitted use, you will need to obtain permission directly from Proyectos de I+D “Excelencia” y Proyectos de I+D+I “Retos Investigación” BFU2014- the copyright holder. To view a copy of this license, visit http://creativecommons.org/ 55090-P awarded to T.M.-B. and O.F., Centro de Excelencia Severo Ochoa 2013–2017 licenses/by/4.0/. and Centro de Excelencia Maria de Maeztu 2016–2019. We acknowledge the support from the CERCA Programme of the Generalitat de Catalunya, institutional support from the Spanish Ministry of Economy, Industry and Competitiveness (MEIC) through the © The Author(s) 2019 8 NATURE COMMUNICATIONS | (2019)10:4 | https://doi.org/10.1038/s41467-018-07885-5 | www.nature.com/naturecommunications