Re^3: awk 2 perl oneliner

The only prob I have now is, when applied to my real data a ">" sign is missing at the beginning ... and one is too much at the end.

The reason for this is the way the records are being split. The first record (and only the first) will have a leading ">", and all records except for the last will have a trailing ">". Now, if you happen to not print out the first record (because it doesn't match) the initial ">" will be missing from the output file. Similarly, if you happen to not print the last record, there will be a stray ">" from the previous record.

The proper way to handle those corner cases would be to always remove any leading or trailing ">"s, and always add a new one for output:

...
# scan second file
while (<>) {
    s/^>//;  # remove leading  '>'
    s/>$//;  # remove trailing '>'

    print ">$_" if /^([\w\d]+)/ && $names{$1};
    #  add ^
}
[download]

(This way, you also no longer need to allow a leading ">" in the regex (the >? part). )

Comment on Re^3: awk 2 perl oneliner Select or Download Code

Replies are listed 'Best First'.
Re^4: awk 2 perl oneliner by space_agent (Acolyte) on Mar 07, 2010 at 12:57 UTC
Thank you very many, almut ;-). Now I have an efficient way to extract sequences from a fasta file (doing bioinformatics stuff). It is more than 10times faster than the awk oneliner from my first post. Here is almuts code slighly adjusted. It can be used to extract sequences from a multifasta file. `usage extract-seq.pl file1 file2 cat file1 AV12349 J234054 . . cat file2 >AV12349 AAACCGTACCTAGAACAGTA ACAGTACA >J434594 ACCGTTAGTACGATACA . .` [download] file1 contains ID to extract, file2 contains sequences in multifasta format. output will be written to the file extracted.seqs `#!/usr/bin/perl $output = "extracted.seqs"; open (OUT,">$output") or die; # create names lookup table from first file while (<>) { chomp; $names{$_} = 1; last if eof; } $/ = "\n>"; # set input record separator # scan second file while (<>) { s/^>//; #remove leading ">" s/>$//; #remove trailing ">" if (/^([\w\d]+)/ && $names{$1}) { print OUT ">$_";} }` [download]	[reply] [d/l] [select]

Replies are listed 'Best First'.

Re^4: awk 2 perl oneliner
by space_agent (Acolyte) on Mar 07, 2010 at 12:57 UTC

Thank you very many, almut ;-). Now I have an efficient way to extract sequences from a fasta file (doing bioinformatics stuff). It is more than 10times faster than the awk oneliner from my first post.
Here is almuts code slighly adjusted. It can be used to extract sequences from a multifasta file.

usage extract-seq.pl file1 file2
cat file1
AV12349
J234054
   .
   .
cat file2
>AV12349
AAACCGTACCTAGAACAGTA
ACAGTACA
>J434594
ACCGTTAGTACGATACA
    .
    .
[download]

file1 contains ID to extract, file2 contains sequences in multifasta format. output will be written
to the file extracted.seqs

#!/usr/bin/perl

$output = "extracted.seqs";

open (OUT,">$output") or die;

# create names lookup table from first file
while (<>) {
    chomp;
    $names{$_} = 1;
    last if eof;
}
$/ = "\n>";  # set input record separator

# scan second file
while (<>) {
     s/^>//;  #remove leading ">"
     s/>$//;  #remove trailing ">"

     if (/^([\w\d]+)/ && $names{$1}) { print OUT ">$_";}

}
[download]

[reply]
[d/l]
[select]