text: user4712: 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 <id.txt > 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.
  
  user298: 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
  
  
  
  user776: 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.
  
  user104: 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
  
  
original_text: 4712: <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>
   298: <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>
   776: <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>
   104: <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>
  