Difficulty recognizing BAM file?
Closed this issue · 3 comments
I was using tldr to interrogate some output BAM files from the methylation caller Megalodon, I'm using the BAM files Megalodon (a methylation caller that uses sequence alignment to improve its methylation calling, or at least that is my understanding) can give as outputs (so these files should include methylation data, which might change how the data is formatted?). I wanted to use these rather than from a typical aligner like minimap2 to preserve the methylation information on them, but was wondering if that was causing my issue. Anyways here is the error output:
Traceback (most recent call last):
File "/home1/malonema/.local/bin/tldr", line 2128, in
main(args)
File "/home1/malonema/.local/bin/tldr", line 1807, in main
bam = pysam.AlignmentFile(bams[0])
File "pysam/libcalignmentfile.pyx", line 751, in pysam.libcalignmentfile.AlignmentFile.cinit
File "pysam/libcalignmentfile.pyx", line 1000, in pysam.libcalignmentfile.AlignmentFile._open
ValueError: file has no sequences defined (mode='r') - is it SAM/BAM format? Consider opening with check_sq=False
and here is the script I'm running
tldr -b $bamList -r $reference -o $outBaseName -e none --detail_output --extend_consensus 100 --procs 32
Yes I've noticed that the .bams output by megalodon don't seem to work directly with tldr despite their being generated internally via minimap2. If you sort and index them it resolves the error but yields no output. I'll investigate further.
So I've tracked this down, the megalodon docs here: https://nanoporetech.github.io/megalodon/file_formats.html#modified-base-mapping, the "reads" in mod_mappings.bam aren't actually the nanopore reads but are sequences from the reference genome. This is a bit unfortunate as it makes mod_mappings.bam unsuitable for most downstream applications dealing with sequence variation (variant detection, phasing, etc). The other bam file, mappings.bam may be suitable for TLDR but I'm not sure as it seems to involve mappy, which to my understanding doesn't output soft-clipped reads. You're probably better off re-aligning the .fastq output by megalodon with minimap2 and using that output for TLDR.
closing due to likely resolution