This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Batch Fetching Fasta Sequences From Bed File

Hi,

I have a number of bed files, for which I would like to extract the corresponding fasta sequences. What is the best way to do this in batch (i.e. without using e.g. the UCSC table browser), since UCSC does not provide (to my knowledge) any web service.. any Bioperl solution ? Or some other web service ? If possible, I would like to avoid installing the genome locally (otherwise the bedtools would be a solution).

Thx!

fasta bed
• 22,900 views • link
updated by Maximilian Haeussler Maximilian Haeussler • written by Carl • Carl

How many sequences are you talking? a dozen? thousands? millions? The solution may very well differ depending on your needs.

• 0 views • link
written by Chris Miller Chris Miller

make sure there is a reason for avoiding to download the genome - I've never seen an analysis that did not need to be repeated many times over and that can quickly backfire and you'll end up with far slower solutions

• 0 views • link
written by Istvan Albert Istvan Albert

4 answers

You can do this with bedtools if you download the sequence data here's a shell script for the entire process if you're using hg19:

BED=your.bed
for chr in `seq 1 22` X Y
do
    wget -O - -q http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr1.fa.gz | gunzip -c >> hg19.fa
done

fastaFromBed -fi hg19.fa -bed $BED -fo $BED.fasta

EDIT: argh. Somehow I missed that you explicitly said you didn't want to download the sequence data or bedtools would be a solution.

• 128 views • link
written by brentp brentp

".../chr1.fa.gz ..." produced a fasta file with the chr1 sequence concatenated 24 times, I needed to put the variable name in place of the 1 ".../chr1$chr.fa.gz..." to get this to work:

BED=$1
for chr in `seq 1 22` X Y
do
    wget -O - -q http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr$chr.fa.gz | gunzip -c >> hg19.fa
done

fastaFromBed -fi hg19.fa -bed $BED -fo $BED.fasta

(NB - "BED=$1" - bedfile name is first argument passed to script)

• 0 views • link
written by rjactonspsfcf rjactonspsfcf

I have a solution for hg18/b36, but not for hg19/b37.

awk '{print $1":"($2+1)"-"$3}' in.bed | xargs samtools faidx ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/pilot_data/technical/reference/human_b36_male.fa.gz > out.fa

For razip compressed and faidx indexed fasta files, samtools is able to retrieve subsequences over network.

In addition, for a couple of regions, you can still use UCSC/Ensembl URL to retrieve sequences. For many regions, I believe Ensembl has APIs, though I have not used that.

• 129 views • link
written by lh3 lh3

so all that's required for hg19 is an razip'ed file on a public server, yes?

• 0 views • link
written by brentp brentp

And it is accompanied by a ".fai" file, the fasta index: two files on a public HTTP/FTP server.

• 0 views • link
written by lh3 lh3

Yes, there is a web service (DAS) to fetch the genomic DNA sequence: e.g.:

http://genome.ucsc.edu/cgi-bin/das/hg19/dna?segment=chr1:10001000,10003000

see this other post for some solutions related to your problem: Perl To Retrieve Sequences From Ucsc Genome Browser

• 157 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

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/

If you're not on the hg19 genome, you have to index your .fa file first with faToTwoBit

• 0 views • link
updated by Ram Ram • written by Maximilian Haeussler Maximilian Haeussler

Log in to answer this question.

More posts like this