首页 文章

获得以FASTA中特定氨基酸开始的蛋白质序列的 Headers 行

提问于
浏览 次
0

嗨，大家好，所以我一直试图使用PERL只打印从FASTA文件以“MAD”或“MAN”（前3个aa）开头的蛋白质序列的 Headers （整个> 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的 Headers 是我正在寻找的

>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
perl fasta

2 回答

  • 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，位置3中的第4行，等等 . 因为've already incremented the index, you'重新打印 Headers 行之后的行 . 通过在循环结束时递增 $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++;
}

    我希望这有用而且清晰！

    回复于
  • 0

    你的打印声明有问题 . 应该是：

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

    但是，通常最好只逐行处理文件，除非有特定的原因需要啜饮整个文件：

    #!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]
    回复于

相关问题