Beefy Boxes and Bandwidth Generously Provided by pair Networks
Problems? Is your data what you think it is?

Re^2: Reverse Complement

by stamp1982 (Novice)
on Jun 30, 2013 at 20:54 UTC ( #1041650=note: print w/replies, xml ) Need Help??

in reply to Re: Reverse Complement
in thread Reverse Complement

Do you have to use STDIN in this case or do I have to store it a file or just run the code as it is? How did you get your outputs? I run the code and it did not allow any user input

please two sequences as arguments on the command line. Press any key to continue

what am I missing and why cant I run the code?


Replies are listed 'Best First'.
Re^3: Reverse Complement
by AnomalousMonk (Chancellor) on Jun 30, 2013 at 21:34 UTC

    As posted, the code is a cut-and-paste of a command line session. It is intended to show source code and invocation switches, pragmata, etc., command line input (to @ARGV), and command line output. The source code is everything between the pair of  " (double-quotes). (Note that these are the only double-quotes in the posted code section.)

    I don't know what you are doing to run the code, but try this:

    1. Download the posted code section using the  [download] link just after it.
    2. Edit out the first " and everything before it, and also the last " and everything after it. What remains should be the entire source code.
    3. Save the edited file.
    4. Invoke the edited file with the samealmost the same command line invocation used originally (Update: do not include the  -e switch as I originally wrote below!) and with the sample input you want, e.g.:
          perl -wMstrict -l  GGGGaaaaaaCatttatatat  atatataaatttttATtcccc
    5. Examine output. Experiment with various different input data. Play with variations on the source code. Enjoy!

    In fact, here's the edited section of code that should be all ready to go, just use the  [download] link and save to a file of your choice:

    die 'Please enter exactly two sequences' unless @ARGV == 2; ;; my ($seq1, $seq2) = map uc, @ARGV; die 'Please enter only ATCG sequences' if grep /[^ATCG]/, $seq1, $seq2 +; ;; die 'Sequence lengths differ' if length $seq1 != length $seq2; ;; (my $comp = $seq2) =~ tr/ATCG/TAGC/; my $mask = $seq1 ^ reverse $comp; ;; if ($mask !~ m{ [^\000] }xms) { print 'The sequences are reverse-complementary'; exit; } ;; $mask =~ tr{\000-\377}{_^}; print 'Sequences are non-reverse-complementary at these locations:'; print qq{@ARGV}; print qq{$mask }, ''.reverse $mask;

    Update: I just downloaded and saved the code block immediately above in this reply and went through the steps listed above (except it doesn't have to be edited), and it works. Be sure not to use the  -e switch when invoking a file; it is used for command-line source code statements. HTH.

      So how can I get this code to read from a .txt file? Because it is asking me to enter a sequence at line 1
        So how can I get this code to read from a .txt file?

        You have to write code to open a file, read data from it (perhaps line-by-line with readline), and then operate on the data using techniques shown you by rjt and others and developed by yourself, and then write the transformed data as appropriate.

        I'm sorry to be so vague, but this is a very general question! Perhaps something like chromatic's freely-downloadable Modern Perl would help.

        So how can I get this code to read from a .txt file?

        This thread seems to be turning into a "How do I process a data file?" kind of discussion. If so, it would probably be better to start another thread (perhaps with a link back to this thread for reference) in SoPW with a request of that nature. Be sure to include info on the structure of the data in the file you wish to process.

      Can you please explain this in your comment?

      and with the sample input you want, e.g.: perl -wMstrict -l GGGGaaaaaaCatttatatat atatataaatttttATtcccc

        The command line invocation example
            perl -wMstrict -l  GGGGaaaaaaCatttatatat  atatataaatttttATtcccc
        was intended to show how to invoke a source code file (prepared as described in previous replies) and followed by two sequence strings  GGGGaaaaaaCatttatatat and  atatataaatttttATtcccc to be compared for reverse-complementarity. (BTW: I've tried this exact command line and it works.)

      My apologies for being so slow. At the moment this is what I have.
      perl -wMstrict -l GGGGaaaaaaCatttatatat atatataaattttt +ATtcccc die "Please enter two sequences as arguments on the command line\n" unless @ARGV == 2; my ($orig, $comp) = map { uc } @ARGV; die "Please enter only ATCG sequences\n" if grep { /[^ATCG]/ } $orig, +$comp; die "Sorry, the two sequences you have just entered are of different l +engths." . "\nPlease try again on the command line.\n" if length $orig != length $comp; # $comp will transform back to $orig iff it is the reverse complement $comp =~ y/ATCG/TAGC/; if ($orig eq reverse $comp) { print "Yes, the two sequences are reverse-complement of each other +.\n"; exit; } else { die "Unfortunately, the two sequences are not reverse-complement.\ +n"; }

      comments I get are : unquoted strings may clash with future reserved words and compilation error and syntax error at wMstrict -l

        ... this is what I have. ...

        I'm guessing, but the only thing that occurs to me is that you actually have the line
            perl -wMstrict -l  GGGGaaaaaaCatttatatat  atatataaatttttATtcccc
        as the first non-blank line in your source code file. If so, this line should not be in the source code file. This line is intended as the command line invocation of the source code file. Please see previous examples.

        The source code file (which I have named in these examples) should have everything from the first
            die "Please enter two sequences as arguments on the command line\n"
                 unless @ARGV == 2;
        statement to the very last  } closing curly brace on the final if-statement.

        If this is not the problem, please include exact copies of all warnings or error messages in future posts.

Log In?

What's my password?
Create A New User
Node Status?
node history
Node Type: note [id://1041650]
and all is quiet...

How do I use this? | Other CB clients
Other Users?
Others musing on the Monastery: (3)
As of 2017-12-15 03:26 GMT
Find Nodes?
    Voting Booth?
    What programming language do you hate the most?

    Results (416 votes). Check out past polls.