calculate denominator for all sites and list homopolymer distance.
the denominator is the total count of overlapping bases in each bq_bin we need to report this so users can get an error rate by location in the bed file.
we also want to report distance to homopolymer where 0 means it's in, -n means it's before the read has passed through the homopolymer and +n means it's after the read has passed through the homopolymer. This will depend on the read orientation.
this will be a separate command:
fraguracy denominator \
--min-homopolymer-size 3 \ # must be 2 or more
--max-homopolymer-distance 25 \
--fasta $fasta \
--min-mapping-quality 10 \
--output-prefix $sample \ # will write $sample.fraguracy.denominator.bed.gz
$bam_path
$sample.fraguracy.denominator.bed.gz will look like:
#chrom start end bq-bin distance-to-hp denominator
chr1 0 123 0-10 >25 135
chr1 130 131 0-10 25 2
chr1 131 132 0-10 24 1
...
adding zip here for Jason to test fraguracy_dev.zip
Jason pointed out that it's better to have an extra partition of homopolymer distance (before and after). This will make counts.txt very sparse, but we can use pandas or some other group-by operation to sum all read position or sum all homopolymer distances to get the desired plots.
also note that a reference homopolymer of length $m may only partly overlap a read, if it overlaps less than min-homopolymer-bases of the read then we would not count it as a homopolymer.