awk for motif splicing from an amino acid alignment

Question

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.

Answer 1

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.

Answer 2

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