Get Snp Position From A Python Interface.

Question

[Peixe]

Hi,

I have a Python script which takes as input a list of SNPs, then calls tabix to retrieve .VCF file and does some stuff with it. The problem is that i need to manually add (using UCSC genome browser web*) the coordinates of the SNP for tabix to work properly.

Does any one know any simple python API / function / module to easy add this information without having to perform the search manually?

Thanks!

Answer 1

[Pierre Lindenbaum]

You could use an external database to find the coordinates of your rs## (not tested, please check the +0/+1 for the coordinates)

for R in rs25 rs26 rs27 
do
    mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -Dhg19 -N -e "select concat(chrom,':',chromStart+1,'-',chromEnd) from snp135 where name='$R';" |\
    xargs tabix my.vcf.gz
done

[Peixe]

Yes, this fits perfect! Coordinates are right. Thanks, @Pierre Lindenbaum!

[Shicheng Guo]

ERROR 2013 (HY000): Lost connection to MySQL server at 'reading initial communication packet', system error: 104
open: No such file or directory
[tabix] was bgzip used to compress this file? my.vcf.gz

[_quantum_girl_]

Pierre Lindenbaum why "from snp135"?

[Pierre Lindenbaum]

because 11 years ago that was the name of the snpTable.

now for hg38: https://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/snp151.sql

Answer 2

[brentp]

You can do this with cruzdb

from cruzdb import Genome
import sys

fname = sys.argv[1]
hg19 = Genome('hg19')
snp135 = hg19.snp135

for rs in (l.split()[0].rstrip() for l in open(fname) if l.startswith("rs")):
    print snp135.filter_by(name=rs).first()

which takes a file with the rs id's in the first columns and prints out a bed file with the chromosome locations from dbsnp 135.

[Peixe]

Wow! This is amazing too! I didn't know about the existence of this module. I think i will really get much profit out of it in the future. Thanks, @brentp!

[epigene]

hi @brentp, I run into your tool and noted that the following line is giving an error.

snp135 = hg19.snp135

I realized that it's because snp135 is no longer the latest database. How do you determine which version of dbsnp to use without knowing what's currently available? PS: I tried snp141 and it worked..

[forever]

When I try to run it, I got the below error:

Traceback (most recent call last):
  File "hellowworld.py", line 1, in <module>
    import Genome
ImportError: No module named 'Genome'

What shall I do?

[WouterDeCoster]

Have you installed cruzdb first?

[forever]

Yes I installed it and now get another error:

Traceback (most recent call last):
  File "hellowworld.py", line 1, in <module>
    from cruzdb import Genome
  File "C:\Python34\lib\site-packages\cruzdb\__init__.py", line 5, in <module>
    from . import soup
  File "C:\Python34\lib\site-packages\cruzdb\soup.py", line 1, in <module>
    from . import sqlsoup
  File "C:\Python34\lib\site-packages\cruzdb\sqlsoup.py", line 458
    except KeyError, ke:
                   ^
SyntaxError: invalid syntax

[WouterDeCoster]

I guess this synthax error is because you use Python3.4 and the module was written for python2, but that's just my first guess.

[forever]

When I run it, I got

fname = sys.argv[1]
IndexError: list index out of range

What is the problem?

[forever]

I got

fname = sys.argv[1]
IndexError: list index out of range

[WouterDeCoster]

You need to save the code as a script and use the input file as a argument when executing the script using your python interpreter.