This script takes as input a list of SNP (see first line for format) and get the flanking region and output it as fasta file. Dependencies : Bioconductor : BSgenome and BSgenome.Hsapiens.UCSC.hg18 (replace with the genome you want to use) to have an idea

SNP-Flank.R

#rs10012775	chr4	122896328	0.45819	T	C	0.358

#biocLite("BSgenome.Hsapiens.UCSC.hg18")
library("BSgenome.Hsapiens.UCSC.hg18")

con <- file("SNPS_FILE")
Lines <- readLines(con)
N <- length(Lines)
for (i in 1:N){

str <- strsplit(Lines[i],"\t")

chr <- unlist(str)[2]
snp <- unlist(str)[1]
position <- as.numeric(unlist(str)[3])
if (unlist(str)[3]=="N/A") next

major <- unlist(str)[5]
minor <-  unlist(str)[6]
offset <- 25


nameseqmaj <- paste(">",snp,"-",chr,"-","MAJOR","-",major,sep='')
seqMaj <- paste(getSeq(Hsapiens,chr,position-offset,position-1),
             major,
             getSeq(Hsapiens,chr,position+1,position+offset),
             sep='')

nameseqmin <- paste(">",snp,"-",chr,"-","MINOR","-",minor,sep='')
seqMin <- paste(getSeq(Hsapiens,chr,position-offset,position-1),
             minor,
             getSeq(Hsapiens,chr,position+1,position+offset),
             sep='')

write(nameseqmaj,"FlankSnpHg18",append=T)
write(seqMaj,"FlankSnpHg18",append=T)
write(nameseqmin,"FlankSnpHg18",append=T)
write(seqMin,"FlankSnpHg18",append=T)

}

close(con)