I have `VCF` file (containing diallelic variants) and reference genome in `Fasta` of some non-model plant. I'd like to extract `SNPs` flanking sequences. I've found that `bedtools` and `samtools faidx` could be to some extent useful, but apparently don't solve the issue. I need to get output (`fasta` or tabular file) with following SNPs representation :
TCTCTGCCAATCACTAGAGGCCGCTTTCGCTTTTA[A/G]TTTGTGTGTGGTCAGAGTTCTTCCGGACTTT
Here a little python solution:

    import pysam
    
    
    # open vcf file
    vcf = pysam.VariantFile("input.vcf")
    # open fasta file
    genome = pysam.FastaFile("genome.fa")
    # define by how many bases the variant should be flanked
    flank = 50
    
    # iterate over each variant
    for record in vcf:
        # extract sequence
        #
        # The start position is calculated by subtract the number of bases
        # given by 'flank' from the variant position. The position in the vcf file
        # is 1-based. pysam's fetch() expected 0-base coordinate. That's why we
        # need to subtract on more base.
        #
        # The end position is calculated by adding the number of bases
        # given by 'flank' to the variant position. We also need to add the length
        # of the REF value and subtract again 1 due to the 0-based/1-based thing.
        #
        # Now we have the complete sequence like this:
        # [number of bases given by flank]+REF+[number of bases given by flank]
        seq = genome.fetch(record.chrom, record.pos-1-flank, record.pos-1+len(record.ref)+flank)
    
        # print out tab seperated columns:
        # CRHOM, POS, REF, ALT, flanking sequencing with variant given in the format '[REF/ALT]'
        print(
            record.chrom,
            record.pos,
            record.ref,
            record.alts[0],
            '{}[{}/{}]{}'.format(seq[:flank], record.ref, record.alts[0], seq[flank+len(record.ref):]),
            sep="\t"
        )

fin swimmer

**EDIT**:

I rewrote the code. The version before had an off-by-one error.