I have sequenced a set of bacterial clones, derived from the same reference strain, and identified a few thousand mutations. Sometimes I have more than one mutation in the same gene and I would like to know whether those cases are expected or not, given a random distribution of mutations.

I was thinking that I could do this by simulating genomes of my reference strain with mutations following the same mutation spectrum as the one I have in my vcf file. Does anyone know of a tool that can achieve this?