I'm currently adding a few SNPs randomly into a FASTA within python using BioPython. In the following example from BioPython, I add an SNV at location "5"
http://biopython.org/DIST/docs/api/Bio.Seq.MutableSeq-class.html
from Bio.Seq import MutableSeq
from Bio.Alphabet import generic_dna
my_seq = MutableSeq("ACTCGTCGTCG", generic_dna)
my_seq[5] = "A"
print(my_seq)
## outputs ACTCGACGTCG
In my case, I'm using the reference hg38 as the input FASTA to manipulate.
After randomly inserting these SNPs, I would like to find a region whereby there are no SNPs within (let's say) +/- 50Kb, or the next nearest thing.
Is it possible to align the FASTA with variants against the reference within python? Otherwise, is there a quick way to check the variant FASTA against the reference FASTA, finding a region with no variants? What would be an efficient way to execute this task?