Could not retrieve index file for primary.bam

[eboileau]

Description

I wish to use the latest NanoSim v3.0.1 (supports reading .gz sequence files, and bam files), but read_analysis.py does not complete. The primary.bam file is not indexed, an that might be an issue for pysam, but there seem to be more, this might be related to #129.

Error

2021-07-14 12:59:14: Processing alignment file: bam
[W::hts_idx_load3] The index file is older than the data file: analysis/Nanopore.bam.bai
021-07-14 12:59:16: Aligned reads analysis
[E::idx_find_and_load] Could not retrieve index file for 'analysis/nanosim_model/sim_primary.bam'

and further down

2021-07-14 12:59:17: match and error models
[E::idx_find_and_load] Could not retrieve index file for 'analysis/nanosim_model/sim_primary.bam'
Traceback (most recent call last):
  File "/beegfs/homes/eboileau/.miniconda3/envs/scNapBar-dev/bin/besthit_to_histogram.py", line 318, in hist
    cs_string = alnm.get_tag('cs')
  File "pysam/libcalignedsegment.pyx", line 2399, in pysam.libcalignedsegment.AlignedSegment.get_tag
  File "pysam/libcalignedsegment.pyx", line 2438, in pysam.libcalignedsegment.AlignedSegment.get_tag
KeyError: "tag 'cs' not present"

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/eboileau/.miniconda3/envs/scNapBar-dev/bin/read_analysis.py", line 720, in <module>
    main()
  File "/home/eboileau/.miniconda3/envs/scNapBar-dev/bin/read_analysis.py", line 710, in main
    error_model.hist(prefix, alnm_ext)
  File "/beegfs/homes/eboileau/.miniconda3/envs/scNapBar-dev/bin/besthit_to_histogram.py", line 320, in hist
    cs_string = get_cs(alnm.original_sam_line.split()[5], alnm.get_tag('MD'))
AttributeError: 'pysam.libcalignedsegment.AlignedSegment' object has no attribute 'original_sam_line'

Expected behavior

read_analysis.py completes successfully, and generates the model to be used as input for simulator.py.

To reproduce

read_analysis.py genome -i Nanopore.fq.gz -ga Nanopore.bam -o nanosim_model/sim

and this also occurs with all uncompressed input files (I guess this is expected, since NanoSim now outputs compressed files anyway)

read_analysis.py genome -i Nanopore.fq -ga Nanopore.sam -o nanosim_model/sim

However, using NanoSim 2.5.0, the latest command is successful.

Environment

Python 3.7.6
conda 4.9.2
NanoSim 3.0.1 ( but Note that the version has not been updated in some scripts, e.g. read_analysis.py --version return NanoSim 3.0.0, although I am using 3.0.1 )
pysam 0.16.0.1 (samtools 1.10, using htslib 1.10.2)

[cheny19]

Hi @eboileau , this is a different error which I cannot reproduce with your command. My gut feeling is this error happens with this line The index file is older than the data file: analysis/Nanopore.bam.bai. Could you try running read_analysis.py without your own alignment file? If you have to, make sure it contains either 1) CS tag or 2) both MD tag and cigar string.

[kmnip]

The primary.bam file is not indexed, an that might be an issue for pysam

Pysam can read BAM files line by line without an index file just fine. The index file is only required for random access on the BAM file. So, the message from Pysam, [E::idx_find_and_load] Could not retrieve index file, can be safely ignored. You would see this message in a successful run of NanoSim.

[eboileau]

Thanks for your prompt reply. I can reproduce this error with different input files. The error actually arises from the primary.bam file that is generated by NanoSim.

I can reproduce this error with v3.0.1.
I also downloaded v3.0.0, as well as v3.0.0-beta, and I get the same error (but using uncompressed files as input).
With v2.6.0 and v2.5.0, I have no error, and read_analysis.py finishes successfully, so this seems to be due to changes in v3.0.0-beta.

[eboileau]

@kmnip right, without an index, random access via fetch() and pileup() is disabled.

[eboileau]

I attach a small test set of about 200 Nanopore reads, mapped with

minimap2 -v1 -t 6 -ax splice --MD -ub /biodb/genomes/homo_sapiens/GRCh38_96/GRCh38_96.fa sample.fastq.gz | samtools sort - -@6 -T /scratch/global_tmp/ -o sample.bam

The NanoSim command to reproduce the error (v3.0.0-beta, v3.0.0, and v3.0.1) is

read_analysis.py genome -i sample.fastq.gz -ga sample.bam -o sim

test.tar.gz

[eboileau]

Actually, I was using NanoSim version release v3.0.0-beta, v3.0.0, and v3.0.1 (tags). But if I clone master, and run read_analysis.py, I have no error, only the warning [E::idx_find_and_load] Could not retrieve index file for 'sim_primary.bam.
@cheny19 this might explain why you could not reproduce the error, if you were using the master branch.

It would be nice to be able to use the latest version release without error, and eventually be able to install it via conda. Thanks.

[cheny19]

I see. It seems there's no much difference with the code in master branch and v3.0.1, so I guess the error has something to do with conda installation.

[eboileau]

v3.0.1 is not in Anaconda. I tested all releases (v3.0.0-beta, v3.0.0, v3.0.1, and v2.0.6) by just downloading the tag, and running the python script read_analysis.py. I did the same for master, I cloned, and ran the script. The rest of the environment remained the same for all tests. Only v2.0.5 was installed via conda, but v.2.05, v2.0.6, and master all were fine. I doubt this has to do with conda... And since dependencies remained the same, this must be something in NanoSim...

In any case, if master has no error, do you think it would be possible to create a new release soon? I don't know what are your development plans, but that would be great. Since I use NanoSim in my workflow, I would rather pin a specific release for reproducibility.

[eboileau]

By the way, were you able to reproduce and trace the error with the small dataset I uploaded, using releases v3.0.0-beta, v3.0.0, v3.0.1 ?

[cheny19]

Oh, I finally figured out what went wrong! It's the pysam error when retrieving cs tag and MD tag. It was mentioned in another issue too, and fixed by the most recent commit, which is why the master branch works now. I can quickly make a release so you can tag the current master branch for your pipeline.

[eboileau]

Awesome, thanks! I'm waiting to see the next release.

[eboileau]

Thanks for updating bioconda with v3.0.2.