How to fetch sequence from chromosome coordinates

[kumardeep]

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 I have huge file containing chr1 to chr21 and chr X and chrY. Is there any platform where I can give the coordinate file and that give me the desired nucleotide sequences ?
Is there any difference between using hg18 human reference genome or hg19 for the same purpose. Because the experiment has used hg18 assembly to generate the coordinates (which they have given in supplementary material).

Thanks in advance.
KD

[marloch]

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 I have huge file containing chr1 to chr21 and chr X and chrY. Is there any platform where I can give the coordinate file and that give me the desired nucleotide sequences ?
Is there any difference between using hg18 human reference genome or hg19 for the same purpose. Because the experiment has used hg18 assembly to generate the coordinates (which they have given in supplementary material).

Thanks in advance.
KD

Hello,

you can use the faidx command of samtools.
First you need to create an indexed version of your FASTA:

samtools faidx genome.fa

Then you can retrieve any region, using the syntax chromosome:start-stop

samtools faidx genome.fa chr18:5233922-5235189

Cheers,
Marco

[kumardeep]

Thanks Marloch...

[JQL]

I am also trying to retrieve sequences from chromosome coordinates, but I have 50, not just one. I hope I don't have to run samtools 50 times.

Anybody has an easy way doing this?

I tried in R using system(), but there are some problems.
if I run >system('samtools faidx genome.fa chr2:1000:2000'), it prints out the sequence I wanted in console, but if I put in a variable x<-"chr2:1000:2000", it only prints out "x", but no sequences. So how to get the system() to recognize a varialbe?

Also, can't seems to get the output to a file.
>system('samtools faidx genome.fa chr2:1000:2000 >out.txt') doesn't work

thanks in advance.

[krespim]

I am also trying to retrieve sequences from chromosome coordinates, but I have 50, not just one. I hope I don't have to run samtools 50 times.

Anybody has an easy way doing this?

I tried in R using system(), but there are some problems.
if I run >system('samtools faidx genome.fa chr2:1000:2000'), it prints out the sequence I wanted in console, but if I put in a variable x<-"chr2:1000:2000", it only prints out "x", but no sequences. So how to get the system() to recognize a varialbe?

Also, can't seems to get the output to a file.
>system('samtools faidx genome.fa chr2:1000:2000 >out.txt') doesn't work

thanks in advance.

You can use the Biostring package I think. Check out this tutorial (end of page).

[crazyhottommy]

it is an old post, I have an old post addressing this http://crazyhottommy.blogspot.com/20...inates-to.html

[maximilianh]

UCSC version

You can use a command like this:

twoBitToFa http://hgdownload.cse.ucsc.edu/gbdb/hg19/hg19.2bit -bed=input.bed test.fa

or for a single region:

twoBitToFa http://hgdownload.cse.ucsc.edu/gbdb/hg19/hg19.2bit test.fa -seq=chr21 -start=1 -end=10000

Requires the UCSC tool twoBitToFa, available from http://hgdownload.cse.ucsc.edu/admin/exe/