Skip to content

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

Open

Description

@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

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?

Activity

  1. andrewdavidsmith commented on May 12, 2025

    @andrewdavidsmith
    Collaborator

    @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.

  2. YuanfengZhang commented on May 12, 2025

    @YuanfengZhang
    Author

    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]

  3. YuanfengZhang commented on May 12, 2025

    @YuanfengZhang
    Author
    ❯ 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?

  4. andrewdavidsmith commented on May 12, 2025

    @andrewdavidsmith
    Collaborator

    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.

  5. YuanfengZhang commented on May 12, 2025

    @YuanfengZhang
    Author

    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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions