1
$\begingroup$

I'm looking to splice out specific amino acid alignment positions as a continuous motif from an alignment (~15 000 sequences). The code is written in Perl (I agree not many developers now) but splicing method is awk which a lot of you know. Please note the hash is "tied" (not shown), this just means the order of the sequences is preserved. I can get the awk code below to return single amino acid positions.

I am struggling to get a loop working in awk.

I agree the Perl wrapper makes it more complicated, but its part of a huge data pipeline (all written in Perl).

#!/usr/local/ActivePerl-5.20/bin/env perl
#use strict;
use warnings;

# I loaded the alignment into a hash, where key = sequence and value = protein
my %dna = (
    virus1 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF", 
    virus2 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF", 
    virus3 => "LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF"
);

# I define the start and end of the alignment positions I am splicing
my $n = 1;
my $z = 5;

# I turn the alignment into a string with separated by carriage returns 
my @array = map {"$dna{$_}"} keys %dna;
my $string = join "\n", @array;

# This is the awk loop, note in Perl the "``" means issue a shell command 

my $splice = `echo "$string" | awk 'BEGIN {FS=""} {x = $n; while (x < $z); print \$ x; ++x}' 2>1`;
print $splice;
__DATA__

Expected output

LGVLV
LGVLV
LGVLV

Observed output I get the original sequences returned

I am familiar with Python/Perl solutions of line by line processing but awk is more efficient IMO.

$\endgroup$
  • $\begingroup$ Could you edit and show the output you're expecting? I speak both Perl and awk, but I don't get what output you want. Also, I really doubt that calling awk will be more efficient than doing it natively in perl. You don't seem to be using anything that awk is particularly good at. $\endgroup$ – terdon Feb 25 at 0:44
  • $\begingroup$ Amended, thanks. I will post the solution I got working, but this is because I don't understand awk $\endgroup$ – Michael G. Feb 25 at 0:49
  • $\begingroup$ The main problem is that your while only encloses the first statement after it, and not the second. You need { } for that. If you could explain what that awk is supposed to be doing (I'm pretty sure you've copied it wrong, is there really a space between \$ and x?), I should be able to post an answer. Do you just want the first 5 characters of each sequence? Don't you want positions $n to $z (3-7)? $\endgroup$ – terdon Feb 25 at 0:54
1
$\begingroup$

I think you are trying to use awk to get characters 3 through 7 (inclusive) of each input sequence. If so, a working version of your approach would be this:

my $splice = `echo "$string" | awk 'BEGIN {FS=\"\"} {x = $n; while (x < $z){ printf \"%s\",\$x; ++x} print ""}' 2>1`;

That, when added to your perl script, produces:

$ foo.pl
VLVIL
VLVIL
VLVIL

That isn't the output you show, but the output you show is just the first 5 letters of each sequence, which doesn't seem to correspond to what your script is doing.

Now, the reason your version didn't work, is that the while was like so:

while (x < $z); print \$ x; ++x}

Assuming the \$ x was a typo, the problem is that the ++x isn't part of the while loop. So that should never even have exited. Seeing as we're in the same time zone, and given the time, I am guessing you're tired and that you miscopied something when trying to give us a reproducible example. In any case, to get the loop right, you need {}:

while (x < $z){print \$x; ++x}

However, if all you want is a substring, a loop is a very inefficient way of doing that. You can instead use awk's substr command:

  substr(s, m[, n ])
             Return  the at most n-character substring of s that begins at position
             m, numbering from 1. If n is omitted, or if n specifies  more  charac‐
             ters than are left in the string, the length of the substring shall be
             limited by the length of the string s.

So, in your script, that would be:

my $splice = `echo "$string" | awk '{print substr(\$0,$n,5)}'

Since you're already using a language that excels in string manipulation, I really don't understand why you would try to use awk instead. It will not be faster or more efficient than the same thing in perl:

print substr($_,$n,5) . "\n" for @array;

By calling awk you are also calling the shell which will be slower by definition, especially if you're echoing thousands of lines. In fact, I tried this with 15000 sequences (just repetitions of your first one) and using either the awk loop I gave above, or the awk+substr, or the native perl substr:

  1. awk loop, average time taken to run the same thing 100 times: 0.0525
  2. awk substr, average time taken to run the same thing 100 times: 0.0522s
  3. perl substr, average time taken to run the same thing 100 times: 0.0464

And you can make it even faster by skipping the array and iterating over the keys directly:

print substr($dna{$_},$n,5) . "\n" for keys(%dna);

That averaged at 0.0376s over 100 runs.

$\endgroup$
  • $\begingroup$ Huge thanks @terdon. You are right I should have performed a runtime check :-( such as use Timer::Runtime. I worked on the implicit assumption awk was faster for this operation. David Cross "lectured" on runtime in detail in "Data Munging for Perl": There is no denying that 0.0525 vs. 0.0376s is 30% faster. $\endgroup$ – Michael G. Feb 25 at 14:51
  • 1
    $\begingroup$ @MichaelG. in my experience, awk is blindingly fast when your data are organized in fields because of its implicit splitting. If not, however, it doesn't give any real advantage over a language like Perl which is simply brilliant for text parsing. $\endgroup$ – terdon Feb 25 at 15:29
1
$\begingroup$

If you have the sequences in a seperate file like this:

$ cat input
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVAHLALIAAFKVRPALLVSFIF
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGDVXXXXXXAAFKVRPALLVSFIF
LGVLVILLMVQEGLKKRMTTKIIISTSMAVLVAMILGGFSMSDLAKLAILMGATFAEMNTGGZZZZZZLIAAFKVRPALLVSFIF

you can use a simple cut to fetch a substring:

$ cut -c1-5 input
LGVLV
LGVLV
LGVLV
$\endgroup$

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.