How to reverse complement the DNA sequences for given inverse/reverse coordinates?

Question

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.

Answer 1

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:

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

   fastaexplode genome.fasta 
   

2. 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.

3. 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
   

4. Finally, remove the per-contig fasta files

   rm Contig*fa
   

Answer 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

Answer 3

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.