BAM indexes for very large contigs

[mvdbeek]

from https://biostar.usegalaxy.org/p/29884/

Dear community,

is there a way to modify the default way Galaxy creates the index *.BAI files? We work with a large salamander genome, therefore, quite some coordinates in the BAM files are beyond 2^30. Therefore, the index files have to be created using either biobambam2:bamSort or samtools index -c -m 14. However, I cannot really figure out if there is a simple way to "teach" Galaxy how to do that. We have our own institute's Galaxy instance, therefore, I can change pretty much anything there...

I am sorry if a similar question was posted before and would be thankful for any advice.

Best regards, Sergej

[mvdbeek]

I guess there's a couple of options.

• The easiest and lamest is if we change set_meta to not fail if an index can't be generated. That means users can generate BAM files by any tool without further adjustment, but tools or display applications (genome browser etc.) that require a bam index will fail.

• We could just replace the .bai index with a .csi index. I don't think that is a good idea as they are slightly different, though htslib seems to not care about this. Not everyone uses htslib, so that's risky, certainly if we can't exclude that we silently skip some data.

• We could generate .bai and .csi indexes (The entire file needs to be read, so this isn't without cost), or if we generate .csi files if the contig size exceeds the maximum for .bai indexes (we can check this in the bam header, don't need to read the entire file).

I think for now it's probably okay to add an additional optional csi_index metadata element and store the csi_index when a bam index can't be generated. That means tool authors will need to check if a csi_index or a bam_index exists. I haven't looked into how display_applications interact with the bam_index, but I guess we can use whatever index exists.

[SergejN]

I would suggest generating bai if possible and csi, otherwise, and maybe renaming csi to bai (I'm not quite sure about it, but many tools seem to accept that). The rationale behind it is that htslib can handle it. And tools that don't use htslib to read bam/bai files, would fail for large contigs anyway.
But I do agree that a more generic solution would be better. in this regard, I like the last paragraph of your post - a metadata element that can be checked for in the tool is quite handy.

[mvdbeek]

Started looking into this, but it appears there are some more issues with coordinate beyond 2**30, some of which may be pysam or htslib issues (pysam-developers/pysam#741).
I've got some WIP in mvdbeek@e386b8d, may need to fall back to parsing the text representation of the header.

[ieguinoa]

Hi @mvdbeek ,

I would like to comment something on this issue. I'm also experiencing it with the wheat genome (https://wheat-urgi.versailles.inra.fr/Seq-Repository/Assemblies).
In my case though I found that I could build the standard bai index of the corresponding bam files using samtools (outside of Galaxy) but pysam fails so I get an error in the output bam when running set_metadata. I'm not sure if the latest changes in pysam (pysam-developers/pysam#732) would solve this.
Just to note that there are some cases in between where maybe the bai index can still be built and there is no need for a more complicated workaround building the .csi one.

Ignacio

[mvdbeek]

What is the longest contig in your assembly @ieguinoa ?
What error do you see in pysam ?

[ieguinoa]

What is the longest contig in your assembly @ieguinoa ?
What error do you see in pysam ?

The longest contig is 844676927 bp.

The error is:

[E::hts_idx_push] Region 536882295..536882446 cannot be stored in a bai index. Try using a csi index with min_shift = 14, n_lvls >= 6
Traceback (most recent call last):
File "index_bam_pysam.py", line 4, in
pysam.index(sys.argv[1], sys.argv[2])
File "miniconda2/envs/pysam/lib/python2.7/site-packages/pysam/utils.py", line 75, in call
stderr))
pysam.utils.SamtoolsError: 'samtools returned with error 1: stdout=, stderr=samtools index: failed to create index for "NG-18597_lib311030_6332_1_hisat2_mapping.bam": Numerical result out of range\n'

[AjitPS]

@rob123king and I also have seen issues trying to use mappers like hisat2 (that generate bam files using a BAI index) for the full wheat genome. Would be useful to add or see how we can add an option to the tool/mapper (hisat2, etc.) to use CSI index instead.

i.e., to change pysam for mappers like hisat2 to generate .csi (instead of the default .bai index) when the .bam file is too big.

As currently, hisat2 + featureCounts will work but downstream SNP calling will not be correct due to the hisat2 output bam losing the end exceeding the bai index limit.
@mvdbeek

[nsoranzo]

Unfortunately this needs to wait for a new release of pysam, the latest is 0.15.2 and does not contain pysam-developers/pysam#755 .

[FredericBGA]

The wheat genome is also is this case, and needs .csi indexes.

0.15.3 is here now, and the PR of @mvdbeek has been merged!

[FredericBGA]

I can do that and submit a PR with a little help!

[mvdbeek]

It may make sense to start from mvdbeek@e386b8d, but haven't looked into this since then. Let me know if you need any help.

[mvdbeek]

There are also no wheels yet for 0.15.3, so that'd be the first step, pysam-developers/pysam#831 is the problem here

[FredericBGA]

thank you, but I 'm afraid that I won't be able to overcome pysam-developers/pysam#831 :(