GoKawiilCentromere mapping and annotation
To identify completely and accurately assembled centromeres in whole-genome assemblies from diverse humans, we developed a pipeline called CenMAP (v.0.4.3.1; https://github.com/logsdon-lab/CenMAP). In brief, CenMAP aligns each whole-genome assembly to the centromeric regions of the CHM13 reference genome7 (v.2.0) using minimap2 (ref. 46) (v.2.28) with the following parameters: -x asm20 --secondary=no -s 25000 -K 15G. The resulting BAM files are then converted into BED files using rustybam (v.0.1.34; https://github.com/vollgerlab/rustybam) and the BED files are intersected with a BED file containing coordinates of unique regions within the pericentromeres that correspond to each chromosome using BEDTools47 (v.2.31.1). Contigs spanning the centromeres are subsequently identified and renamed to indicate the chromosome they most accurately align to based on the greatest number of alignments to the pericentromeric region. To confirm that the centromeric contigs contain α-satellite sequences, CenMAP runs dna-brnn48 (v.0.1-r65), which annotates α-satellite repeats within each centromeric contig. Neighbouring α-satellite repeat annotations are iteratively merged if they are within 10 kb of each other, yielding the coordinates of the centromeric α-satellite HOR array(s). The coordinates of the α-satellite HOR array(s) are extended on either side by approximately 500 kb, producing the coordinates of the centromeric regions within each contig.
Evaluation of contiguity and accuracy within centromeric regions
CenMAP evaluates the contiguity and accuracy of the centromeric regions by first extracting the regions using BEDTools and annotating them with RepeatMasker49 (v.4.1.0) using the following parameters: -engine rmblast -species human -dir {output_directory} -pa {threads} {input_contig}. CenMAP uses these annotations along with a custom tool we developed, CenStats (v.0.0.10; https://github.com/logsdon-lab/CenStats), to check the completeness of the contigs and reorient misoriented contigs from p to q arm relative to the CHM13 reference genome (v.2.0). This is achieved with the subcommand status, which compares the repeat sequence edit distance and Jaccard index relative to the reference genome, selecting the best matching orientation. Completely assembled centromeric contigs are subsequently reoriented to match the orientation of centromeres within the CHM13 reference genome (v.2.0) and then evaluated for the presence of assembly errors. To detect assembly errors in the centromeric regions, we developed NucFlag (v.0.3.3), an upgraded version of NucFreq50 that uses peak detection in the first- and second-highest read pileup nucleotide frequency to detect collapsed and misjoined sequences within the contigs. With NucFlag, PacBio HiFi reads are first aligned to the whole-genome assemblies using pbmm2 (v.1.13.1; https://github.com/PacificBiosciences/pbmm2) with the following parameters: --log-level DEBUG --preset SUBREAD --min-length 5000. Alignments are filtered using SAMtools51 (v.1.22) and the filter flag 2308, to remove secondary, supplementary and low-quality alignments. Only regions annotated as ‘ALR/Alpha’ by RepeatMasker are evaluated, and regions with nucleotide frequencies that are <95% of the mean of the highest nucleotide frequency are considered candidate misjoins, while regions with >225% of the mean highest nucleotide frequency are considered collapses. If the heterozygous site ratio exceeds 20% of the mean second nucleotide frequency, these regions are considered smaller collapses. Only centromeric regions that are free of misjoins, collapses, or other assembly errors are considered accurately assembled and retained for downstream analyses.
Annotation of centromeric regions
To determine the sequence composition and organization of all centromeric regions, CenMAP (v.0.4.3.1) uses a custom tool that we developed (Snakemake-HumAS-SD v.0.1.0; https://github.com/logsdon-lab/Snakemake-HumAS-SD) that leverages the α-satellite monomer HMM libraries from HumAS-HMMER (https://github.com/fedorrik/HumAS-HMMER_for_AnVIL) with the rapid annotation of StringDecomposer52 to efficiently annotate the centromeric regions. Snakemake-HumAS-SD first constructs a comprehensive α-satellite monomer library for each chromosome based on an existing HumAS-HMMER HMM profile (https://github.com/fedorrik/HumAS-HMMER_for_AnVIL/blob/main/AS-HORs-hmmer3.4-071024.hmm). Chromosomes with shared α-satellite monomers, such as chromosomes 1/5/19 as well as chromosomes 13/14/21/22, had their libraries merged. We developed an α-satellite HOR variation detection method that generates BED-format annotations of α-satellite HOR variants for each centromere based on the previous SV approach (https://github.com/fedorrik/stv). The method takes α-satellite monomer sequences as input and segments them based on three criteria: (1) adjacent monomers are separated by gaps >160 bp; (2) adjacent monomers are on opposite strands; or (3) monomer IDs are not monotonically increasing (on the forward strand) or decreasing (on the reverse strand). Compared with the original SV method, our approach can detect previously unresolvable α-satellite HOR structures, such as the newly identified 8-mer (S1C16H1L.2-9) expansion on chromosome 16, and annotates divergent α-satellite HORs that were not annotated using the previous SV method.
Estimation of the length of active α-satellite HOR arrays within centromeric regions
To estimate the length of active α-satellite HOR arrays, we ran CenStats (v.0.0.10) with the length subcommand, which leverages the annotations generated by Snakemake-HumAS-SD (described above) to output length statistics of each active α-satellite HOR array. CenStats first uses Snakemake-HumAS-SD α-satellite HOR SV BED files to identify the coordinates of active α-satellite HOR monomers (those with names containing ‘L’). Adjacent active α-satellite HOR monomers are then merged into one block if they are within around 1.5 α-satellite monomer lengths (256 bp) of each other, and blocks containing fewer than two α-satellite monomers are removed. To account for interruptions of stretches of α-satellite HORs caused by MEIs, adjacent blocks of active α-satellite HORs separated by less than 8 kb are iteratively merged into a single active α-satellite HOR array. Finally, the coordinates of α-satellite HOR arrays that consist of at least five α-satellite HOR monomers and have at least 90% of their total length composed of active α-satellite HORs is outputted as a BED file. To determine the orientation of active α-satellite HOR arrays, we followed a similar procedure, but incorporated both strand orientation and distance to group active α-satellite HOR monomers into arrays.
Detection of CDRs within centromeric regions
In addition to annotating centromeric sequences, CenMAP (v.0.4.3.1) uses methylation information in native long-read sequencing data to annotate the CpG methylation profile across each centromeric region. CenMAP first aligns native ONT long-read sequencing data from each genome to the relevant whole-genome assembly using minimap2 (ref. 46) (v.2.28) with the following command: minimap2 -x lr:hqae -s 4000 -y -I 8 G {ref} {ont_reads}. CenMAP filters the alignments using SAMtools51 (v.1.9) -F 2308 to retain only primary alignments and remove partial and secondary alignments. Finally, CenMAP uses CDR-Finder31 (v.1.0.4), a tool to detect hypomethylated regions termed CDRs22, spanning from the first hypomethylated base to the last hypomethylated base and inclusive of the sequence in between. We manually curated the CDR calls by checking each region outputted in the BED file using IGV53, removing false positives, recovering missed calls and merging calls that were within 150 kb of each other.
... continue reading