使用Perl解析蛋白质结构域

2001年11月16日 by

Perl编程语言因其实用性而受到生物学家们的喜爱。在我的书《《生物信息学入门Perl》》中，我展示了生物学家想要编写的许多程序都可以轻松甚至愉快地用Perl完成。

我的书教生物学家如何用Perl编程，即使他们以前从未编程过。这篇文章将使用我在书中中后期的章节中找到的Perl水平，在了解了一些基础知识之后。然而，即使不知道任何编程的生物学家也可以阅读这篇文章。他们应该能够浏览这篇文章中的程序代码，只阅读注释，以了解Perl在实际应用中使用的情况，使用真实生物数据。

计算机上的生物数据通常是结构化的ASCII平面文件——也就是说，是纯文本文件——或者是在关系数据库中。这两种数据源都可以用Perl程序轻松处理。对于这篇文章，我将讨论一个平面文件数据源，即Prosite数据库，它包含有关蛋白质结构域的宝贵生物信息。我将演示如何使用Perl提取和使用蛋白质结构域信息。在《生物信息学入门Perl》中，我还展示了如何处理其他几个类似的数据源，包括GenBank（基因数据库）、PDB（蛋白质数据银行）、BLAST（基本局部比对搜索工具）输出文件和REBASE（限制性内切酶数据库）。

什么是Prosite？

Prosite代表“蛋白质位点和模式词典”。要了解更多关于Prosite背后迷人生物学的信息，请访问Prosite用户手册。以下是用户手册中对Prosite的简介描述：

“Prosite是一种确定从基因组或cDNA序列翻译的未知蛋白质功能的方法。它包含一个生物意义显著的位点和模式数据库，通过适当的计算工具可以快速且可靠地识别新序列属于哪个已知的蛋白质家族（如果有的话）。”

在某些情况下，未知蛋白质的序列与已知结构的任何蛋白质相比关系太远，以至于不能通过整体序列比对检测到其相似性。然而，可以通过其序列中特定残基类型的出现来识别，这些特定残基类型被称为模式、基序、签名或指纹。这些基序的出现是由于蛋白质特定区域的特定结构要求，这些区域可能对它们的结合性质或酶活性很重要。

Prosite作为一组纯文本文件提供数据以及文档。Prosite的主页提供了一个用户界面，允许您查询数据库并查看文档。您还可以从Prosite ftp站点获取数据库以供本地安装。对于非商业用户，其使用是免费的。这里涉及一些引人入胜且重要的生物学；接下来程序中的有趣和有用的Perl编程技巧。请参阅Prosite用户手册了解生物学背景，以及《生物信息学入门Perl》了解编程背景。或者继续阅读，以了解当将编程技能与生物数据相结合时可以做到什么程度。

Prosite数据

Prosite 数据可以下载到您的电脑上。它位于名为 prosite.dat 的 ASCII 平文件中，大小超过 4MB。为此篇文章创建的小版本文件名为 prosmall.dat，可以从此处获取。这个版本的数据仅包含完整文件中的前几条记录，这使得您更容易下载和测试，并且这是我们将在本文后面讨论的代码中使用的文件。

Prosite 还提供了一份配套的数据文件，prosite.doc，其中包含了 prosite.dat 中所有记录的文档。尽管我们不会在本文中使用它，但我还是建议您看看它，并考虑如何使用这里展示的代码以及 Prosite 中的信息。

James Tisdall 将于 2002 年 1 月 28 日至 31 日在亚利桑那州的图森参加 O'Reilly 首届生物信息学技术会议。更多详情请访问生物信息学会议网站。

位于 prosite.dat（或我们的小型测试文件 prosmall.dat）中的 Prosite 数据以“记录”形式组织，每个记录由几行组成，并始终包含一个 ID 行和包含“//”的终止行。Prosite 行均以两个字符的代码开始，指定该行上的数据类型。以下是来自 Prosite 用户手册中所有可能记录类型的详细说明。

** ID 标识（每个条目开始；每条目一个）

AC 访问号（每条目一个）

DT 日期（每条目一个）

DE 简短描述（每条目一个）

PA 模式（≥0 每条目）

MA 矩阵/轮廓（≥0 每条目）

RU 规则（≥0 每条目）

NR 数值结果（≥0 每条目）

CC 注释（≥0 每条目）

DR 与 SWISS-PROT 的交叉引用（≥0 每条目）

3D 与 PDB 的交叉引用（≥0 每条目）

DO 指向文档文件的指针（每条目一个）

// 终止行（结束每个条目；每条目一个）

这些行类型中的每种类型都有某些类型的信息，这些信息以特定方式格式化，详细说明在 Prosite 文档中。

Prosite 模式

让我们具体看看 Prosite 模式。这些模式以一种迷你语言的形式呈现，描述了一组简短的蛋白质片段，这些片段可能是已知生物活性的区域。以下是来自 Prosite 用户手册中的模式“语言”描述

PA（模式）行包含 Prosite 模式的定义。模式使用以下约定进行描述

使用标准的 IUPAC 单字母编码表示氨基酸。

