Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

mpileup error messages #1652

Closed
hd2326 opened this issue May 4, 2022 · 17 comments · Fixed by samtools/htslib#1430
Closed

mpileup error messages #1652

hd2326 opened this issue May 4, 2022 · 17 comments · Fixed by samtools/htslib#1430

Comments

@hd2326
Copy link

hd2326 commented May 4, 2022

Are you using the latest version of samtools and HTSlib? If not, please specify.

(run samtools --version)

(built as a singularity container)

samtools 1.15
Using htslib 1.15
Copyright (C) 2022 Genome Research Ltd.

Please describe your environment.

  • machine architecture (run uname -m on Linux/Mac OS or wmic os get OSArchitecture on Windows)

x86_64 (architecture of the singularity container)

  • compiler (run gcc --version or clang --version)

gcc (Debian 10.2.1-6) 10.2.1 20210110
Copyright (C) 2020 Free Software Foundation, Inc.

(gcc of the singularity container)

Please specify the steps taken to generate the issue, the command you are running and the relevant output.

Greetings!

As I am running the mpileup module across multiple bam files, some bam files might produce one of the following errors:

  1. malloc(): invalid size (unsorted)
  2. malloc(): invalid next size (unsorted)
  3. malloc(): corrupted top size
  4. free(): invalid next size (fast)
  5. corrupted size vs. prev_size while consolidating

The code I ran:

for file in *.bam; do
samtools view -h -b -F 4 $file > $file.mapped.bam
samtools mpileup --output-mods --positions dna.bed $file.mapped.bam > $file.pileup
done

It seems that mpileup breaks when there is an error, and produces a truncated pileup file. Any insights on what's going on? Thank you very much in advance for your help!

@jkbonfield
Copy link
Contributor

I'm guessing this is a bug in the --output-mods part as the rest has been more extensively tested.

Is it possible to do this on successively smaller bam files (easiest done in SAM) to find a small test case that triggers it? I'll explore some more at my end, but getting public data with mods in it is challenging, so most of my testing was with made up values rather than real data.

If you're confident, it's also possible to build with more debugging information which would spot memory corruptions. Not for production usage, but this can pinpoint errors more precisely (given a sufficiently modern compiler; can use clang too):

make clean
make CC="gcc -fsanitize=address" CFLAGS=-g

jkbonfield added a commit to jkbonfield/htslib that referenced this issue May 4, 2022
This attempted to grow memory by the maximum amount of space a base
modification would take up, but due to a misunderstanding of kstring
it kept adding this to the original size rather than actually growing
the allocated size.

(Probably) fixes samtools/samtools#1652
@jkbonfield
Copy link
Contributor

I've found a bug in htslib which likely triggers this problem. See samtools/htslib#1430

It may undergo other improvements yet, but it would be helpful if you could confirm whether this change cures your issue. It's a small two line change to htslib/sam.c:

diff --git a/sam.c b/sam.c
index dd1b7d9..865b55f 100644
--- a/sam.c
+++ b/sam.c
@@ -5306,6 +5306,7 @@ int bam_plp_insertion_mod(const bam_pileup1_t *p,
                 hts_base_mod mod[256];
                 if (m && (nm = bam_mods_at_qpos(p->b, p->qpos + j - p->is_del,
                                                 m, mod, 256)) > 0) {
+                    int o_indel = indel;
                     if (ks_resize(ins, ins->l + nm*16+3) < 0)
                         return -1;
                     ins->s[indel++] = '[';
@@ -5329,6 +5330,7 @@ int bam_plp_insertion_mod(const bam_pileup1_t *p,
                                              qual);
                     }
                     ins->s[indel++] = ']';
+                    ins->l += indel - o_indel; // grow by amount we used
                 }
             }
             break;

@hd2326
Copy link
Author

hd2326 commented May 4, 2022

@jkbonfield

Thank you so much for the quick response!

Just want to make sure that I need to change the bam_plp_insertion_mod function in ./htslib-1.15/sam.c.

