Background: The CFTR gene has a Poly TG/Poly T region in intron 8 that has [varying clinical][1] consequences. 
I'm using targeted NGS for this region (150bp, paried end) and I'm trying to call the PolyTG/PolyT genotypes from fastq files. My first attempt was bbmap run with default parameters and a reference including (TG)9-12(T)5,7,9:

    bbmap.sh ref=$ref in1=$R1 in2=$R2 scafstats=$out

39/40 samples could be correctly called based on the bbmap stats (header shown below). Heterozygous calls were made from ref_name's with the highest unambiguous reads. Homozygotes could be figured based on the ratio of the top two. 
  
      #sample	ref_name	per_unambiguousReads	unambiguousMB	per_ambiguousReads	ambiguousMB	unambiguousReads	ambiguousReads	assignedReads	assignedBases


I expanded the reference to catch rare genotypes (e.g. T6 and TG13). The expanded reference has (TG)8-13(T)2-10. Running bbmap with this reference results in many incorrect calls.

I'm wondering if there are parameters I should be tweaking in bbmap or if there is a better solution to calling the Poly TG/Poly T genotypes that I'm not thinking of.

Thanks,
Bob

  [1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2699602/
Just to wrap this up. I found this [paper][1]. I contacted Dr. Pagin and he explained their script used grep on fastqs and counted different alleles to make a call per sample. I used that idea in this [script][2]. If anyone uses this script in the future, the HET_CUTOFF (line 554) will likely vary depending on the panel you use to sequence CFTR; adjust according to your data.


  [1]: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0149426
  [2]: https://github.com/rjsicko/scripts/blob/master/polyTG_table_v1.0.sh