SAM file is not sorted ERROR even if the input file is sorted

[YuanfengZhang]

Describe the bug
SAM file is not sorted ERROR even if the input file is sorted

cd result/BL/trim-galore/no-pre-dedup/no-pre-calibrate/segemehl/samtools/no-calibrate
mkdir -p dnmtools
samtools sort -@ 16 -o BL.sorted.bam BL.bam
dnmtools counts \
    -c /data/BL.fa \
    -t 16  \
    -o dnmtools/BL.dnmtools.tmp \
    BL.sorted.bam
cd dnmtools

awk -F'\t' '$6 > 4' \
    BL.dnmtools.tmp \
    | pigz  -p 16 --best > BL.dnmtools.gz

stdout:

[bam_sort_core] merging from 0 files and 16 in-memory blocks...
[error: 0][ERRNO: 0][Success][SAM file is not sorted]

I'm sure that this error is from dnmtools instead of samtools

dnmtools/src/analysis/methcounts.cpp

Line 464 in e18ae72

if (tid < prev_tid) throw dnmt_error("SAM file is not sorted");

I tried samtools index -n but it failed too. I wonder what exact sorting tool and method does dnmtools require?

[andrewdavidsmith]

@YuanfengZhang This should be easy to fix. First, most times when we have a sorting error it's due to locale when sorting. I think the docs mention something about this, but I suggest doing a quick search of locale and sorting. Iirc we suggest to use LC_ALL=C, but if the sorting is done inside samtools then I'm not sure that would be the issue.

Since you've identified the place where the error is caught, are you able to make a change to the code and print the tid and prev_tid? That might give more clues. I'm happy to debug if you can send me any test data that allows me to reproduce the issue.

Edit: I just looked at the code and that error only arises at bounds between chromosomes. So if you run it with -progress it might show which chromosome transition is causing the problem, and give us some hints.

[YuanfengZhang]

processing chr1
processing chr10
processing chr11
processing chr12
processing chr13
processing chr14
processing chr15
processing chr16
processing chr17
processing chr18
processing chr19
processing chr2
processing chr20
processing chr21
processing chr22
processing chr3
processing chr4
processing chr5
processing chr6
processing chr7
processing chr8
processing chr9
processing chrM
processing chrX
processing chrY
processing lambda
processing pUC19
[error: 0][ERRNO: 0][Success][SAM file is not sorted]

[YuanfengZhang]

❯ samtools view -H BL.sorted.bam | rg 'pUC19'         
@SQ     SN:pUC19        LN:2686

This is no records on pUC19 contig in the .bam file, which means it should be skipped instead of raising an error. Am I right? OR it's because dnmtools cannot handle spike-in contigs?

[andrewdavidsmith]

It should work as-is. The target sequence is in the SAM header, so it will be treated like a chromosome. There should be no problem if a target sequence has no mapped reads. Is there any way you can give me a failing test case? I don't mind big downloads. You can email me directly. Also, I should mention that the progress is printed after that error condition is tested, so the issue would happen at the transition from pUC19 to the reference sequence (TID) that comes after it. I'm signing off but I can look at this further tomorrow.

[YuanfengZhang]

Thank you very much. I find out that this error is not about the content of the .bam file nor aligners, because it can be reproduced with minimal ones from bwa-meme / bwa-meth / segemehl output .bam files. I have sent (from zhangyuanfeng1997@foxmail.com) the paierd .fq.gz, .sam (for convenience) and spike-in.fa to your email. Please check them.

As for the mapping, deduplication and sorting command history, please check the header of the .bam file.
To prepare the BL.fa, please run cat GRCh38_no_alt.fa spike-in.fa > BL.fa.