Following is my current bam_plp_insertion_mod function, just want to make sure we are on the same page. If possible, would you please give me the entire updated bam_plp_insertion_mod function? Thank you very much!

int bam_plp_insertion_mod(const bam_pileup1_t *p,
hts_base_mod_state *m,
kstring_t *ins, int *del_len) {
int j, k, indel, nb = 0;
uint32_t *cigar;

if (p->indel <= 0) {
    if (ks_resize(ins, 1) < 0)
        return -1;
    ins->l = 0;
    ins->s[0] = '\0';
    return 0;
}

if (del_len)
    *del_len = 0;

// Measure indel length including pads
indel = 0;
k = p->cigar_ind+1;
cigar = bam_get_cigar(p->b);
while (k < p->b->core.n_cigar) {
    switch (cigar[k] & BAM_CIGAR_MASK) {
    case BAM_CPAD:
    case BAM_CINS:
        indel += (cigar[k] >> BAM_CIGAR_SHIFT);
        break;
    default:
        k = p->b->core.n_cigar;
        break;
    }
    k++;
}
nb = ins->l = indel;

// Produce sequence
if (ks_resize(ins, indel+1) < 0)
    return -1;
indel = 0;
k = p->cigar_ind+1;
j = 1;
while (k < p->b->core.n_cigar) {
    int l, c;
    switch (cigar[k] & BAM_CIGAR_MASK) {
    case BAM_CPAD:
        for (l = 0; l < (cigar[k]>>BAM_CIGAR_SHIFT); l++)
            ins->s[indel++] = '*';
        break;
    case BAM_CINS:
        for (l = 0; l < (cigar[k]>>BAM_CIGAR_SHIFT); l++, j++) {
            c = seq_nt16_str[bam_seqi(bam_get_seq(p->b),
                                      p->qpos + j - p->is_del)];
            ins->s[indel++] = c;
            int nm;
            hts_base_mod mod[256];
            if (m && (nm = bam_mods_at_qpos(p->b, p->qpos + j - p->is_del,
                                            m, mod, 256)) > 0) {
                if (ks_resize(ins, ins->l + nm*16+3) < 0)
                    return -1;
                ins->s[indel++] = '[';
                int j;
                for (j = 0; j < nm; j++) {
                    char qual[20];
                    if (mod[j].qual >= 0)
                        sprintf(qual, "%d", mod[j].qual);
                    else
                        *qual=0;
                    if (mod[j].modified_base < 0)
                        // ChEBI
                        indel += sprintf(&ins->s[indel], "%c(%d)%s",
                                         "+-"[mod[j].strand],
                                         -mod[j].modified_base,
                                         qual);
                    else
                        indel += sprintf(&ins->s[indel], "%c%c%s",
                                         "+-"[mod[j].strand],
                                         mod[j].modified_base,
                                         qual);
                }
                ins->s[indel++] = ']';
            }
        }
        break;
    case BAM_CDEL:
        // eg cigar 1M2I1D gives mpileup output in T+2AA-1C style
        if (del_len)
            *del_len = cigar[k]>>BAM_CIGAR_SHIFT;
        // fall through
    default:
        k = p->b->core.n_cigar;
        break;
    }
    k++;
}
ins->s[indel] = '\0';
ins->l = indel; // string length

return nb;      // base length

}

@jkbonfield
Copy link
Contributor

jkbonfield commented May 5, 2022

The easiest way of editing it is simply to cut and paste the quoted diff above and save to e.g. file p and use patch on it:

$ patch < p
patching file sam.c

You can also see the whole function in sam.c file from the git PR with the fix in it:

https://github.com/jkbonfield/htslib/blob/a253c0af2de9e8659ed7d5d9ab1d215df0c9b749/sam.c#L5242-L5352

jkbonfield added a commit to jkbonfield/htslib that referenced this issue May 5, 2022
This attempted to grow memory by the maximum amount of space a base
modification would take up, but due to a misunderstanding of kstring
it kept adding this to the original size rather than actually growing
the allocated size.

(Probably) fixes samtools/samtools#1652
@hd2326
Copy link
Author

