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

Posted 2023-03-24

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

【英文标题】：Get the header lines of protein sequences that start with specific amino acid in FASTA

【发布时间】：2014-11-28 15:47:15

【问题描述】：

大家好，我一直在尝试使用 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

【问题讨论】：

print "@lines [$index-1] ... 没有多大意义。打印整个数组？ [$index-1] 应该是尝试打印上一行，或者实际上是在括号中打印 $index-1，例如如果您在第 10 行，则实际打印 [9]?

你在第 10 行，如果第 10 行满足要求，想打印第 9 行

那么你会想要$prev = $index - 1; print "$lines[$prev]"。

【参考方案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] 

【讨论】：

奇怪我把@lines[$index-1]中的@改成$，还是没有给出标题行

也许实际上共享了一些示例数据？您所说的标题行实际上是在匹配的数据之前吗？

###这是一小部分，第一个seq的标题是我要找的>gi|16128078|ref|NP_414627.1| UDP-N-乙酰胞壁酰-L-丙氨酰-D-谷氨酸：内消旋二氨基庚二酸连接酶[大肠杆菌str。 K-12 substr。 MG1655] MADRNLRDLLAPWVPDAPSRALREMTLDSRVAAAGDLFVAVVGHQADGRRYIPQAIAQGVAAIIAEAKDE ATDGEIREMHGVPVIYLSQLNERLSALAGRFYHEPSDNLRLVGVTGTNGKTTTTQLLAQWSQLLGEISAV MGTVGNGLLGKVIPTENTTGSAVDVQHELAGLVDQGATFCAMEVSSHGLVQHRVAALKFAASVFTNLSRD HLDYHGDMEHYEAAKWLLYSEHHCGQAIINADDEVGRRWLAKLPDAVAVSMEDHINPNCHGRWLKATEVN 跨度>

编辑您的问题以添加该信息。我认为如果下一行以 MAD 或 MAN 开头，您想要打印 >gi 行吗？

感谢您的帮助，但我只是想知道我的代码出了什么问题，希望从错误中吸取教训。

【参考方案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++;

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

【讨论】：

别忘了检查open是否失败。