How to use Bed file to extract sequence from FASTA file?

Question: How to use Bed file to extract sequence from FASTA file?

4.1 years ago by

I tried bedtools getfasta and i get the errors that chromosome was not found in fasta file but I have triple checked it there is no blank space the chromosome name in bed file is exactly the same as in fasta file. I would like to know is there any alternatives other than using bedtools getfasta in order to extract the sequence.

bed chip-seq sequence fasta software error • 5.7k views

ADD COMMENT • link •

modified 4.1 years ago by Maximilian Haeussler • 1.3k • written 4.1 years ago by allyson1115ar • 20

samtools faidx extracts subsequence from indexed reference sequences (http://samtools.sourceforge.net/samtools.shtml)

Also use search feature to look for similar posts on this forum.

ADD REPLY • link modified 4.1 years ago • written 4.1 years ago by Ashutosh Pandey ♦ 11k

I tried. But can't solve it.

ADD REPLY • link written 4.1 years ago by allyson1115ar • 20

Solved well. Thanks.

ADD REPLY • link written 4.1 years ago by allyson1115ar • 20

There is one problem in this tool. The result is not accurate. The sequence extracted is not the same as in bed file coordinates

ADD REPLY • link written 4.1 years ago by allyson1115ar • 20

It's more likely that you made an error than that samtools faidx did...

ADD REPLY • link written 4.1 years ago by Devon Ryan ♦ 90k

This is a part of the output i get from samtools faidx:

>chr1:179757197-179758470:1251-1255
TGAGT
>chr1:201237463-201238874:41-45

>chN
>chr1:201237463-201238874:62-80
238874
TGCCACAGCTGN
>chr1:201237463-201238874:62-81
238874

ADD REPLY • link written 4.1 years ago by allyson1115ar • 20

You appear to have used a heavily mistyped command (e.g., chN, also the ranges are non-sensical). Post the exact command that you used and mention where you got the genome.

ADD REPLY • link written 4.1 years ago by Devon Ryan ♦ 90k

4.1 years ago by

Matt Shirley ♦ 9.0k

Cambridge, MA

Matt Shirley ♦ 9.0k wrote:

You might try the "faidx" command from https://github.com/mdshw5/pyfaidx. There is a --default-seq parameter that allows filling in nonexistent sequence, as well as a --lazy parameter that disables bounds checking. You can pass a bed file using --bed.

ADD COMMENT • link written 4.1 years ago by Matt Shirley ♦ 9.0k

4.1 years ago by

Maximilian Haeussler • 1.3k

UCSC

Maximilian Haeussler • 1.3k wrote:

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

ADD COMMENT • link modified 4.1 years ago • written 4.1 years ago by Maximilian Haeussler • 1.3k

Please log in to add an answer.

Similar posts • Search »