Skip to content

Error in modkit pileup --phased #552

@jkh00

Description

@jkh00

Hi,

I'm running modkit pileup v0.6.0 on a sorted phased bam:

modkit pileup --phased --modified-bases 5mC --prefix test --ref genome.fasta sorted.phased.bam out

and I ran into the following error:

discarded 0 contigs with zero aligned reads
parsed 1 base modification(s). Base modifictions other than 'C:m' will be counted as 'N_other'.
producing phased output
attempting to sample 10042 reads
Using filter threshold 0.8652344 for C.
using optimized workers for A,C all-context
[00:00:00] ---------------------------------------- 0/40001 genome positions 0/s 0s
0 rows written
0 ~records errored
thread '' panicked at modkit-core/src/pileup/mod.rs:39:9:
assertion left == right failed
left: 0
Done, processed 0 rows.
thread 'main' panicked at modkit-core/src/pileup/subcommand.rs:1379:34:
worker thread 0 paniced: Any { .. }

when I ran it without --phased, it completed without error:

discarded 0 contigs with zero aligned reads
parsed 1 base modification(s). Base modifictions other than 'C:m' will be counted as 'N_other'.
adding single-base motif: 'C 0'
attempting to sample 10042 reads
Using filter threshold 0.8652344 for C.
using optimized workers for A,C all-context
Done, processed 1706 rows.

This is a preview of the first 10 rows of the output bedmethyl

chr22	440	441	m	4	+	440	441	255,0,0	4	0.00	0	4	0	0	0	0	0
chr22	441	442	m	4	-	441	442	255,0,0	4	0.00	0	4	0	0	0	0	0
chr22	620	621	m	3	+	620	621	255,0,0	3	0.00	0	3	0	1	0	0	0
chr22	621	622	m	3	-	621	622	255,0,0	3	0.00	0	3	0	0	0	0	1
chr22	684	685	m	3	+	684	685	255,0,0	3	0.00	0	3	0	0	0	0	1
chr22	685	686	m	3	-	685	686	255,0,0	3	0.00	0	3	0	0	0	0	1
chr22	859	860	m	3	+	859	860	255,0,0	3	0.00	0	3	0	0	0	0	1
chr22	860	861	m	3	-	860	861	255,0,0	3	0.00	0	3	0	0	0	0	1
chr22	884	885	m	2	+	884	885	255,0,0	2	0.00	0	2	0	0	0	0	2
chr22	903	904	m	1	-	903	904	255,0,0	1	0.00	0	1	0	0	0	0	3

The input bam I'm using is a small subset of reads from human genome, and it does contained the haplotype tags.
When I checked the number of HP tags in the bam file with samtools view sorted.phased.bam | grep -c "HP:i:2" and | grep -c "HP:i:1", there are seven haplotag 1 and four haplotag 2.

I'm not sure how to troubleshoot based on the error message that I've gotten. Could you let me know what is causing the bug here?

Thanks again in advance!

Metadata

Metadata

Assignees

No one assigned

    Labels

    bugSomething isn't workingbuild-availablecustom build produced for fix.

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions