getflank.r

# Load the BSgenome package and download the hg19 reference sequence 
# For documentation see http://www.bioconductor.org/packages/release/bioc/html/BSgenome.html
source("http://www.bioconductor.org/biocLite.R")
biocLite("BSgenome")
biocLite("BSgenome.Hsapiens.UCSC.hg19") #installs the human genome (~850 MB download).

library('BSgenome.Hsapiens.UCSC.hg19')

# Function to pull flanking sequence. Defaults to +/- 10 bp
getflank <- function(position, alleles="[N/N]", chr="chr12", offset=10) {
    leftflank  <- getSeq(Hsapiens,chr,position-offset,position-1)
    rightflank <- getSeq(Hsapiens,chr,position+1,position+offset)
    paste(leftflank,alleles,rightflank,sep="")
}

# Try this as an example. should return "TGGCAACTCC[C/A]TTCCATTTGC", 
#   which is exactly what you'd see if you searched db SNP: 
#   http://www.ncbi.nlm.nih.gov/sites/entrez?db=snp&cmd=search&term=rs1520218
getflank(102741896, alleles="[C/A]")