<p>The <code>exonerate</code> tool comes with some very handy utilities for sequence manipulation. Here, we are interested in <code>fastasubseq</code> and <code>fasterevcomp</code>. First, install <code>exonerate</code> as described <a href="https://www.ebi.ac.uk/about/vertebrate-genomics/software/exonerate" rel="nofollow noreferrer">here</a>. If you're using Debian or Ubuntu or one of their derivatives, it might be enough to run</p>

<pre><code>sudo apt install exonerate
</code></pre>

<p>Then, follow these steps:</p>

<ol>
<li><p>Create one fasta file for each sequence in <code>genome.fasta</code></p>

<pre><code>fastaexplode genome.fasta 
</code></pre></li>
<li><p>A little preprocessing</p>

<pre><code><span class="math-container">$ awk -F[:-] '{ 
                if($2&gt;$3){
                    start=$3-1; 
                    len=$2-start" -"
                }
                else{
                    start=$2-1; 
                    len=$3-start
                } 
                print $</span>1,start,len}' id.txt 
Contig3 7 8 -
Contig2 5 5
Contig1 3 9 -
</code></pre>

<p>So, that will print the target contig's name, the start coordinate of the target sequence, the length of the target sequence and a <code>-</code> if it is on the reverse strand. The <code>fastasubseq</code> program I will use below starts counting from 0 and not from 1, so we also need to subtract 1 from the coordinates to get them right.</p></li>
<li><p>Now, feed that into <code>fastafetch</code> and <code>fastarevcomp</code></p>

<pre><code>tmpFile=<span class="math-container">$(mktemp); 
    &gt; subseqs.fa
    awk -F[:-] '{ 
                if($2&gt;$3){
                    start=$3-1; 
                    len=$2-start" -"
                }
                else{
                    start=$2-1; 
                    len=$3-start
                } 
                print $1,start,len}' id.txt |
    while read cont start len rev; do 
        fastasubseq "$</span>cont".fa <span class="math-container">$start $</span>len &gt; <span class="math-container">$tmpFile; 
        if [[ -n $</span>rev ]]; then 
            fastarevcomp "<span class="math-container">$tmpFile" &gt;&gt; subseqs.fa; 
        else 
            cat "$</span>tmpFile" &gt;&gt; subseqs.fa; 
    fi
    done
</code></pre>

<p>This will produce the following file when run on your data:</p>

<pre><code>$ cat subseqs.fa 
&gt;Contig3:subseq(6,9) [revcomp]
GATCAGCGC
&gt;Contig2:subseq(4,6)
AGCGCT
&gt;Contig1:subseq(2,10) [revcomp]
CTTGATGGCC
</code></pre></li>
<li><p>Finally, remove the per-contig fasta files</p>

<pre><code>rm Contig*fa
</code></pre></li>
</ol>

CC BY-SA 4.0
2019-06-14T06:29:02.870Z
Thank you so much terdon for your time and help. Your script works awesome. However, I did a mistake in my question. Instead of mentioning Contig3:15-7, Contig1:12-3, I have mentioned Contig3:15-7 and Contig2:5-10.              
                                                                                                                                       I need to reverse complement the sequences, only if it in reverse order/inverse order. Which ever sequences are not in reverse order (Properly order sequences), only to be extracted as its not to be reverse complemented.
Kumar
CC BY-SA 4.0
2019-06-14T06:31:16.100Z
Please accept my sincere apologies for the same. Could you please do the needful.
Kumar
CC BY-SA 4.0
2019-06-14T08:21:04.610Z
@K.Dineshkumar I'm sorry, but I don't understand the problem. The script will only reverse complement if the start position is smaller than the end position. Isn't that what you want? In the example above, `Contig3:15-7` and `Contig1:12-3` have been reverse complimented, but `Contig2:5-10` has not.
terdon
CC BY-SA 4.0
2019-06-14T08:39:44.767Z
@K.Dineshkumar please edit your question and add this. Show us the expected output. It's very hard to read in the comments.
terdon
CC BY-SA 4.0
2019-06-14T08:49:05.093Z
@K.Dineshkumar 
 Never mind, I found the problem. The `fastasubseq` tool uses 0-based coordinates (the first position of a sequence is 0), but you have 1-based (the first position of a sequence is 1). Try the updated answer.
terdon
CC BY-SA 4.0
2019-06-14T09:32:43.913Z
@terdon Please see my edited question, wherein I have included my expected outcome.
Kumar
CC BY-SA 4.0
2019-06-14T09:38:14.107Z
@terdon, The revised/updated script printing one base less in the end of the sequence. For example,  `Contig2:5-10` says, print the sequence from `5 to 10` But the script print only `5-9`.
Kumar
CC BY-SA 4.0
2019-06-14T09:56:44.997Z
@K.Dineshkumar argh! Yes, of course, sorry. Fixed now, try the (last, I hope) update. And thanks for adding the output to your question, now I can check :)
terdon
CC BY-SA 4.0
2019-06-14T11:19:14.567Z
@terdon, Now it works perfect. I really appreciate your help. Thank you.
Kumar
CC BY-SA 4.0
2019-06-14T13:41:39.000Z
@K.Dineshkumar you're welcome. If this answers solves your issue, please take a moment to [accept it](//$$SITEURL/help/accepted-answer) by clicking on [the checkmark on the left](https://meta.stackexchange.com/a/5235/203101). That will mark the question as answered and is the way that thanks are conveyed on the Stack Exchange sites. You don't need to do it now, you might want to wait in case anyone can give a better solution, but please remember to accept an answer since that marks the question as answered.
terdon
CC BY-SA 4.0
2019-06-15T17:40:43.953Z
@K.Dineshkumar thanks for accepting! Feel free to accept one of the other answers though, I won't take it the wrong way :) They give more efficient solutions.
terdon
CC BY-SA 4.0
2019-06-13T15:00:15.163Z
2019-06-14T09:56:17.523Z
2019-06-14T09:56:17.523Z
terdon
<p>Since you tagged python, here is a solution that solves this in a single script using Biopython:</p>

<pre><code>import sys
from Bio import SeqIO

d = {}
with open('id.txt') as f:
    for line in f:
        contig, coords = line.split(':')
        start, end = map(int, coords.split('-'))
        d[contig] = start, end

for record in SeqIO.parse('genome.fasta', 'fasta'):
    start, end = d[record.id]
    if start &gt; end:
        record.seq = record.seq[end - 1: start].reverse_complement()
    else:
        record.seq = record.seq[start - 1: end]
    SeqIO.write(record, sys.stdout, 'fasta')
</code></pre>

<p>Output on your example:</p>

<pre><code>&gt;Contig1
CTTGATGGCC
&gt;Contig2
AGCGCT
&gt;Contig3
GATCAGCGC
</code></pre>

CC BY-SA 4.0
2019-06-15T18:03:11.653Z
Thank you @Chris_rands
Kumar
CC BY-SA 4.0
2019-07-06T13:55:24.417Z
I try to use it for a bacterial genome. But it shows error as follows,          
 `Traceback (most recent call last):
  File "/home/ga/Documents/commands/scripts/amphora/reverse_comple_seq_extrac.py", line 12, in <module>
    start, end = d[record.id]
KeyError: '1'`
Kumar
CC BY-SA 4.0
2019-07-08T07:37:58.193Z
@K.Dineshkumar try using `d[record.description]` instead, i can't know what will work without seeing your data though
Chris_Rands
CC BY-SA 4.0
2019-07-08T08:40:00.920Z
My genome of interest can be obtained from the following link `https://www.ncbi.nlm.nih.gov/Traces/wgs/CAHO01?display=contigs`
Kumar
CC BY-SA 4.0
2019-07-08T08:40:26.903Z
My coordinate is as follows,                                                                                                  `CAHO01000012.1:92353-95295
CAHO01000017.1:27362-27949
CAHO01000014.1:11449-12012
CAHO01000018.1:51156-52928
CAHO01000030.1:58079-56889
CAHO01000001.1:60348-62039
CAHO01000026.1:45206-44508
`
Kumar
CC BY-SA 4.0
2019-07-08T08:57:45.717Z
When I tried with `d[record.description]` instead `of d[record.id]`, I end up the error as follows,                                                                                                                                      `Traceback (most recent call last):
  File "/home/ga/Documents/commands/scripts/amphora/reverse_comple_seq_extrac.py", line 12, in <module>
    start, end = d[record.description]
KeyError: 'CAHO01000001.1 Xanthomonas citri pv. mangiferaeindicae LMG 941 WGS project CAHO00000000 data, contig XMINcontig00001, whole genome shotgun sequence'`
Kumar
CC BY-SA 4.0
2019-06-14T14:05:20.187Z
2019-06-14T14:05:20.187Z
Chris_Rands
<p>If you're okay with using a Perl script, you could convert your headers into stranded BED6 and use the following script <code>bed2faidxsta.pl</code> with <code>samtools</code>-indexed FASTA to get back stranded sequence:</p>

<p><a href="https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640" rel="nofollow noreferrer">https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640</a></p>

<p>This script looks at the sixth column of your BED file of chromosome, start, and stop positions, to decide whether to do a reverse complement of the queried sequence. </p>

<p>You could use <code>awk</code> to quickly build such a BED6 file:</p>

<pre><code><span class="math-container">$ awk -vFS='[:-]' '($</span>0~/^&gt;/){ if (<span class="math-container">$2 &lt; $</span>3) { strand = "+"; } else if (<span class="math-container">$2 &gt; $</span>3) { strand = "-"; } else { printf("Coordinates must differ\n"); exit 1; } printf("%s\t%d\t%d\t.\t0\t%c\n", <span class="math-container">$1, $</span>2, $3, strand; }' in.fa &gt; out.bed6
</code></pre>

<p>To set up indexed FASTA, here's a generic shell-based approach, which uses [Kent utilities][1] <code>fetchChromSizes</code> to get a list of chromosomes for the build <code>mm9</code>:</p>

<pre><code><span class="math-container">$ mkdir mm9
$</span> cd mm9
<span class="math-container">$ while read chr; do echo "$</span>chr"; wget -qO- "http://hgdownload.cse.ucsc.edu/goldenpath/mm9/chromosomes/<span class="math-container">$chr.fa.gz" | gunzip -c - &gt; $</span>chr.fa; samtools faidx $chr.fa; done &lt; &lt;(fetchChromSizes mm9 | cut -f1)
</code></pre>

<p>This can be adjusted for <code>mm10</code>, <code>hg19</code>, <code>hg38</code>, and so on.</p>

<p>Once you have indexed FASTA and the above Perl script, you could do something as easy as the following:</p>

<pre><code>$ ./bed2faidxsta.pl --fastaDir="/path/to/indexed/mm9" &lt; in.bed &gt; out.fa
</code></pre>

<p>If you're working with a reference genome, the nice thing about this approach is that <code>awk</code>, Perl, and indexed lookups are fast. You can reuse indexed FASTA for other projects that also require sequence lookups by genomic interval.</p>

CC BY-SA 4.0
2019-06-15T18:03:29.453Z
Thank you @Alex Reynolds
Kumar
CC BY-SA 4.0
2019-06-14T18:41:01.200Z
2019-06-14T18:41:01.200Z
Alex Reynolds
<p>I have the series of coordinates in <code>id.txt</code> file, whose coordinates sequences are in <code>genome.fasta</code> file.
The coordinates of id.txt file is shown below,</p>

<pre><code>Contig3:15-7
Contig2:5-10
Contig1:12-3
</code></pre>

<p>The genome.fasta file is shown below,</p>

<pre><code>&gt;Contig1
AAGGCCATCAAGGACGTGGATGAGGTCGTCAAG
&gt;Contig2
ACGAAGCGCTCGCCAAGGCCGAAGAAGAAGGCC 
&gt;Contig3
GCTGCGGCGCTGATCCTGGCGGCCCGCGCCGAG
</code></pre>

<p>I have used the following codes to extract the sequences from <code>genome.fasta</code> based on the coordinates in <code>id.txt</code> by using the following command <code>xargs samtools faidx genome.fasta &lt;id.txt &gt; result.fasta</code>
It has extracted only <code>Contig2:5-10</code>sequence. Because, the coordinates of <code>Contig2:5-10</code>is in proper order. But <code>Contig3:15-7</code> and <code>Contig1:12-3</code> coordinates are in inverse/reverse order, So samtools could not fetch those sequences.</p>

<p>I need to extract those sequences from <code>Contig1:12-3</code> and <code>Contig3:15-7</code> coordinates and also I want to <code>reverse complement</code> them. I have  larger fasta file and plenty of coordinates sequences to be extracted, So I need to automate this process. So, please help me to automate this process.</p>

<p>The expected output is</p>

<pre><code>&gt;Contig3    
GATCAGCGC
&gt;Contig2
AGCGCT
&gt;Contig1
CTTGATGGCC
</code></pre>

<p>Thank you in advance.</p>

CC BY-SA 4.0
2019-06-13T13:12:25.550Z
Welcome to the site. Do you have any preference about the tool to be used? What else have you tried to set the coordinates in the right order? Usually this kind of problems are solved by working on every tiny problem one after the other.
llrs
CC BY-SA 4.0
2019-06-13T13:15:17.153Z
I have used linux command the following command to keep the coordinates in right order,                                                                                                                                           "awk '$2 > $3 {printf( "%s\t%s\t%s\n", $1, $3, $2); next} id.csv".
Kumar
CC BY-SA 4.0
2019-06-13T13:17:20.050Z
So now you have the right order and the positions. What is stopping you from extracting the sequences?
llrs
CC BY-SA 4.0
2019-06-13T13:18:20.610Z
By using the above mentioned command I formatted the coordinates in right order and made it compatible for samtools. But still, one more thing I have to do is reverse complement of those right ordered sequences. I have jumbled coordinates, so the command changes the order wherever it needed.
Kumar
CC BY-SA 4.0
2019-06-13T13:22:47.757Z
I need to extract the sequence based on the coordinates specified above (even though it is mixed up with inverse ordered and right ordered). But, Now I need to reverse complement the right ordered (right ordered by using above command) sequences. So that I could go for my further analysis. Moreover, I need to automate this process. The reverse complement implementation to be taken place, wherever the inverse coordinates are there.
Kumar
CC BY-SA 4.0
2019-06-13T14:30:18.180Z
What do you know biopython? It was designed to make reverse-complement sequences and translate if needed
llrs
<python><fasta><samtools><shell>
How to reverse complement the DNA sequences for given inverse/reverse coordinates?
# How to reverse complement the DNA sequences for given inverse/reverse coordinates?

(Asked by: username_0 on 2019-06-13)

I have the series of coordinates in id.txt file, whose coordinates sequences are in genome.fasta file.
The coordinates of id.txt file is shown below,
Contig3:15-7
Contig2:5-10
Contig1:12-3

The genome.fasta file is shown below,
>Contig1
AAGGCCATCAAGGACGTGGATGAGGTCGTCAAG
>Contig2
ACGAAGCGCTCGCCAAGGCCGAAGAAGAAGGCC 
>Contig3
GCTGCGGCGCTGATCCTGGCGGCCCGCGCCGAG

I have used the following codes to extract the sequences from genome.fasta based on the coordinates in id.txt by using the following command xargs samtools faidx genome.fasta  result.fasta
It has extracted only Contig2:5-10sequence. Because, the coordinates of Contig2:5-10is in proper order. But Contig3:15-7 and Contig1:12-3 coordinates are in inverse/reverse order, So samtools could not fetch those sequences.
I need to extract those sequences from Contig1:12-3 and Contig3:15-7 coordinates and also I want to reverse complement them. I have  larger fasta file and plenty of coordinates sequences to be extracted, So I need to automate this process. So, please help me to automate this process.
The expected output is
>Contig3    
GATCAGCGC
>Contig2
AGCGCT
>Contig1
CTTGATGGCC

Thank you in advance.

**Tags:** python, fasta, samtools, shell
> **Comments:**
> **username_1**: Welcome to the site. Do you have any preference about the tool to be used? What else have you tried to set the coordinates in the right order? Usually this kind of problems are solved by working on every tiny problem one after the other.  
> **username_0**: I have used linux command the following command to keep the coordinates in right order,                                                                                                                                           "awk '$2 > $3 {printf( "%s\t%s\t%s\n", $1, $3, $2); next} id.csv".  
> **username_1**: So now you have the right order and the positions. What is stopping you from extracting the sequences?  
> **username_0**: By using the above mentioned command I formatted the coordinates in right order and made it compatible for samtools. But still, one more thing I have to do is reverse complement of those right ordered sequences. I have jumbled coordinates, so the command changes the order wherever it needed.  
> **username_0**: I need to extract the sequence based on the coordinates specified above (even though it is mixed up with inverse ordered and right ordered). But, Now I need to reverse complement the right ordered (right ordered by using above command) sequences. So that I could go for my further analysis. Moreover, I need to automate this process. The reverse complement implementation to be taken place, wherever the inverse coordinates are there.  
> **username_1**: What do you know biopython? It was designed to make reverse-complement sequences and translate if needed  

---

## 3 answers

### Answer by username_2 on 2019-06-13. Score: 1 | Accepted answer

The exonerate tool comes with some very handy utilities for sequence manipulation. Here, we are interested in fastasubseq and fasterevcomp. First, install exonerate as described here. If you're using Debian or Ubuntu or one of their derivatives, it might be enough to run
sudo apt install exonerate

Then, follow these steps:

Create one fasta file for each sequence in genome.fasta
fastaexplode genome.fasta 

A little preprocessing
$ awk -F[:-] '{ 
                if($2>$3){
                    start=$3-1; 
                    len=$2-start" -"
                }
                else{
                    start=$2-1; 
                    len=$3-start
                } 
                print $1,start,len}' id.txt 
Contig3 7 8 -
Contig2 5 5
Contig1 3 9 -

So, that will print the target contig's name, the start coordinate of the target sequence, the length of the target sequence and a - if it is on the reverse strand. The fastasubseq program I will use below starts counting from 0 and not from 1, so we also need to subtract 1 from the coordinates to get them right.
Now, feed that into fastafetch and fastarevcomp
tmpFile=$(mktemp); 
    > subseqs.fa
    awk -F[:-] '{ 
                if($2>$3){
                    start=$3-1; 
                    len=$2-start" -"
                }
                else{
                    start=$2-1; 
                    len=$3-start
                } 
                print $1,start,len}' id.txt |
    while read cont start len rev; do 
        fastasubseq "$cont".fa $start $len > $tmpFile; 
        if [[ -n $rev ]]; then 
            fastarevcomp "$tmpFile" >> subseqs.fa; 
        else 
            cat "$tmpFile" >> subseqs.fa; 
    fi
    done

This will produce the following file when run on your data:
$ cat subseqs.fa 
>Contig3:subseq(6,9) [revcomp]
GATCAGCGC
>Contig2:subseq(4,6)
AGCGCT
>Contig1:subseq(2,10) [revcomp]
CTTGATGGCC

Finally, remove the per-contig fasta files
rm Contig*fa

> **Comments:**
> **username_0**: Thank you so much terdon for your time and help. Your script works awesome. However, I did a mistake in my question. Instead of mentioning Contig3:15-7, Contig1:12-3, I have mentioned Contig3:15-7 and Contig2:5-10.                                                                                                                                                      I need to reverse complement the sequences, only if it in reverse order/inverse order. Which ever sequences are not in reverse order (Properly order sequences), only to be extracted as its not to be reverse complemented.  
> **username_0**: Please accept my sincere apologies for the same. Could you please do the needful.  
> **username_2**: @K.Dineshkumar I'm sorry, but I don't understand the problem. The script will only reverse complement if the start position is smaller than the end position. Isn't that what you want? In the example above, `Contig3:15-7` and `Contig1:12-3` have been reverse complimented, but `Contig2:5-10` has not.  
> **username_2**: @K.Dineshkumar please edit your question and add this. Show us the expected output. It's very hard to read in the comments.  
> **username_2**: @K.Dineshkumar   Never mind, I found the problem. The `fastasubseq` tool uses 0-based coordinates (the first position of a sequence is 0), but you have 1-based (the first position of a sequence is 1). Try the updated answer.  
> **username_0**: @terdon Please see my edited question, wherein I have included my expected outcome.  
> **username_0**: @terdon, The revised/updated script printing one base less in the end of the sequence. For example,  `Contig2:5-10` says, print the sequence from `5 to 10` But the script print only `5-9`.  
> **username_2**: @K.Dineshkumar argh! Yes, of course, sorry. Fixed now, try the (last, I hope) update. And thanks for adding the output to your question, now I can check :)  
> **username_0**: @terdon, Now it works perfect. I really appreciate your help. Thank you.  
> **username_2**: @K.Dineshkumar you're welcome. If this answers solves your issue, please take a moment to [accept it](//$$SITEURL/help/accepted-answer) by clicking on [the checkmark on the left](https://meta.stackexchange.com/a/5235/203101). That will mark the question as answered and is the way that thanks are conveyed on the Stack Exchange sites. You don't need to do it now, you might want to wait in case anyone can give a better solution, but please remember to accept an answer since that marks the question as answered.  
> **username_2**: @K.Dineshkumar thanks for accepting! Feel free to accept one of the other answers though, I won't take it the wrong way :) They give more efficient solutions.  

