From genomic coordinate to Reference base

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!

genome

• 4,185 views • link

updated 4.4 years ago by Ram

• written 11.5 years ago by 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!

• 0 views • link

updated 4.4 years ago by Ram

• written 11.5 years ago by Nicola Casiraghi

3 answers

samtools faidx and some shell scripting should work.

• 0 views • link

written 11.5 years ago by Devon Ryan

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

• 0 views • link

written 11.5 years ago by Lars

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

• 0 views • link

written 11.5 years ago by Devon Ryan

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

• 0 views • link

updated 4.4 years ago by Ram

• written 11.5 years ago by Lars

You may also find the bedtools getfasta command useful.

• 0 views • link

updated 4.4 years ago by Ram

• written 11.5 years ago by Devon Ryan

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

• 0 views • link

updated 4.4 years ago by Ram

• written 11.5 years ago by Alex Reynolds

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)

• 0 views • link

updated 4.4 years ago by Ram

• written 11.5 years ago by David Langenberger

More posts like this

Codis marker exact genomic coordinates
written 6.9 years ago by nmmsv •

I'm looking for a complete list of genomic positions for codis markers. I find a lot of resources online about the genomic region ( i.e., …
Retrieve PTM (protein post translational modifications) for a list of protein coordinates (specific…
written 7.1 years ago by Lalla •

Dear all I am analyzing a set of alternative spliced, protein-coding, exons. I have been able to map the exons coordinates to the corresponding protein …
Calculating GC content of all CpG islands
written 7.3 years ago by Baby Bioinformatician

Hi Biostars, I have two files containing the genomic ranges of all CpG islands and promoter sequences of Human genome hg19. I would like to …
Aligning short reads to a human specific gene instead of reference genome
written 8.3 years ago by seta

Hi all, My focus of study is on the two human homologous genes related to blood groups, which have also hybrid forms due to gene-rearrangement. …
how to get corresponding transcript/gene sequence from aligned reads to human genome?
written 9.5 years ago by seta

Hi all, I have two set Illumina reads from mRNA sequencing, control and treatment. I separately mapped them to the human reference genome and got …
The importance of LB field in @RG tag for BreakDancer analysis
written 11.8 years ago by Nicola Casiraghi

I have WGS paired-end data and I want to use BreakDancer to call structural variants. Both my case and controls bam files do not have, …
BICseq R package problem
written 12.1 years ago by Nicola Casiraghi

Hi all, I am interested in using the tool BICseq for a CNV detection analysis. Unfortunately, I am writing you because I cannot solve a …
problems in hg19 and b37 compatibility
written 12.2 years ago by Nicola Casiraghi

Hi everybody, A bam file has been aligned using hg19 reference genome. Thus, the chromosome notation is [chrM, chr1, chr2, chr3, chr4, ..., chrX,chrY]. I …
How To Fetch Sequence From Chromosome Coordinates
written 13.9 years ago by monukmr98

Hi all I have chromosome coordinates (e.g. 5233922 to 5235189 of chr18). how can I get the sequence corresponding to this region ? Like this …
Genomic Change To Aa Change?
written 15.4 years ago by Biomed

I have a set of genomic coordinates and genomic sequence changes (i.e g.Agenomicposition_C ). I need to find in which transcripts and what corresponding amino …