hd2326 commented May 5, 2022

Thank you @jkbonfield for the instructions!

I fixed the function as you suggested, but still get the errors (including some new errors):

  1. corrupted size vs. prev_size
  2. double free or corruption (!prev)
  3. free(): corrupted unsorted chunks
  4. free(): invalid next size (fast)
  5. malloc(): corrupted top size
  6. malloc(): invalid next size (unsorted)
  7. malloc(): invalid size (unsorted)
  8. malloc(): mismatching next->prev_size (unsorted)
  9. mremap_chunk(): invalid pointer
  10. munmap_chunk(): invalid pointer
  11. samtools: malloc.c:2396: sysmalloc: Assertion `(old_top == initial_top (av) && old_size == 0) || ((unsigned long) (old_size) >= MINSIZE && prev_inuse (old_top) && ((unsigned long) old_end & (pagesize - 1)) == 0)' failed.
  12. Segmentation fault (core dumped)

Just FYI, I'm processing some long nanopore sequencing reads. Maybe that causes the problem?

@jkbonfield
Copy link
Contributor

Is any of your data public? It's hard to diagnose and fix the issue without having data to trigger?

It may only need 1 read also. How soon does it die? Is it possible to use head/tail on a SAM version of your input file to gradually reduce it down to get a minimal data set that triggers the issue?

@hd2326
Copy link
Author

hd2326 commented May 10, 2022

Hey @jkbonfield the data I analyzed is private. So maybe I upload the data to Google drive and share with you? May I have your gmail account?

Since almost all the reads gave error messages and the messages are different, maybe we should look into the merged bam? Any thoughts?

P.S. It seems that mpileup gives error messages when analyzing reads with abundant modifications (Mm tag).

@jkbonfield
Copy link
Contributor

The messages being different is probably irrelevant, as any malloc related error message is most likely the same bug showing up in different ways. If you can reproduce it with a single read then that makes for an easy to debug test case.

You can email me direct using my gmail.com email, which has the same username as my github handle.

Thanks

@hd2326
Copy link
Author

hd2326 commented May 11, 2022

Just emailed you the example read:)

@hd2326
Copy link
Author

hd2326 commented May 11, 2022 via email

@jkbonfield
Copy link
Contributor

I've downloaded your test data (thank you) and can reproduce the issue before I apply samtools/htslib#1430, but not after.

I also tried running it under Valgrind and built using "gcc -fsanitize=address". Neither tool identifies an issue. So I'm confident that the bug I fixed was the same one you are seeing.

Did you recompile both htslib and samtools after applying the fix? Are you certain you are running the correct samtools binary?

@hd2326
Copy link
Author

hd2326 commented May 12, 2022

Ahhh I see! I didn't re-compile samtools and that might be the problem. Will try and keep you posted!

@hd2326
Copy link
Author

hd2326 commented May 13, 2022

Yes it works! I did extensive testing and got no error! Thank you so much for the help @jkbonfield !

@hd2326 hd2326 closed this as completed May 13, 2022
@jkbonfield
Copy link
Contributor

Please to have resolved it. Thanks for the aid in debugging this.

@hd2326
Copy link
Author

hd2326 commented May 13, 2022

Sure!

BTW, will the fix be included in the next release? When would be the next release? Thank you very much!

@jkbonfield
Copy link
Contributor

It's not merged in yet, but as it's an obvious bug with a simple fix I don't see any blockers to it going through review and merge so I'm sure it'll be in the next release.

Our typical release cycle is every 3-4 months, but it can vary a little bit sometimes. Based on that "summer" sounds sufficiently broad to be an accurate guess on the release. I can't say more accurately.

@hd2326
Copy link
Author

hd2326 commented May 13, 2022

Got it! Thank you so much @jkbonfield !

whitwham pushed a commit to samtools/htslib that referenced this issue May 18, 2022
This attempted to grow memory by the maximum amount of space a base
modification would take up, but due to a misunderstanding of kstring
it kept adding this to the original size rather than actually growing
the allocated size.

(Probably) fixes samtools/samtools#1652
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants