代码之家 › 专栏 › 技术社区 › Tom

编写一个接受fasta并反转所有序列的Perl脚本(没有BioPerl)?

fasta bioinformatics perl

0

Tom · 技术社区 · 8 年前

-这个问题-

>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 = $1; ## 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 = $1;
    print $OUT "$line\n";
    }
    print $OUT "$protein\n"; ## print last portion of reversed protein
}

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

2 回复 | 直到 8 年前

1

2

Borodin 8 年前

这将按你的要求进行

%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

2

1

mbethke 8 年前

我不知道这是否只是一个使用玩具数据集或实际研究FASTA的类,其大小可以是千兆字节。如果是后者,最好不要像您的程序和Borodin的程序那样将整个数据集保存在内存中,而是一次读取一个序列,将其反向打印出来,然后忘记它。下面的代码实现了这一点,并且还处理了FASTA文件 may have asterisks as sequence-end markers 只要他们开始 > 不 ; .

#!/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);