Biopython is a great tool for interacting with biological databases. I use it to retrieve records from NCBI’s Entrez databases including Pubmed.
Unfortunately – one notable database biopython has trouble working with is the SNP database. This is due to the Bio.Entrez
parser being unable to handle the XML returned from this database. One solution is to use a built in Python XML parser, but I thought I’d try to come up with an easier solution.
To solve this problem – I wrote a function for retrieving SNP data, and parsing it into an array.
from pprint import pprint as pp
from Bio import Entrez
Entrez.email = "YOUR@EMAIL.HERE"
def pull_line(var_set, line):
"""
This function parses data from lines in one of three ways:
1.) Pulls variables out of a particular line when defined as "variablename=[value]" - uses a string to find the variable.
2.) Pulls variables based on a set position within a line [splits the line by '|']
3.) Defines variables that can be identified based on a limited possible set of values - [categorical variable, specified using an array]
"""
line_set = {}
for k,v in var_set.items():
if type(v) == str:
try:
line_set[k] = [x for x in line if x.startswith(v)][0].replace(v,'')
except:
pass
elif type(v) == int:
try:
line_set[k] = line[v]
except:
pass
else:
try:
line_set[k] = [x for x in line if x in v][0]
except:
pass
return line_set
def pull_vars(var_set,line_start,line,multi=False):
"""
Delegates and compiles data from entrez flat files dependent on whether
the type of data trying to be pulled is contained in unique vs. non-unique lines.
For example - the first line of the flat file is always something like this:
rs12009 | Homo Sapiens | 9606 | etc.
This line is unique (refers to RefSnp identifier)- and only occurs once in each flat file. On the other hand, lines
beginning with "ss[number]" refer to 'submitted snp' numbers and can appear multiple times.
"""
lineset = [x.split(' | ') for x in line if x.startswith(line_start)]
if len(lineset) == 0:
return
# If the same line exists multiple times - place results into an array
if multi == True:
pulled_vars = []
for line in lineset:
# Pull date in from line and append
pulled_vars.append(pull_line(var_set,line))
return pulled_vars
else:
# Else if the line is always unique, output single dictionary
line = lineset[0]
pulled_vars = {}
return pull_line(var_set,line)
def get_snp(q):
"""
Takes as input an array of snp identifiers and returns
a parsed dictionary of their data from Entrez.
"""
response = Entrez.efetch(db='SNP', id=','.join(q), rettype='flt', retmode='flt').read()
r = {} # Return dictionary variable
# Parse flat file response
for snp_info in filter(None,response.split('\n\n')):
print snp_info
# Parse the First Line. Details of rs flat files available here:
# ftp://ftp.ncbi.nlm.nih.gov/snp/specs/00readme.txt
snp = snp_info.split('\n')
# Parse the 'rs' line:
rsId = snp[0].split(" | ")[0]
r[rsId] = {}
# rs vars
rs_vars = {"organism":1,
"taxId":2,
"snpClass":3,
"genotype":"genotype=",
"rsLinkout":"submitterlink=",
"date":"updated "}
# rs vars
ss_vars = {"ssId":0,
"handle":1,
"locSnpId":2,
"orient":"orient=",
"exemplar":"ss_pick=",
}
# SNP line variables:
SNP_vars = {"observed":"alleles=",
"value":"het=",
"stdError":"se(het)=",
"validated":"validated=",
"validProbMin":"min_prob=",
"validProbMax":"max_prob=",
"validation":"suspect=",
"AlleleOrigin":['unknown',
'germline',
'somatic',
'inherited',
'paternal',
'maternal',
'de-novo',
'bipaternal',
'unipaternal',
'not-tested',
'tested-inconclusive'],
"snpType":['notwithdrawn',
'artifact',
'gene-duplication',
'duplicate-submission',
'notspecified',
'ambiguous-location;',
'low-map-quality']}
# CLINSIG line variables:
CLINSIG_vars = {"ClinicalSignificance":['probable-pathogenic','pathogenic','other']}
# GMAF line variables
GMAF_vars = {"allele":"allele=",
"sampleSize":"count=",
"freq":"MAF="}
# CTG line variables
CTG_vars = {"groupLabel":"assembly=",
"chromosome":"chr=",
"physmapInt":"chr-pos=",
"asnFrom":"ctg-start=",
"asnTo":"ctg-end=",
"loctype":"loctype=",
"orient":"orient="}
# LOC line variables
LOC_vars = {"symbol":1,
"geneId":"locus_id=",
"fxnClass":"fxn-class=",
"allele":"allele=",
"readingFrame":"frame=",
"residue":"residue=",
"aaPosition":"aa_position="}
# LOC line variables
SEQ_vars = {"gi":1,
"source":"source-db=",
"asnFrom":"seq-pos=",
"orient":"orient="}
# Pull out variable information:
r[rsId]['rs'] = pull_vars(rs_vars,"rs",snp)
r[rsId]['ss'] = pull_vars(ss_vars,"ss",snp,True)
r[rsId]['SNP'] = pull_vars(SNP_vars,"SNP",snp)
r[rsId]['CLINSIG'] = pull_vars(CLINSIG_vars,"CLINSIG",snp)
r[rsId]['GMAF'] = pull_vars(GMAF_vars,"GMAF",snp)
r[rsId]['CTG'] = pull_vars(CTG_vars,"CTG",snp,True)
r[rsId]['LOC'] = pull_vars(LOC_vars,"LOC",snp,True)
r[rsId]['SEQ'] = pull_vars(SEQ_vars,"SEQ",snp,True)
return r
snp = get_snp(["12009"])
print pp(snp)