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"
from Bio.Seq import MutableSeq from Bio.Alphabet import generic_dna my_seq = MutableSeq("ACTCGTCGTCG", generic_dna) my_seq = "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?