All avgoutcoverage = 0

[mrmilhaven]

Hello,

I am trying to mutate a rather large genome, here is the command:
python ~/bamsurgeon/bin/addsnv.py -v input_mutations.bed -f realvariants_Ayumu.sorted.golden.newheader.bam -r panTro6.autosomal.fa -o realvariants_Ayumu_genome0_1.newheader.sorted.golden.newheader.mutated.bam

It run's without an issue until the end in which I get this error:
ERROR 2021-12-10 12:39:50,440 no succesful mutations

Looking at each individual mutation being generated, they all have this issue:

INFO 2021-12-10 12:39:48,728 haplo_chr7_43034975_43034975 avgincover: 102.666667, avgoutcover: 0.000000
WARNING 2021-12-10 12:39:48,729 haplo_chr7_43034975_43034975 dropped for outcover/incover < 0.9

Essentially all avgoutcover is zero and they get dropped. The average coverage across these genomes is 100X, so I am confused why there are so little output reads. I have also had an issue with the organization of the reference genome and the BAM file, and I have since fixed that issue by reorganizing the header of the BAM using samtools reheader.

Best,
Mark

[mrmilhaven]

Hello,

Just following up on this issue. Is there any headway?

[adamewing]

Sorry - difficult to get to bamsurgeon stuff lately. When I've seen this in the past it's usually due to the realignment step failing somewhere. If the original .bam was created via (e.g.) bwa mem, make sure the same aligner is used for realignment (--aligner mem in the example case). The default, which is what you're using based on the command line args above, is bwa aln which is super outdated at this point (as is bamsurgeon in general, IMO).

[mrmilhaven]

Thank you! It is starting to mutate successfully. However, some mutations are being skipped because is states there are no reads in the region. But when I look at those specific sites, they have adequate coverage (i.e. 95X). Do you know why this could happen?

Results:
INFO 2022-01-21 10:47:24,918 haplo_chr1_177131995_177131995 creating tmp bam: addsnv.tmp/haplo_chr1_177131995_177131995.tmpbam.a43d8ea8-cd4f-4590-80d7-8314ae
e35aea.bam
INFO 2022-01-21 10:47:53,631 haplo_chr1_177131995_177131995 len(readlist): 0
WARNING 2022-01-21 10:47:53,633 haplo_chr1_177131995_177131995 too few reads in region (0) skipping...

samtools depth -a

chr1 177131994 96
chr1 177131995 95
chr1 177131996 92

[mrmilhaven]

Hello! Any updates on this issue?