Perl code for reverse complement of DNA sequence

Question

I was trying to make a perl code to get the reverse complement of fasta sequences of DNA in a .fna file format. the sequence02C.fna file contains 100s of DNA sequences :

>adbca3e
TGCTCCCCACGCTTGCCTCTCCAGTACTCAACCAAAGCAGTCTCTAGAAAAACAGTTTCCAACGCAATACGATGGAATTCCACTTCCCAAATATCTC
>4c2a958
TCCCCACGCTTTCGCGCTTCAGCGTCAGTATCTGTCCAGTGAGCTGACTTCTCCATCGGCATTCCTACACAGTACTCTAGAAAAACAGTTTCTGCTC
>0639b5b
TCGCGCCTCAGTGTCCAACGCAATACGAGTTGCAGACCAGGACACATGGAATTCCACTTCCCTCTCCAGTACTCAACCAAAGCAGTCTCTAGAAAAG

I have used the following command which can open the file and make reverse but does not show the sequence ID (eg: >adbca3e) in the output.

The code is:

#!/usr/local/perl

open (NS, "sequence02C.fna");
while (<NS>) {
    if ($_ =~ tr/ATGC/TACG/) {print $_;}
}

output file is only the complementary of the sequence but not reverse. Additionally, it does not contain the sequence IDs ">adbca3e"

Could anyone please suggest the appropriate code to do this reverse complementary of this sequence at once and getting the result into an output file?

Maybe you should explain the difference between "complementary" and "reverse". — TLP, Apr 25 '21 at 17:46

ikegami · Answer 1 · 2021-04-26T00:27:13.210

You only print the lines that contained a A, T, G or C. You want to print every line, so the print shouldn't be conditional.

#!/usr/local/perl

use strict;               # Always
use warnings;             # Always

while (<>) {
    if (/^>/) {           # Only modify lines starting with ">".
       tr/ATGC/TACG/;
       $_ = reverse($_);  # You didn't reverse.
    }

    print;                # Print undconditionally.
}

(tr/// and print use $_ by default.)

Note I didn't open the file. You can use the program as follows:

perl program.pl sequence02C.fna >sequence02C_revcomp.fna

or

perl -i~ program.pl sequence02C.fna

The latter modifies the file in place. (Be careful! Test it first. It does make a backup, though.)

TLP · Answer 2 · 2021-04-25T22:14:11.663

1

You say you have a program to "make reverse", but it only gives complementary. Maybe that is a very obvious description for you, but it is not very clear to me.

If by "reverse" you mean to print the string backwards, just use the reverse function. Complementary I assume is taking the corresponding nucleobases, which is what your transliteration is meant to do tr/ATGC/TACG/.

To fix not printing ids, just remove the if condition on the print-statement.

What I would do is just use the diamond operator for a small program like this:

use strict;
use warnings;
use feature 'say';

while (<>) {
    chomp;
    unless (/^>/) {
        tr/ATGC/TACG/;            # transliterate non-ids
        my $reverse = reverse;    # reverse $_
        say $reverse;             # do something with $reverse
    }
    say;          # print current line
}

Then you can use this program like this:

$ perl program.pl sequence02C.fna > output.txt

edited Apr 25 '21 at 22:14

answered Apr 25 '21 at 18:03

TLP

66,756
10
92
149

Thanks @TLPfor the reply with codes. I have used the following program, just with the open command adding to the suggested program : open (NS, "sequence02C.fna"); while () { chomp; but found error output: [siddiquee@intron2 realtask]$ perl rc.pl syntax error at rc.pl line 13, near "my " Global symbol "$reverse" requires explicit package name at rc.pl line 13. Global symbol "$reverse" requires explicit package name at rc.pl line 14. Execution of rc.pl aborted due to compilation errors. – mashuk Apr 25 '21 at 20:12
Global symbol .. is when you do not declare a variable using `my`. In your case, sounds like you missed ending a line with semi-colon, so that `my` got corrupted. The error is not visible in the code you show in the comment however. – TLP Apr 25 '21 at 21:09
@mashuk There was actually a semi-colon missing in my code, added it now. – TLP Apr 25 '21 at 22:14

Perl code for reverse complement of DNA sequence

2 Answers2