How to batch save short prespecified sequences from the human genome

Question

[jamesebarrett]

I would like to save sequences from, say, chr1:100-150, chr1:300-350 etc to a file. I have several hundred such genome regions that I want to do this for using GRCh37. I would be grateful if someone could point me to a quick way of doing this.

Edit: samtools faidx ref.fa chr1:100-150 works nicely. Many thanks to answers below.

Answer 1

[Alex Reynolds]

If you have a BED file of intervals (in.bed) and a directory of FASTA for your genome assembly of interest (e.g., hg19), you could use a script like bed2faidxsta.pl to write data to a FASTA file (out.fa):

$ ./bed2faidxsta.pl --fastaDir=/path/to/hg19/fasta < in.bed > out.fa

It relies on samtools and indexed FASTA files.

Here's the script:

[jamesebarrett]

Brilliant thank you. I am currently downloading chromFa.tar.gz from https://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/. Could you please tell me how I can index these?

Answer 2

[Jean-Karim Heriche]

Use the Ensembl perl API with something like this:

use strict;
use warnings;
use Bio::EnsEMBL::Registry;


my $input = $ARGV[0];

my ($chr,$coords) = split(/:/,$input);
my ($start,$end) = split(/-/,$coords);

my $registry = "Bio::EnsEMBL::Registry";
$registry->load_registry_from_db(-host => 'ensembldb.ensembl.org',
                             -port => 3337,
                             -user => 'anonymous',
                             -passwd => undef,
                             -db_version => '99');
my $slice_adaptor = $registry->get_adaptor( 'Homo sapiens', 'Core', 'Slice' );
# Obtain a slice covering the region
my $slice = $slice_adaptor->fetch_by_region( 'chromosome', $chr, $start, $end);
print ">$chr:$start-$end\n",$slice->seq(),"\n";