How to get protein change HGVS data for a given dbSNP RS ID, using NCBI EUtils

protein_hgvs_for_snp_id.py

"""Python 3 script to get partial protein HGVS given NCBI dbSNP ID

Example:

$ python3 protein_hgvs_for_snp_id.py 334
snp_id:
334
gene:
{'name': 'HBB', 'gene_id': '3043'}
protein_change_hgvs:
Glu7Ala

"""
from urllib import request
import json
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("snp_id", metavar="I", type=str,
                    help="dbSNP ID, e.g. 6025 (for rs6025)",
                    nargs="?",  # Makes this argument optional
                    default="6025")

snp_id = parser.parse_args().snp_id

eutils_url = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?retmode=json&db=snp&id=" + snp_id
f = request.urlopen(eutils_url)
content = f.read()

eutils_json = json.loads(content.decode("utf-8"))
snp_json = eutils_json["result"][snp_id]
gene = snp_json["genes"][0]

snp_docsum_list = snp_json["docsum"].split(";")

protein_change_hgvs = ""
for item in snp_docsum_list:
    # If this item contains a known RefSeq protein, e.g.
    # NP_056953.2:p.Pro12Ala,XM_011533840.1:c.-2-28078C&gt
    if item[2:5] == "NP_":
        protein_change_hgvs = item.split("p.")[1].split(",")[0]

print("snp_id:")
print(snp_id)

print("gene:")
print(gene)

print("protein_change_hgvs:")
print(protein_change_hgvs)