bug counting alignments in cram files

[matthdsm]

Hi,

I noticed a discrepancy when loading a cram vs bam file.
BAM:

• pysam v0.16.0.1

reads_file = pysam.AlignmentFile(infile, "rb")
reads_file.mapped
176197222
reads_file.unmapped
2958898

• samtools v1.12

# mapped reads
samtools view -c -F 260 sample.bam 
176197222

# unmapped reads
samtools view -c -f 4 sample.bam 
2958898

CRAM:

• pysam v0.16.0.1

reads_file = pysam.AlignmentFile(infile, "rc", reference_filename=reference)
reads_file.mapped
0
reads_file.unmapped
6073771974947137635

• samtools v1.12

# mapped reads
samtools view -c -F 260  -T reference.fa sample.cram 
176197222

# unmapped reads
samtools view -c -f 4   -T reference.fa sample.cram 
2958898

Is there a way to fix this? HTSlib seems to parse everything correctly

Cheers
M

[jmarshall]

CRAM indexes do not contain the mapped/unmapped counts that are being used to calculate these total count statistics. So the reasonable thing to do would be to leave the totals as 0 / None / raise an exception in this case.

[matthdsm]

Hi John,

Ok, makes sense, but what about the huge number of unmapped reads?
It's way bigger than the total number of reads.

Thanks
M

[jmarshall]

It's almost as though it's returning uninitialised memory…
This is a minor HTSlib bug; see samtools/htslib#1335.