Extract flanking sequence of SNP call

[JdeBruin]

Hi all.

I have mapped sequencing data to a reference using BWA and called SNP using samtools mpileup and generated a vcf using bcftools. What I would like to achieve is to extract the SNP along with a 35base flanking sequences. I would appreciate any help with this.

Thanks

[westerman]

bedtools always comes to mind when asking questions like "how to extract sequence". It should be easy to convert the VCF single coordinate into the start & end coordinates for a BED file (Excel or other spreadsheet program or a simple Perl/Python script). Then use fastaFromBed.

[swbarnes2]

samtools faidx, when given genomic coordinates, will make a fasta of the desired region. That another way, using software you've already got.