Question: Finding a variant's flanking sequence
0
gravatar for angrypigeon
13 days ago by
angrypigeon80
angrypigeon80 wrote:

I am looking to improve the speed with which I find genomic sequence that flanks a given variant. Currently for a variant of interest, I use python to pull its chromosome and location, then query the UCSC genome browser for a given number of downstream and upstream bases like this:

check_output('wget -qO- http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=%s:%s,%s' % (chrom,low,high), stderr=STDOUT, shell=True)

This works just fine, especially when only working with a handful of variants. But when working with thousands of variants, this solution becomes quite slow because of the constant queries. So, is there a faster way of accomplishing the same task? I assume there is a more elegant way to use an offline copy of hg19 or something that would eliminate the need to constantly probe the genome browser?

sequencing sequence • 122 views
ADD COMMENTlink modified 13 days ago by Chris Miller19k • written 13 days ago by angrypigeon80
2
gravatar for Chris Miller
13 days ago by
Chris Miller19k
Washington University in St. Louis, MO
Chris Miller19k wrote:

The command you're looking for is samtools faidx

Example: samtools faidx hg19.fa 1:6403804-6403874

ADD COMMENTlink written 13 days ago by Chris Miller19k

This seems to be what others use, but it doesn't seem to be working for me. I do have hg19.fa.fai within the same directory as hg10.fa, but running samtools faidx hg19.fa 1:6403804-6403874 just outputs >1:6403804-6403874. Am I missing something?

ADD REPLYlink written 13 days ago by angrypigeon80
2

Check if you chromosomes' name start with 'chr' :

head hg19.fa
ADD REPLYlink written 13 days ago by michael.ante1.9k

Beautiful, samtools faidx hg19.fa chr1:6403804-6403874 works perfectly.

ADD REPLYlink written 13 days ago by angrypigeon80
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1496 users visited in the last hour