I have bam files for 10 individuals mapped onto a reference genome. I'd like to extract SNPs and their flanking regions (50 bp upstream and 50 bp downstream) based on a coverage threshold. F.e. minimum coverage of 20 for all 10 bam files combined. I'd also prefer to have not more than 1 polymorhpic site in the flanking regions.
Any suggestions on how to approach this? I assume the first step would be to merge the bam files?