获取 FASTA 中以特定氨基酸开头的蛋白质序列的标题行

Question

大家好，我一直在尝试使用 PERL 仅打印 FASTA 文件中以“MAD”或“MAN”（前 3 个 aa）开头的蛋白质序列的标题（整个 >gi 行）。但我无法弄清楚哪个部分出了问题。 提前致谢！

#!usr/bin/perl
use strict;

my $in_file = $ARGV[0];
open( my $FH_IN, "<", $in_file );    ###open to fileholder
my @lines = <$FH_IN>;
chomp @lines;
my $index = 0;

foreach my $line (@lines) {
    $index++;
    if ( substr( $line, 0, 3 ) eq "MAD" or substr( $line, 0, 3 ) eq "MAN" ) {
        print "@lines [$index-1]\n\n";
    } else {
        next;
    }
}

这是FASTA文件的一小部分，第一个seq的头是我要找的

>gi|16128078|ref|NP_414627.1| UDP-N-acetylmuramoyl-L-alanyl-D-glutamate:meso-diaminopimelate ligase [Escherichia coli str. K-12 substr. MG1655] MADRNLRDLLAPWVPDAPSRALREMTLDSRVAAAGDLFVAVVGHQADGRRYIPQAIAQGVAAIIAEAKDE ATDGEIREMHGVPVIYLSQLNERLSALAGRFYHEPSDNLRLVGVTGTNGKTTTTQLLAQWSQLLGEISAV MGTVGNGLLGKVIPTENTTGSAVDVQHELAGLVDQGATFCAMEVSSHGLVQHRVAALKFAASVFTNLSRD HLDYHGDMEHYEAAKWLLYSEHHCGQAIINADDEVGRRWLAKLPDAVAVSMEDHINPNCHGRWLKATEVN

Answer 1

您的打印语句有问题。应该是：

print "$lines[$index-1]\n\n";

但是，除非有特定原因需要对整个文件进行 slurp，否则最好逐行处理文件：

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

my $file = shift;

#open my $fh, "<", $in_file;
my $fh = \*DATA;

while (<$fh>) {
    print if /^>/ && <$fh> =~ /^MA[DN]/;
}

__DATA__
>gi|16128078|ref|NP_414627.1| UDP-N-acetylmuramoyl-L-alanyl-D-glutamate:meso-diaminopimelate ligase [Escherichia coli str. K-12 substr. MG1655] 
MADRNLRDLLAPWVPDAPSRALREMTLDSRVAAAGDLFVAVVGHQADGRRYIPQAIAQGVAAIIAEAKDE
ATDGEIREMHGVPVIYLSQLNERLSALAGRFYHEPSDNLRLVGVTGTNGKTTTTQLLAQWSQLLGEISAV
MGTVGNGLLGKVIPTENTTGSAVDVQHELAGLVDQGATFCAMEVSSHGLVQHRVAALKFAASVFTNLSRD
HLDYHGDMEHYEAAKWLLYSEHHCGQAIINADDEVGRRWLAKLPDAVAVSMEDHINPNCHGRWLKATEVN
–

输出：

>gi|16128078|ref|NP_414627.1| UDP-N-acetylmuramoyl-L-alanyl-D-glutamate:meso-diaminopimelate ligase [Escherichia coli str. K-12 substr. MG1655] 

Answer 2

由于您想知道如何改进您的代码，这里是您的程序的注释版本，其中包含一些关于如何更改它的建议。

#!/usr/bin/perl
use strict;

您还应该添加use warnings pragma，它会启用警告（如您所料）。

my $in_file = $ARGV[0];

最好检查$ARGV[0] 是否已定义，如果未定义，则给出适当的错误消息，例如

my $in_file = $ARGV[0] or die "Please supply the name of the FASTA file to process";

如果没有定义$ARGV[0]，Perl 将执行die 语句。

open( my $FH_IN, "<", $in_file );  # open to fileholder

您应该检查脚本是否能够打开输入文件；您可以通过添加die 语句来使用与上一条语句类似的结构：

open( my $FH_IN, "<", $in_file ) or die "Could not open $in_file: $!";

特殊变量$! 保存有关文件无法打开的错误消息（例如，文件不存在、没有读取权限等）。

my @lines = <$FH_IN>;
chomp @lines;
my $index = 0;

foreach my $line (@lines) {
    $index++;
    if ( substr( $line, 0, 3 ) eq "MAD" or substr( $line, 0, 3 ) eq "MAN" ) {
         print "@lines [$index-1]\n\n";

这是脚本中的问题点。首先，访问数组中的项目的正确方法是使用$lines[$index-1]。其次，数组中的第一项位于索引 0，因此文件的第 1 行将位于@lines 中的第 0 位，第 4 行位于第 3 位等。因为您已经增加了索引，所以您正在打印 标题行之后的行。通过在循环结束时增加 $index 可以轻松解决此问题。

    }
    else {
       next;
    }

这里没有必要使用next，因为else 语句后面没有代码，所以告诉Perl 跳过循环的其余部分没有任何好处。

固定的代码如下所示：

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

my $in_file = $ARGV[0] or die "Please supply the name of the FASTA file to be processed";
open( my $FH_IN, "<", $in_file ) or die "Could not open $in_file: $!";
my @lines = <$FH_IN>;
chomp @lines;

my $index = 0;
foreach my $line (@lines) {
    if ( substr( $line, 0, 3 ) eq "MAD" or substr( $line, 0, 3 ) eq "MAN" ) {
        print "$lines[$index-1]\n\n";
    }
    $index++;
}

我希望这是有帮助和明确的！