编写一个接受 fasta 并反转所有序列的 Perl 脚本（没有 BioPerl）？

Question

我不知道这是否只是 Stawberry Perl 的一个怪癖，但我似乎无法做到运行。我只需要做一个 fasta 并反转其中的每个序列。

-问题-

我有一个 multifasta 文件：

>seq1
ABCDEFG
>seq2
HIJKLMN

预期输出为：

>REVseq1
GFEDCBA
>REVseq2
NMLKJIH

脚本在这里：

$NUM_COL = 80; ## set the column width of output file
$infile = shift; ## grab input sequence file name from command line
$outfile = "test1.txt"; ## name output file, prepend with “REV”
open (my $IN, $infile);
open (my $OUT, '>', $outfile);
$/ = undef; ## allow entire input sequence file to be read into memory
my $text = <$IN>; ## read input sequence file into memory
print $text; ## output sequence file into new decoy sequence file
my @proteins = split (/>/, $text); ## put all input sequences into an array


for my $protein (@proteins) { ## evaluate each input sequence individually
    $protein =~ s/(^.*)\n//m; ## match and remove the first descriptive line of
    ## the FATA-formatted protein
    my $name = ; ## remember the name of the input sequence
    print $OUT ">REV$name\n"; ## prepend with #REV#; a # will help make the
    ## protein stand out in a list
    $protein =~ s/\n//gm; ## remove newline characters from sequence
    $protein = reverse($protein); ## reverse the sequence

    while (length ($protein) > $NUM_C0L) { ## loop to print sequence with set number of cols

    $protein =~ s/(.{$NUM_C0L})//;
    my $line = ;
    print $OUT "$line\n";
    }
    print $OUT "$protein\n"; ## print last portion of reversed protein
}

close ($IN);
close ($OUT);
print "done\n";

Answer 1

这会按照你的要求去做

它从 FASTA 文件中构建一个散列 %fasta，保持数组 @keys 以保持序列的顺序，然后打印出散列的每个元素

序列的每一行在添加到哈希之前使用reverse反转，使用unshift以相反的顺序添加序列的行

程序期望输入文件作为命令行参数，并将结果打印到STDOUT，可在命令行重定向

use strict;
use warnings 'all';

my (%fasta, @keys);

{
    my $key;

    while ( <> ) {

        chomp;

        if ( s/^>\K/REV/ ) {
            $key = $_;
            push @keys, $key;
        }
        elsif ( $key ) {
            unshift @{ $fasta{$key} }, scalar reverse;
        }
    }
}

for my $key ( @keys ) {
    print $key, "\n";
    print "$_\n" for @{ $fasta{$key} };
}

输出

>REVseq1
GFEDCBA
>REVseq2
NMLKJIH

更新

如果您更喜欢重新包装序列以便短行在末尾，那么您只需要重写转储哈希的代码

该方案以原文件中最长行的长度为限，将反转后的序列重新包装成相同的长度。更清楚的是，指定一个显式长度而不是计算它会很简单

您需要在程序顶部添加 use List::Util 'max'

my $len = max map length, map @$_, values %fasta;

for my $key ( @keys ) {
    print $key, "\n";
    my $seq = join '', @{ $fasta{$key} };
    print "$_\n" for $seq =~ /.{1,$len}/g;
}

给定原始数据，输出与上述解决方案相同。我用这个作为输入

>seq1
ABCDEFGHI
JKLMNOPQRST
UVWXYZ
>seq2
HIJKLMN
OPQRSTU
VWXY

有了这个结果。所有行都换行到 11 个字符 - 原始数据中最长 JKLMNOPQRST 行的长度

>REVseq1
ZYXWVUTSRQP
ONMLKJIHGFE
DCBA
>REVseq2
YXWVUTSRQPO
NMLKJIH

Answer 2

我不知道这是否仅适用于使用玩具数据集的 class 或大小可能达到千兆字节的实际研究 FASTA。如果是后者，那么不像您的程序和 Borodin 那样将整个数据集保存在内存中是有意义的，而是一次读取一个序列，将其打印出来并忘记它。以下代码执行此操作，还处理 may have asterisks as sequence-end markers 的 FASTA 文件，只要它们以 > 开头，而不是 ;.

#!/usr/bin/perl
use strict;
use warnings;

my $COL_WIDTH = 80;

my $sequence = '';
my $seq_label;

sub print_reverse {
    my $seq_label = shift;
    my $sequence = reverse shift;
    return unless $sequence;
    print "$seq_label\n";
    for(my $i=0; $i<length($sequence); $i += $COL_WIDTH) {
        print substr($sequence, $i, $COL_WIDTH), "\n";
    }
}

while(my $line = <>) {
    chomp $line;
    if($line =~ s/^>/>REV/) {
        print_reverse($seq_label, $sequence);
        $seq_label = $line;
        $sequence = '';
        next;
    }
    $line = substr($line, 0, -1) if substr($line, -1) eq '*';
    $sequence .= $line;
}
print_reverse($seq_label, $sequence);

编写一个接受 fasta 并反转所有序列的 Perl 脚本（没有 BioPerl）？

Write a Perl script that takes in a fasta and reverses all the sequences (without BioPerl)?

perl

bioinformatics

fasta

输出

更新