Support for references larger than 4.29 Gbp

[xscapintime]

Hi, I am trying to use ABISMAL with a 6.06-Gbp phased diploid T2T reference. The index is created, but loading it fails with:

[Tue Aug 11 19:49:33 2026] loading index /home/user/genome/idx/abismal.h9/H9_T2T.abismalidx
terminate called after throwing an instance of 'std::length_error'
  what():  cannot create std::vector larger than max_size()

I noticed that genome-wide coordinates appear to use std::uint32_t, whose maximum value is about 4.29 billion. Could this be caused by coordinate overflow? Would you consider adding 64-bit coordinate support for large diploid references?

Thank you!

[andrewdavidsmith]

@xscapintime Your guess is likely correct. I'll look into it. You can try it yourself, but likely that specific uint32_t is interdependent on another one somewhere else. At one point we were keeping a lot of 32-bit values, and changing those to 64-bit would obviously double the space. But I think we moved past that implementation detail, so it might not matter anymore. Also, there might be an easy solution for your situation: you basically have the whole diploid, which we might be able to do with two indices. And I think there should be a lot of redundancy in the diploid so we might be saving space in the index due to that.

Can you post the link to the downloadable reference so I'm working with the exact same one? I have a T2T somewhere, but I think it was the first published one, and the same size as hg38.

[xscapintime]

Thank you for your reply! The diploid T2T assembly I am using is T2T-H9. I downloaded the two haplotypes from NCBI:

wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/054/883/195/GCA_054883195.1_H9_T2T.hap1/GCA_054883195.1_H9_T2T.hap1_genomic.fna.gz
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/054/883/265/GCA_054883265.1_H9_T2T.hap2/GCA_054883265.1_H9_T2T.hap2_genomic.fna.gz

I renamed the chromosomes as chr1_hap1, chr1_hap2, etc., and concatenated the two FASTA files, with all hap1 chromosomes followed by all hap2 chromosomes:

cat H9_T2T.hap1.renamed.fa H9_T2T.hap2.renamed.fa \
    > H9_T2T.renamed.fa

The combined reference contains 46 nuclear chromosomes and has a total length of 6,058,641,370 bp.

By “two indices,” do you mean mapping separately to hap1 and hap2 and comparing the results afterward, or having Abismal search two haplotype indices jointly during a single mapping run? Joint searching would be ideal because I would like the two haplotypes to compete for the best alignment.

In the meantime, I tried mapping directly from the combined FASTA using -g, without loading the prebuilt index:

time abismal -v      -t 8  -c 50000 -g $ref -s ${out}_abismal_c50k_se_stats.txt  -o ${out}_abismal_c50k_se.bam $p1

This completed without an error, but the alignments appeared to be truncated:

total_reads: 155956136
mapped:
    num_mapped: 85777503
    num_unique: 73425257
    num_ambiguous: 12352246
    percent_mapped: 55.001
    percent_unique: 47.0807
    percent_ambiguous: 7.92033
    unique_error:
        edits: 70514658
        total_bases: 7211163938
        error_rate: 0.00977854
num_unmapped: 70178633
num_skipped: 1110325
percent_unmapped: 44.999
percent_skipped: 0.711947

samtools idxstats ${out}_abismal_c50k_se_format.sorted.bam
chr1_hap1	250382078	9906830	0
chr2_hap1	242438398	10532509	0
chr3_hap1	200940782	8630654	0
chr4_hap1	193253728	8228797	0
chr5_hap1	186212366	7778117	0
chr6_hap1	173051588	7351106	0
chr7_hap1	161716607	7142411	0
chr8_hap1	145465483	6003222	0
chr9_hap1	135918162	5208648	0
chr10_hap1	134385439	3152072	0
chr11_hap1	133922055	0	0
chr12_hap1	133376810	0	0
...
chr21_hap2	42759301	0	0
chr22_hap2	48623973	0	0
chrX_hap2	153704768	0	0
*	0	0	0

Mapping to the haploid T2T-CHM13 reference worked fine previously, so I checked the genome sizes:

cat chm13v2.0.fa.fai |awk '{sum+=$2} END {print sum}'
3117292070

cat H9_T2T.renamed.fa.fai | head -n 11 | awk '{sum+=$2} END {print sum}'
1957686686

The alignments seemed to stop at chr11_hap1, although the cumulative genome size at this point had not yet exceeded the uint32_t limit.

I wonder what you make of this and would appreciate any suggestions you may have. Thank you very much!

[andrewdavidsmith]

@xscapintime Thanks, this will help! My comment was really just about how I might be able to make this work internally without too many changes. I would still want to have the mapping identify the best among both haplotypes, as you mention. I might have been thinking out loud a bit.

[andrewdavidsmith]

@xscapintime There's now a branch named typedefs-for-genome-indexing-values which has the following in AbismalIndex.hpp:

using genome_pos_t = std::uint32_t;  // Position in a genome

If you checkout the branch, change this value to be std::uint64_t, then build the code, it should (might?) do what you need. I was only able to confirm that the 64-bit values pass the tests setup for the 32-bit values, but obviously the binary index file differs. I think it's about 2x the size as expected.

I also simulated reads from the genome the way you explained it, and then mapped them back. They had no sequencing errors, no methylation, and very high bisulfite conversion. The fraction of reads mapping uniquely in my test was about 55% I think. So I assume that's what you might expect with real data. Again, my tests were very limited. Even if the code doesn't help you, it gave me the opportunity to make sure variables used for different purposes are now marked according to their use.

[xscapintime]

Thank you so much!

I changed:

using genome_pos_t = std::uint32_t;  // Position in a genome

to:

using genome_pos_t = std::uint64_t;

The new index is larger:

6.7G Aug 14 14:07 H9_T2T.64.abismalidx
3.7G Aug 11 19:42 H9_T2T.abismalidx

Mapping to the new index completed without an error, but the resulting alignments are still truncated:

samtools idxstats ${out}_abismal_c50k_se_format.sorted.bam
chr1_hap1	250382078	246562	0
chr2_hap1	242438398	256154	0

...
chr22_hap1	47079544	1208581	0
chrX_hap1	155153153	4325848	0
chr1_hap2	243198671	246903	0
chr2_hap2	242142825	256985	0
chr3_hap2	202683022	223478	0
chr4_hap2	191171759	234347	0
chr5_hap2	182409348	206285	0
chr6_hap2	172871300	210293	0
chr7_hap2	161232845	47891	0
chr8_hap2	145868755	0	0
chr9_hap2	135032811	0	0
chr10_hap2	135648429	0	0
...
chr21_hap2	42759301	0	0
chr22_hap2	48623973	0	0
chrX_hap2	153704768	0	0
*	0	0	0

The alignments appear to stop partway through chr7_hap2, near the cumulative position of 4.3 Gbp:

cat H9_T2T.renamed.fa.fai | head -n 29 | awk '{sum+=$2} END {print sum}'
4274795030
cat H9_T2T.renamed.fa.fai | head -n 30 | awk '{sum+=$2} END {print sum}'
4436027875

This cutoff is very close to the uint32_t limit of 4,294,967,295. I was wondering whether there might be other variables in the code that also need to be changed to std::uint64_t.
I would appreciate any thoughts or suggestions. Thank you again!

[andrewdavidsmith]

Yes, I see the problem, and other symptoms pointing in the same direction. It might take me some time to get back to this. In the meantime, why not map to both haplotypes, and use samtools sort based on read name, then take the better of the two mappers by iterating once through simultaneously in both files? I think that would have similar accuracy, and probably wouldn't take very long unless you have a huge amount of data.