---

### Answer by username_3 on 2019-06-14. Score: 2

Since you tagged python, here is a solution that solves this in a single script using Biopython:
import sys
from Bio import SeqIO

d = {}
with open('id.txt') as f:
    for line in f:
        contig, coords = line.split(':')
        start, end = map(int, coords.split('-'))
        d[contig] = start, end

for record in SeqIO.parse('genome.fasta', 'fasta'):
    start, end = d[record.id]
    if start > end:
        record.seq = record.seq[end - 1: start].reverse_complement()
    else:
        record.seq = record.seq[start - 1: end]
    SeqIO.write(record, sys.stdout, 'fasta')

Output on your example:
>Contig1
CTTGATGGCC
>Contig2
AGCGCT
>Contig3
GATCAGCGC

> **Comments:**
> **username_0**: Thank you @Chris_rands  
> **username_0**: I try to use it for a bacterial genome. But it shows error as follows,            `Traceback (most recent call last):   File "/home/ga/Documents/commands/scripts/amphora/reverse_comple_seq_extrac.py", line 12, in <module>     start, end = d[record.id] KeyError: '1'`  
> **username_3**: @K.Dineshkumar try using `d[record.description]` instead, i can't know what will work without seeing your data though  
> **username_0**: My genome of interest can be obtained from the following link `https://www.ncbi.nlm.nih.gov/Traces/wgs/CAHO01?display=contigs`  
> **username_0**: My coordinate is as follows,                                                                                                  `CAHO01000012.1:92353-95295 CAHO01000017.1:27362-27949 CAHO01000014.1:11449-12012 CAHO01000018.1:51156-52928 CAHO01000030.1:58079-56889 CAHO01000001.1:60348-62039 CAHO01000026.1:45206-44508 `  
> **username_0**: When I tried with `d[record.description]` instead `of d[record.id]`, I end up the error as follows,                                                                                                                                      `Traceback (most recent call last):   File "/home/ga/Documents/commands/scripts/amphora/reverse_comple_seq_extrac.py", line 12, in <module>     start, end = d[record.description] KeyError: 'CAHO01000001.1 Xanthomonas citri pv. mangiferaeindicae LMG 941 WGS project CAHO00000000 data, contig XMINcontig00001, whole genome shotgun sequence'`  

---

### Answer by username_4 on 2019-06-14. Score: 1

If you're okay with using a Perl script, you could convert your headers into stranded BED6 and use the following script bed2faidxsta.pl with samtools-indexed FASTA to get back stranded sequence:
https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640
This script looks at the sixth column of your BED file of chromosome, start, and stop positions, to decide whether to do a reverse complement of the queried sequence. 
You could use awk to quickly build such a BED6 file:
$ awk -vFS='[:-]' '($0~/^>/){ if ($2 < $3) { strand = "+"; } else if ($2 > $3) { strand = "-"; } else { printf("Coordinates must differ\n"); exit 1; } printf("%s\t%d\t%d\t.\t0\t%c\n", $1, $2, $3, strand; }' in.fa > out.bed6

To set up indexed FASTA, here's a generic shell-based approach, which uses [Kent utilities][1] fetchChromSizes to get a list of chromosomes for the build mm9:
$ mkdir mm9
$ cd mm9
$ while read chr; do echo "$chr"; wget -qO- "http://hgdownload.cse.ucsc.edu/goldenpath/mm9/chromosomes/$chr.fa.gz" | gunzip -c - > $chr.fa; samtools faidx $chr.fa; done < <(fetchChromSizes mm9 | cut -f1)

This can be adjusted for mm10, hg19, hg38, and so on.
Once you have indexed FASTA and the above Perl script, you could do something as easy as the following:
$ ./bed2faidxsta.pl --fastaDir="/path/to/indexed/mm9" < in.bed > out.fa

If you're working with a reference genome, the nice thing about this approach is that awk, Perl, and indexed lookups are fast. You can reuse indexed FASTA for other projects that also require sequence lookups by genomic interval.

> **Comments:**
> **username_0**: Thank you @Alex Reynolds  

---