符号 `x’ 表示可以接受任何氨基酸的位置。

通过方括号 `[ ]’ 列出给定位置可接受的可接受氨基酸来指示模糊性。例如：[ALT] 代表 Ala 或 Leu 或 Thr。

通过花括号 `{ }’ 列出给定位置不接受氨基酸来指示模糊性。例如：{AM} 代表除了 Ala 和 Met 之外的所有氨基酸。

模式中的每个元素之间由破折号 `-‘ 分隔。

可以使用跟在元素后面的数字或数字范围来表示模式中元素的重复。例如：x(3) 相当于 x-x-x，x(2,4) 相当于 x-x 或 x-x-x 或 x-x-x-x。

当模式限制在序列的 N 端或 C 端时，该模式以 `<` 或 `>` 符号开始或结束。

句号结束模式。

将 Prosite 模式转换为 Perl 正则表达式的 Perl 子程序

为了在Perl程序中使用这种模式数据，我们需要将Prosite模式转换为Perl正则表达式，这是在Perl中搜索数据中模式的主要方式。为了这篇文章，我将假设您了解基本正则表达式语法。（如果您不了解，请阅读程序注释，并跳过Perl正则表达式。）以下子例程的示例，它将Prosite模式[AC]-x-V-x(4)-{ED}.转换为等效的Perl正则表达式[AC].V.{4}[^ED]。

以下是我们的第一个Perl代码，子例程PROSITE_2_regexp，用于将Prosite模式转换为Perl正则表达式

#
# Calculate a Perl regular expression
#  from a PROSITE pattern
#
sub PROSITE_2_regexp {

  #
  # Collect the PROSITE pattern
  #
  my($pattern) = @_;

  #
  # Copy the pattern to a regular expression
  #
  my $regexp = $pattern;

  #
  # Now start translating the pattern to an
  #  equivalent regular expression
  #

  #
  # Remove the period at the end of the pattern
  #
  $regexp =~ s/.$//;

  #
  # Replace 'x' with a dot '.'
  #
  $regexp =~ s/x/./g;

  #
  # Leave an ambiguity such as '[ALT]' as is.
  #   However, there are two patterns [G>] that need
  #   special treatment (and the PROSITE documentation
  #   is a bit vague, perhaps).
  #
  $regexp =~ s/\[G\>\]/(G|\$)/;

  #
  # Ambiguities such as {AM} translate to [^AM].
  #
  $regexp =~ s/{([A-Z]+)}/[^$1]/g;

  #
  # Remove the '-' between elements in a pattern
  #
  $regexp =~ s/-//g;

  #
  # Repetitions such as x(3) translate as x{3}
  #
  $regexp =~ s/\((\d+)\)/{$1}/g;

  #
  # Repetitions such as x(2,4) translate as x{2,4}
  #
  $regexp =~ s/\((\d+,\d+)\)/{$1}/g;

  #
  # '<' "beginning # $regexp="~" ' '^' ; \< ^ becomes for of s sequence">' becomes '$' for "end of sequence"
  #
  $regexp =~ s/\>/\$/;

  #
  # Return the regular expression
  #
  return $regexp;
}

子例程PROSITE_2_regexp接受Prosite模式并将其各个部分逐步转换为等效的Perl正则表达式，如子例程注释中所述。如果您此时还不了解Perl正则表达式语法，请只阅读注释——即以#字符开始的行。这将让您大致了解子例程，即使您根本不了解Perl。

从O'Reilly的《精通正则表达式：Perl和其他工具的强大技术》了解更多关于正则表达式的信息。

Perl子例程，将Prosite记录解析为其行类型

我们需要完成的另一个任务是解析各种类型的行，以便我们可以轻松获取ID和PA模式行。下一个子例程完成此任务：给定一个Prosite记录，它返回一个哈希，其中每个类型的行都按其两字符“行类型”索引。我们将感兴趣的键是ID键，用于包含标识信息的行；以及PA键，用于包含模式信息的行。

这个“get_line_types”子例程比我们需要的要多。它对所有行类型建立哈希索引，而不仅仅是我们将在这里实际使用的ID和PA行。但那没关系。这个子例程很短且简单，我们可能以后会使用它来处理Prosite记录中的其他类型的行。

通过构建哈希来存储记录的行，我们可以通过给出行类型代码（如ID用于识别号码）来提取记录中我们喜欢的任何数据行。我们可以使用这个哈希来提取我们将感兴趣的两种行类型，即ID标识行和PA模式行。然后，通过将Prosite模式转换为Perl正则表达式（使用我们之前编写的第一个子例程），我们将能够真正地在蛋白质序列中查找所有模式。换句话说，我们将提取模式信息并将其用于我们的Perl程序中，以便在蛋白质序列中搜索模式。

如果您对学习Perl感兴趣，不要错过O'Reilly的热销图书《学习Perl，第3版》，该版本已更新以涵盖Perl 5.6版本，并重写以反映当今学习Perl的程序员的实际需求。有关O'Reilly关于Perl的所有图书的完整列表，请访问perl.oreilly.com。

以下是我们的第二个子例程，它接受一个Prosite记录，并返回一个哈希，其中记录的行按其行类型索引

