2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – llrs Jun 13 at 13:12
  • $\begingroup$ 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". $\endgroup$ – K. Dineshkumar Jun 13 at 13:15
  • $\begingroup$ So now you have the right order and the positions. What is stopping you from extracting the sequences? $\endgroup$ – llrs Jun 13 at 13:17
  • $\begingroup$ 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. $\endgroup$ – K. Dineshkumar Jun 13 at 13:18
  • $\begingroup$ 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. $\endgroup$ – K. Dineshkumar Jun 13 at 13:22
1
$\begingroup$

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
    
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – K. Dineshkumar Jun 14 at 6:29
  • $\begingroup$ Please accept my sincere apologies for the same. Could you please do the needful. $\endgroup$ – K. Dineshkumar Jun 14 at 6:31
  • $\begingroup$ @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. $\endgroup$ – terdon Jun 14 at 8:21
  • $\begingroup$ @K.Dineshkumar please edit your question and add this. Show us the expected output. It's very hard to read in the comments. $\endgroup$ – terdon Jun 14 at 8:39
  • $\begingroup$ @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. $\endgroup$ – terdon Jun 14 at 8:49
1
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Thank you @Chris_rands $\endgroup$ – K. Dineshkumar Jun 15 at 18:03
  • $\begingroup$ 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' $\endgroup$ – K. Dineshkumar Jul 6 at 13:55
  • $\begingroup$ @K.Dineshkumar try using d[record.description] instead, i can't know what will work without seeing your data though $\endgroup$ – Chris_Rands Jul 8 at 7:37
  • $\begingroup$ My genome of interest can be obtained from the following link https://www.ncbi.nlm.nih.gov/Traces/wgs/CAHO01?display=contigs $\endgroup$ – K. Dineshkumar Jul 8 at 8:40
  • $\begingroup$ 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 $\endgroup$ – K. Dineshkumar Jul 8 at 8:40
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you @Alex Reynolds $\endgroup$ – K. Dineshkumar Jun 15 at 18:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.