From genomic coordinate to Reference base

Question

[Nicola Casiraghi]

Hi, I have an extensive set of genomic coordinates (i.e. chr1:9071988) and I need to automatically get the reference (Human Genome hg19 ) base corresponding to that positions.

Could you provide me some hints to solve this issue please?

Thank you in advance!

[Nicola Casiraghi]

Thank you everybody for the useful comments and suggestions, really appreciated! samtools faidx and bedtools getfasta work both perfectly for my purpose. Many thanks again!

Answer 1

[Devon Ryan]

samtools faidx and some shell scripting should work.

[Lars]

What? Sorry, I don't get it, again. samtools? Why samtools? That makes no sense! ;)

[Devon Ryan]

Samtools can be used to extract subsequences from a fasta file. That's what the faidx command does.

[Lars]

Wow! Didn't know that! That is a nice way of doing that! Thank you for that information!

Just looked it up:

samtools faidx hg19.fasta chr1:9071988,9071988

[Devon Ryan]

You may also find the bedtools getfasta command useful.

Answer 2

[Alex Reynolds]

You could easily script the following DAS lookup for each of your coordinates:

$ wget -qO- http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr1:9071988,9071988 | grep -v '^<'
c

Answer 3

[David Langenberger]

You can try fastacmd.

Create DB (run this only once):

​formatdb -i hg19.fa -o T -p F -V

Get nucleotide:

​fastacmd -d hg19.fa -L 9071988,9071989 -s "chr1" (-S 2 if neg. Strang)