#
# Parse a PROSITE record into "line types" hash
#
sub get_line_types {

  #
  # Collect the PROSITE record
  #
  my($record) = @_;

  #
  # Initialize the hash
  #   key   = line type
  #   value = lines
  #
  my %line_types_hash = ();

  #
  # Split the PROSITE record to an array of lines
  #
  my @records = split(/\n/,$record);

  #
  # Loop through the lines of the PROSITE record
  #
  foreach my $line (@records) {

    #
    # Extract the 2-character name
    # of the line type
    #
    my $line_type = substr($line,0,2);

    #
    # Append the line to the hash
    # indexed by this line type
    #
    (defined $line_types_hash{$line_type})
    ?  ($line_types_hash{$line_type} .= $line)
    :  ($line_types_hash{$line_type} = $line);
  }

  #
  # Return the hash
  #
  return %line_types_hash;
}

主程序

现在让我们看看代码是如何工作的。以下程序使用我们刚才定义的子例程逐个从平坦文件prosmall.txt中的数据库读取Prosite记录。然后它将不同类型的行分开（如“PA”用于模式），并使用我们之前编写的PROSITE_2_regexp子例程将模式转换为正则表达式。最后，它在蛋白质序列中搜索正则表达式，并报告匹配模式在序列中的位置。

#!/usr/bin/perl
#
# Parse patterns from the PROSITE database, and
# search for them in a protein sequence
#

#
# Turn on useful warnings and constraints
#
use strict;
use warnings;

#
# Declare variables
#

#
# The PROSITE database
#
my $prosite_file = 'prosmall.dat';

#
# A "handle" for the opened PROSITE file
#
my $prosite_filehandle;

#
# Store each PROSITE record that is read in
#
my $record = '';

#
# The protein sequence to search
# (use "join" and "qw" to keep line length short)
#
my $protein = join '', qw(
MNIDDKLEGLFLKCGGIDEMQSSRTMVVMGGVSG
QSTVSGELQDSVLQDRSMPHQEILAADEVLQESE
MRQQDMISHDELMVHEETVKNDEEQMETHERLPQ
);

#
# open the PROSITE database or exit the program
#
open($prosite_filehandle, $prosite_file)
 or die "Cannot open PROSITE file $prosite_file";

#
# set input separator to termination line //
#
$/ = "//\n";

#
# Loop through the PROSITE records
#
while($record = <$prosite_filehandle>) {

  #
  # Parse the PROSITE record into its "line types"
  #
  my %line_types = get_line_types($record);

  #
  # Skip records without an ID (the first record)
  #
  defined $line_types{'ID'} or next;

  #
  # Skip records that are not PATTERN
  # (such as MATRIX or RULE)
  #
  $line_types{'ID'} =~ /PATTERN/ or next;

  #
  # Get the ID of this record
  #
  my $id = $line_types{'ID'};
  $id =~ s/^ID   //;
  $id =~ s/; .*//;

  #
  # Get the PROSITE pattern from the PA line(s)
  #
  my $pattern = $line_types{'PA'};
  # Remove the PA line type tag(s)
  $pattern =~ s/PA   //g;

  #
  # Calculate the Perl regular expression
  # from the PROSITE pattern
  #
  my $regexp =  PROSITE_2_regexp($pattern);

  #
  # Find the PROSITE regular expression patterns
  # in the protein sequence, and report
  #
  while ($protein =~ /$regexp/g) {
    my $position = (pos $protein) - length($&) +1;
    print "Found $id at position $position\n";
    print "   match:   $&\n";
    print "   pattern: $pattern\n";
    print "   regexp:  $regexp\n\n";
  }
}

#
# Exit the program
#
exit;

此程序以文件 parse_prosite 的形式在线提供。小巧的示例 Prosite 数据库文件名为 prosmall.dat。如果您将这些文件保存在您的计算机（Unix、Linux、Macintosh 或 Windows）上，您可以在命令行提示符中输入以下命令（在保存这两个文件的同一文件夹中）：

% perl parse_prosite

它将产生以下输出：

Found PKC_PHOSPHO_SITE at position 22
   match:   SSR
   pattern: [ST]-x-[RK].
   regexp:  [ST].[RK]

Found PKC_PHOSPHO_SITE at position 86
   match:   TVK
   pattern: [ST]-x-[RK].
   regexp:  [ST].[RK]

Found CK2_PHOSPHO_SITE at position 76
   match:   SHDE
   pattern: [ST]-x(2)-[DE].
   regexp:  [ST].{2}[DE]

Found MYRISTYL at position 30
   match:   GGVSGQ
   pattern: G-{EDRKHPFYW}-x(2)-[STAGCN]-{P}.
   regexp:  G[^EDRKHPFYW].{2}[STAGCN][^P]

如您所见，我们的简短程序逐条记录地遍历 Prosite 数据库，根据其中的行类型解析每条记录。如果记录包含 ID 和模式，它会从中提取它们，将模式创建为 Perl 正则表达式，并在蛋白质序列中搜索该正则表达式，报告找到的模式。