1

我经常使用具有以下格式的生物序列数据 (FASTA)，其中前导左尖括号用作分隔符以指示新的序列标题。这些文件通常有文本换行（标题除外）：

>header1
ACTGACTGACTGACTG
ACTGACTGACTGACTG
>header2
CTGGGACTAGGGGGAG
CTGGGACTAGGGGGAG

通常，我想避免将整个文件读入内存，因为它可能有很多 MB（有时是 GB），所以我尝试专注于 while 循环并逐行读取。但是，我发现自己经常需要添加额外的代码来在文件的顶部或底部做一些独特的事情。例如，今天我想删除某个文件的文本换行，这看起来很简单：

while (my $line = <$inputfasta_fh>) {
    chomp($line);
    if ($line =~ /^>/) {
        print $outputfasta_fh "$line\n";
    }
    else {
        print $outputfasta_fh $line;
    }
}

但是，我意识到我需要在除第一个之外的所有标题之前添加一个换行符（否则它们将被连接到前一个序列的末尾）。所以，这是我粗略的解决方法。

my $switch = 0;
while (my $line = <$inputfasta_fh>) {
    chomp($line);
    if ($line =~ /^>/) {
        if ($switch == 1) {
            print $outputfasta_fh "\n";
        }
        print $outputfasta_fh "$line\n";
        $switch = 1;
    }
    else {
        print $outputfasta_fh $line;
    }
}

以前，我遇到过其他问题，需要对最后一行做一些事情。例如，我有一个脚本可以读取 fasta，存储每个标题，然后开始计算其序列长度（再次逐行），如果它在我指定的范围内，我将其保存到另一个文件中。如果长度超过最大值，计数将中止，但我不知道它是否超过最小值，直到我到达另一个标题或文件末尾。在后一种情况下，我不得不在 while 循环下面重复长度检查子例程。我想避免重复最后一部分。

my $length = 0;
my $header;
my @line_array;

while (my $line = <$inputfasta_fh>) {
    chomp($line);
    if ($line =~ /^>/) {
        # check if previous sequence had a length within range
        if (check_length($length, $minlength, $maxlength) == 1) {
            print $outputfasta_fh "$header\n";
            print $outputfasta_fh join ("\n", @line_array), "\n";
        }
        undef @line_array;
        $header = $line;
        $length = 0;
    }
    else {
        if ($length <= $maxlength) { # no point in measuring any more
            push (@linearray, $line);
            $length += length($line);
        }
    }
}

#and now for the last sequence
if (check_length($length, $minlength, $maxlength) == 1) {
    print $outputfasta_fh "$header\n";
    print $outputfasta_fh join ("\n", @line_array), "\n";
}

sub check_length {
    my ($length, $minlength, $maxlength) = @_;
    if (($length >= $minlength) && ($length <= $maxlength)) {
        return 1;
    }
    else {
        return 0;
    }
}

所以，我的基本问题是如何表明我想在循环中做某事而不诉诸计数器或在循环外重复代码？谢谢你的帮助！

4

2 回答 2

3

以下是您描述的 2 个问题的解决方案。它们是使用BioPerl发行版中的模块解决的。在这种情况下，Bio::SeqIO模块用于打开文件，Bio::Seq模块用于它提供的某些方法（长度、宽度）。您可以看到他们如何简化解决方案！

#!/usr/bin/perl
use strict;
use warnings;
use Bio::SeqIO;

my $in  = Bio::SeqIO->new( -file   => "input1.txt" ,
                           -format => 'fasta');
my $out = Bio::SeqIO->new( -file   => '>test.dat',
                           -format => 'fasta');

while ( my $seq = $in->next_seq() ) {
    $out->width($seq->length); # sequence on 1 line.
    $out->write_seq($seq);
}

my ($minlen, $maxlen) = (40, 1000);

while ( my $seq = $in->next_seq() ){
    my $len = $seq->length;
    out->write_seq($seq) if $minlen <= $len && $len <= $maxlen;
}

看看这些模块是值得的——正如您从这两个示例中看到的那样，生成的代码更加简洁且易于理解。你可以看看BioPerl wiki。HOWTO提供了一些您可以立即使用的示例。

于 2013-07-06T04:27:18.797 回答
1

目前尚不清楚您到底想要实现什么。
但是，如果您确定特殊情况是第一行或最后一行，您有几种方法可以处理它：

不需要常规处理的特殊一线

Process first line
$line = <$INPUT>;
... process line

Regular processing
while(<$INPUT>) {
... process lines
}

也需要定期处理的特殊一线

Process first line
$line = <$INPUT>;
... process line

Regular processing
do {
... process lines
} while(<$INPUT>);

特殊的最后一行，

在这里，您没有办法事先识别最后一行，因此您必须在循环中进行（除非您确切知道有多少行并for为第一个 N-1 使用循环，然后分别处理最后一行)

while(<$INPUT>) {
   break if islastline();
   ... process lines
}
... process last line

或者

while(<$INPUT>) {
   ... process lines
   break if islastline();
}
... process last line

或者

for($i=0; $i<N-1 ; $i++) {
   $line = <$INPUT>;
   ...process lines
}
$line = <$INPUT>
... process last line

您描述的另一种情况，您需要计数并且完成后，循环继续但您不再需要计数是不同的。如果您担心代码看起来“干净”的计数，只需将循环分成两部分：

内部临时处理

first part does the whole package
while(<$INPUT>) {
   ...regular processing
   ...special processing
   break if specialProcessingDone();
}

second part does not need to do special processing anymore
while(<$INPUT>) {
   ...regular processing
}
于 2013-07-06T02:45:12.193 回答