Counts are sometimes higher than coverage (circRNA specific?)

[Raif-fl]

In Marine's output, it is sometimes possible for the counts to be greater than the coverage. To see an example of this behaviour, we can run marine on a small example dataset using the following commands:

module load marine

python /tscc/projects/ps-yeolab3/kflanagan/MARINE-1.0.9-alpha/marine.py --bam_filepath /tscc/lustre/ddn/scratch/kflanagan/Tau_marine_run/one_circ_files/one_circ.bam --output_folder ./results --strandedness 0 --cores 16 --contigs "chr1:145631664|145706593" --bedgraphs "CT" --paired_end --num_intervals_per_contig 2 --keep_intermediate_files

If we look at the final_filtered_site_info.tsv output file, we will see that the edit site at position 74928
has a count of 16 while the coverage is only 11.

This glitch is rare, happening only one other time in the example at position 15509, and does not seem to have a significant effect on downstream analysis. I am unsure of what causes this bug, but this result should be impossible.

I will note that this bug was encountered as part of an ongoing effort to apply MARINE to circular RNA datasets that have been passed through the ciriQUANT pseudo-reference alignment pipeline, and thus may be specific to circular RNA examples.

[byee4]

@howardxu520 discovered a similar issue where counts are greater than coverage. We looked at a few sites with single edit coverage (=1 conversion) and zero coverage. These sites were excluded due to samtools depth excluding softclipped read positions. We're not sure whether or not edit calling counts softclipped positions, it probably shouldn't?

[ekofman]

Yeah probably shouldn't... if you could make a minimal example (for example, a bam containing just one or two reads for which the edit is being erroneously counted), that would be a huge help to me as I check what's going on and make a new integration test to cover that scenario.

[byee4]

Sounds good, we'll get on it!

[byee4]

I've generated a small BAM file on TSCC but can also upload the inputs to here if it's more convenient. Also it looks like I might be wrong, the reads here don't seem to be softclipped:

/tscc/projects/ps-yeolab3/bay001/for_erick/scratch/run_marine.sh

final_filtered_site_info_annotated.txt

[ekofman]

thank you!

[ekofman]

it looks like the bed locations at which coverage is being calculated may be off by one, landing into a soft-clipped region, which I thought I'd addressed at some point but I guess would rarely come up as an issue unless there was an edit at the very end of a read. I'll see if I can fix this in the next few days.