This is a test version of Biostars. For the public version, visit https://www.biostars.org.
fasta - reverse complement sequence

Hi.

I have a fasta file with multiple headers and want to get reverse complement sequences. I usually use FASTX-TOOLKIT, but I want to learn how to do with linux commands.

I tried it using 'awk' and 'if' something like below, but I got a tedious result.

awk '{if($1 ~ />/) print $0; else gsub("[ATCG]","[TAGC]");print $0}' TEST.fa

As you can see, the file has the same string in the header.

How can I get reverse complement sequences without altering the headers?

> ATCG_1
agtcgacATCGATtataggta
> ATCG_2
cgactgcagtcgacATCGACT

Thank you!

fasta reverse complement bash
• 18,485 views • link
updated by dugong • dugong • written by mbk0asis mbk0asis

This is the kind of thing one shouldn't bother doing with awk or similar tools. Sure, you can come up with a solution, but why bother?

• 1 views • link
written by Devon Ryan Devon Ryan

I have this script for quick rev comp, but it assumes pure sequence as input (no headers).

#!/bin/bash
if [ ! -z "$1" ]; then
    echo "$1" | tr "[ATGCatgc]" "[TACGtacg]" | rev
else
    echo ""
    echo "usage: rev_comp_seq DNASEQUENCE"
    echo ""
fi
• 1 views • link
written by 5heikki 5heikki

5 answers

using 'tr', and 'rev', 2 lines per sequence 

 cat input.fa | while read L; do  echo $L; read L; echo "$L" | rev | tr "ATGC" "TACG" ; done
• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thank you, Pierre! It works well. What exactly is the difference between $L and "$L"? How does it work?

• 1 views • link
written by mbk0asis mbk0asis

$L and "$L" there is no difference while read L: the programs reads one line (the header), we print the line; we read the very next line into L; we print the last L, we reverse and translate it.

• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thank you! I understood!

• 1 views • link
written by mbk0asis mbk0asis

Thank! But I also have problem: that command can't work for sequencs with several lines.

• 0 views • link
written by yangjie4546 • yangjie4546

Why does the second 'read L' read the next line? I don't see where the input is coming from.

• 0 views • link
written by 9882788787 • 9882788787

This flips ACGT in the phred scores!

• 0 views • link
written by ryandward • ryandward

Assuming 2 lines per sequence (header and sequence), if not linearize and try something like following or modify more to get how you need

cat fasta.fa | paste - - | awk -F'\t' -vOFS='\t' '{gsub("A", "T", $2); gsub("T", "A", $2); gsub("G", "C", $2); gsub("C", "G", $2); print}' | tr '\t' '\n'
• 1 views • link
written by venu venu

Thank you for the answer! However, the code didn't work because gsubs serially change the sequences. Firstly. all A's are turned to T's, then second 'gsub' change them back to A's. So finally, all the sequences become a pool of A's and G's. Can I use gsub to change multiple strings at once?

• 1 views • link
written by mbk0asis mbk0asis

Perl?

perl -nle'BEGIN {
    @map{ A, a, C, c, G, g, T, t } = ( T, t, G, g, C, c, A, a )
    }
    print /^>/ ?
        $_ :
              join //, map $map{ $_ }, split //, scalar reverse
    ' file.fa
• 1 views • link
written by iraun iraun

Try OpenGene (https://github.com/OpenGene/OpenGene.jl) to get reverse complement very easily with an operator ~

julia> using OpenGene

julia> seq = dna("AAATTTCCCGGGATCGATCGATCG")
dna:AAATTTCCCGGGATCGATCGATCG

julia> ~seq
dna:CGATCGATCGATCCCGGGAAATTT
• 1 views • link
written by chen chen

Based on Pierre's answer here is a working solution for multi-line fasta.

Changes:

  • You do not lose undetermined sequences "N"-s
  • Works on multi line fasta: it escapes transformation on lines starting with '>'

As a 'one'-liner:

cat play.fa | while read L; do if [[ $L =~ ^'>' ]]; then echo $L; else echo $L | rev | tr "ATGC" "TACG" ; fi ; done

As a bash function

function _revcomplement.file { cat $1 | while read L; do if [[ $L =~ ^'>' ]]; then echo $L; else echo $L | rev | tr "ATGC" "TACG" ; fi ; done } ;

Call:

revcomplement.file play.fa

Tried on play.fa:

>MACHU
AGTCACCTTTACCCGGTTTCANNN
AGTCGCCTTTACCCGGTTTCA
CCCCCGGGGGGGGGGGGGGGG
AGTCACCTTTACCCGGTTTCA
AGTCACCTTTACCCGGTTTCA
AGTCGCCTTTACCCGGTTTCA
>PICCHU
ACTGCAGACACAACTACGGGGTTGTGGAGAGCTTCACAGTGCAGCGGCGAGGTGAGCGCGGCGCGGGGCGGGGCCTGAGTCCCTGTGAGCTGGGAATCTGAGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGAGAGAGAGAGAGAGAGAGAGAGAGAGAGACAGAGAGAGAGAGAGCGCCATGTGTG
ACTGCAGACACAACTACGGGGTTGTGGAGAGCTTCACAGTGCAGCGGCGAGGTGAGCGCGGCGCGGGGCGGGGCCTGAGTCCCTGTGAGCTGGGAATCTGAGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGAGAGAGAGAGAGAGAGAGACAGAGAGACAGAGAGAGAGAGAGCGCCATCTGTGAGCATTTAGAATCCTCTCTATCCTGAGCAAGGA
AGTCGCCTTTACCCGGTTTCA

The .fa file has to end with a newline, otherwise the last line is not processed!

• 0 views • link
written by dugong • dugong

Log in to answer this question.

More posts like this