开发用于二硫键研究的生物信息学数据库

蛋白质数据库 (PDB) 生物信息学数据库是世界上最大的由实验确定的蛋白质、核酸和复杂组件结构的存储库。 所有数据均使用 X 射线、光谱学、晶体学、核磁共振等实验方法收集。

本文介绍如何从 PDB 中提取、过滤和清理数据。 反过来，这可以实现文章中解释的分析类型在不同生命领域中蛋白质二硫键的发生：蛋白质数据库中蛋白质的比较，发表于蛋白质工程、设计和选择，第 27 卷，第 3 期， 2014 年 3 月 1 日，第 65-72 页。

PDB 有许多具有不同分辨率、方法、突变等的重复结构。用相同或相似的蛋白质进行实验会在任何组分析中产生偏差，因此我们需要从任何一组重复中选择正确的结构. 为此，我们需要使用一组非冗余 (NR) 蛋白质。

出于规范化的目的，我建议下载用于将原子名称导入使用 3NF 或星型模式和维度建模的数据库的化合物字典。 （我还使用 DSSP 来帮助消除有问题的结构。我不会在本文中介绍它，但请注意，我没有使用任何其他 DSSP 功能。）

本研究中使用的数据包含含有至少一个来自不同物种的二硫键的单单元蛋白质。 为了进行分析，所有二硫键首先按结构域（古细菌、原核生物、病毒、真核生物或其他）以及长度分类为连续的或非连续的。

输出

为了准备好输入到 R、SPSS 或其他工具中，分析师需要将数据放在具有以下结构的数据库表中：

材料和方法

为了实现这一目标，第一步是从 rcsb.org 收集数据。 该站点包含各种格式的可下载实验的 PDB 结构。

尽管数据以多种格式存储，但在此示例中，仅使用格式化的固定空格分隔文本格式 (PDB)。 PDB 文本格式的替代方案是其 XML 版本 PDBML，但它有时包含格式错误的原子位置命名条目，这可能会导致数据分析出现问题。 较旧的 mmCIF 和其他格式也可能可用，但本文不会对其进行解释。

PDB 格式

PDB 格式是一种分段的固定宽度文本格式，例如，可以很容易地被 SQL 查询、Java 插件或 Perl 模块解析。 文件容器中的每种数据类型都表示为以相应标签开头的行——我们将在以下小节中介绍每种标签类型。 行长度小于或等于 80 个字符，其中一个标签占用 6 个或更少的字符加上一个或多个空格，总共占用 8 个字节。 也有标签和数据之间没有空格的情况，通常用于CONECT标签。

TITLE

TITLE标签将一行标记为实验的（部分）标题，包含分子名称和其他相关数据，如特定氨基酸的插入、突变或删除。

 12345678901234567890123456789012345678901234567890123456789012345678901234567890 TITLE A TWO DISULFIDE DERIVATIVE OF CHARYBDOTOXIN WITH DISULFIDE TITLE 2 13-33 REPLACED BY TWO ALPHA-AMINOBUTYRIC ACIDS, NMR, 30 TITLE 3 STRUCTURES

如果TITLE记录有多行，则必须连接标题，按连续编号排序，该编号右对齐放置在字节 9 和 10 上（取决于这些行的数量）。

ATOM

存储在ATOM行中的数据是实验中每个原子的坐标数据。 有时，一个实验有插入、突变、替代位置或多个模型。 这导致多次重复相同的原子。 选择正确的原子将在后面解释。

 12345678901234567890123456789012345678901234567890123456789012345678901234567890 ATOM 390 N GLY A 26 -1.120 -2.842 4.624 1.00 0.00 N ATOM 391 CA GLY A 26 -0.334 -2.509 3.469 1.00 0.00 C ATOM 392 C GLY A 26 0.682 -1.548 3.972 1.00 0.00 C ATOM 393 O GLY A 26 0.420 -0.786 4.898 1.00 0.00 O ATOM 394 H GLY A 26 -0.832 -2.438 5.489 1.00 0.00 H ATOM 395 HA2 GLY A 26 0.163 -3.399 3.111 1.00 0.00 H ATOM 396 HA3 GLY A 26 -0.955 -2.006 2.739 1.00 0.00 H

上面的例子取自实验1BAH 。 第一列标记记录的类型，第二列是原子的序列号。 结构中的每个原子都有自己的序列号。

序列号旁边是原子位置标签，占四个字节。 从该原子位置，可以提取元素的化学符号，该符号并不总是出现在其单独列的记录数据中。

在原子名称之后有一个三字母的残基代码。 在蛋白质的情况下，该残基对应于氨基酸。 接下来，用一个字母对链进行编码。 我们所说的链是指单链氨基酸，有或没有缺口，尽管有时配体可以分配给一条链； 这种情况可以通过氨基酸序列中非常大的间隙检测到，该序列在下一列中。 有时可以扫描包含突变的相同结构，在这种情况下，插入代码在序列列之后的额外列中可用。 插入代码包含一个字母，以帮助区分哪些残留物受到影响。

接下来的三列是每个原子的空间坐标，以埃 (Å) 为单位。 这些坐标旁边是占用列，它说明了原子在那个位置的概率是多少，通常的比例是从零到一。

倒数第二列是温度因子，它携带有关晶体无序的信息，以 Ų 为单位。 大于 60Ų 的值表示混乱，而小于 30Ų 的值表示自信。 它并不总是存在于 PDB 文件中，因为它取决于实验方法。

接下来的列——符号和费用——通常是缺失的。 如上所述，化学符号可以从原子位置列中收集。 当电荷存在时，它作为一个整数后缀在符号后面，后跟+或- ，例如N1+ 。

TER

这标志着链条的末端。 即使没有这条线，也很容易区分链的结束位置。 因此，经常缺少TER线。

MODEL和ENDMDL

MODEL行标记结构模型的开始位置，它包含模型的序列号。

在该模型中的所有原子行之后，它以ENDMDL行结束。

SSBOND

这些系包含半胱氨酸氨基酸之间的二硫键。 二硫键可以存在于其他残基类型中，但在本文中仅分析氨基酸，因此仅包括半胱氨酸。 以下示例来自代码为132L的实验：

 12345678901234567890123456789012345678901234567890123456789012345678901234567890 SSBOND 1 CYS A 6 CYS A 127 1555 1555 2.01 SSBOND 2 CYS A 30 CYS A 115 1555 1555 2.05 SSBOND 3 CYS A 64 CYS A 80 1555 1555 2.02 SSBOND 4 CYS A 76 CYS A 94 1555 1555 2.02

在此示例中，文件中标记了四个二硫键，它们的序列号在第二列中。 所有这些键都使用半胱氨酸（第 3 列和第 6 列），并且它们都存在于链A中（第 4 列和第 7 列）。 在每条链之后有一个残基序列号，表明该键在肽链中的位置。 插入代码位于每个残基序列旁边，但在此示例中它们不存在，因为该区域中没有插入氨基酸。 最后一列之前的两列是为对称操作保留的，最后一列是硫原子之间的距离，以 Å 为单位。

让我们花点时间为这些数据提供一些背景信息。

下面的图片是使用 rcsb.org NGL 查看器拍摄的，显示了实验132L的结构。 特别是，它们显示了一种没有配体的蛋白质。 第一张图片使用棒状表示，CPK 着色以黄色显示硫磺及其键。 V形硫连接代表蛋氨酸连接，而Z形连接是半胱氨酸之间的二硫键。

在下一张图片中，一种称为骨架可视化的简化蛋白质可视化方法用氨基酸着色，其中半胱氨酸是黄色的。 它代表相同的蛋白质，但不包括其侧链，仅包括其肽组的一部分——在这种情况下，是蛋白质骨架。 它由三个原子组成：N 端、C α 和 C 端。 这张图片没有显示二硫键，但它被简化以显示蛋白质的空间排列：

管道是通过将肽键合原子与 C-α 原子连接起来而形成的。 半胱氨酸的颜色与 CPK 着色法中硫的颜色相同。 当半胱氨酸足够接近时，它们的硫会产生二硫键，从而加强结构。 否则，这种蛋白质会结合太多，并且其结构在较高温度下会不太稳定。

CONECT

这些记录用于标记原子之间的连接。 有时这些标签根本不存在，而其他时候所有数据都被填写。在分析二硫键的情况下，这部分可能很有用 - 但不是必需的。 这是因为，在这个项目中，未标记的债券是通过测量距离添加的，所以这将是开销，并且还必须进行检查。

SOURCE

这些记录包含有关从中提取分子的来源生物体的信息。 它们包含子记录以便更容易地在分类中定位，并且具有我们在标题记录中看到的相同的多行结构。

 SOURCE MOL_ID: 1; SOURCE 2 ORGANISM_SCIENTIFIC: ANOPHELES GAMBIAE; SOURCE 3 ORGANISM_COMMON: AFRICAN MALARIA MOSQUITO; SOURCE 4 ORGANISM_TAXID: 7165; SOURCE 5 GENE: GST1-6; SOURCE 6 EXPRESSION_SYSTEM: ESCHERICHIA COLI; SOURCE 7 EXPRESSION_SYSTEM_TAXID: 562; SOURCE 8 EXPRESSION_SYSTEM_STRAIN: BL21(DE3)PLYSS; SOURCE 9 EXPRESSION_SYSTEM_VECTOR_TYPE: PLASMID; SOURCE 10 EXPRESSION_SYSTEM_PLASMID: PXAGGST1-6

NR 格式

这是非冗余 (NR) 链 PDB 集的列表。 它的快照可以在 ftp.ncbi.nih.gov/mmdb/nrtable/ 找到。 其目的是避免由蛋白质相似性引起的不必要的偏差。 NR 具有通过比较所有 PDB 结构创建的具有不同身份 p 值水平的三个集合。 结果被添加到文本文件中，稍后将对其进行解释。 此项目并非需要所有列，因此将仅解释重要的列。

前两列包含唯一的 PDB 实验代码和链标识符，如上述ATOM记录所述。 第 6、9 和 C 列包含有关 p 值代表性的信息，这是由 BLAST 计算的序列的相似性水平。 如果该值为零，则不接受它是集合的一部分； 如果值为 1，则为 1。 提到的列表示接受 p 值截止值分别为 10e-7、10e-40 和 10e-80 的集合。 只有 p 值截止值为 10e-7 的集合将用于分析。

最后一列包含有关结构可接受性的信息，其中a可接受而n不可接受。

 #--------------------------------------------------------------------------------------------------------------------------- # 1 2 3 4 5 6 7 8 9 ABCDEFGHIJKLMNOPQ #--------------------------------------------------------------------------------------------------------------------------- 3F8V A 69715 1 1 1 1 1 1 1 1 1 9427 1 1 0.00 0.00 0.00 0.00 1.08 1 6 5 164 X a 3DKE X 68132 1 2 0 1 2 0 1 2 0 39139 1 1 0.00 0.00 0.00 0.00 1.25 1 11 7 164 X a 3HH3 A 77317 1 3 0 1 3 0 1 3 0 90 1 0 0.00 0.00 0.00 0.00 1.25 1 5 4 164 X a 3HH5 A 77319 1 4 0 1 4 0 1 4 0 90 2 0 0.00 0.00 0.00 0.00 1.25 1 4 4 164 X a

数据库构建和解析数据

现在我们已经知道我们正在处理什么以及我们需要做什么，让我们开始吧。

下载数据

此分析的所有数据都可以在以下三个地址找到：

• ftp.wwpdb.org/pub/pdb/data/structures/divided/pdb/
• ftp.ncbi.nih.gov/mmdb/nrtable/
• ftp.ncbi.nih.gov/pub/taxonomy/taxdmp.zip

前两个链接包含档案列表。 应使用每个文件的最新存档来避免因缺乏分辨率或质量而引起的问题。 第三个链接直接包含最新的分类档案。

解析数据

通常 PDB 文件的解析是由 Java、Perl 或 Python 中的插件或模块完成的。 在本研究中，我编写了一个自定义 Perl 应用程序，而没有使用预先编写的 PDB 解析模块。 原因是在解析大量数据时，根据我的经验，使用实验数据最常见的问题是数据中的错误。 有时，坐标、距离、线长、不应该出现的地方的注释等会出现错误。

处理这个问题的最有效方法是最初将数据库中的所有内容都存储为原始文本。 编写通用解析器来处理完全符合规范的理想数据。 但实际上，数据并不理想，这将在过滤部分进行解释，您将在其中找到 Perl 导入脚本。

数据库建设

在构建数据库时，请注意该数据库是为处理数据而构建的。 稍后的分析将在 SPSS 或 R 中完成。出于我们的目的，这里建议使用至少 8.4 版的 PostgreSQL。

表结构直接从下载的文件中复制而来，只做了一些小改动。 在这种情况下，记录的数量太少，不值得花时间进行规范化。 如前所述，此数据库仅供一次性使用：这些表不是为在网站上提供服务而构建的——它们只是用于处理数据。 一旦完成，它们可以被删除，或者作为补充数据备份，也许是为了让其他研究人员重复这个过程。

在这种情况下，最终结果将是一个表格，然后可以将其导出到一个文件中，以便与 SPSS 或 R 等统计工具一起使用。

表

从ATOM记录中提取的数据必须连接到HEADER或TITLE记录。 数据层次结构如下图所示。

由于这张图片是第三范式 (3NF) 中数据库的简化表示，就我们的目的而言，它包含太多开销。 原因：要计算用于二硫键检测的原子之间的距离，我们需要进行连接。 在这种情况下，我们将有一个表连接到自身两次，并且每次连接到二级结构和一级结构两次，这是一个非常缓慢的过程。 由于并非每个分析都需要二级结构信息，因此如果您需要重用数据或分析更大量的二硫键，建议使用另一种模式：

二硫键不像其他共价键那样频繁，因此不需要仓库模型，尽管可以使用它。 下面的星型模式和维度建模将花费太多时间来开发，并且会使查询更加复杂：

在必须处理所有债券的情况下，我推荐星型模式。

（否则不需要，因为二硫键不像其他键那么常见。在这项工作的情况下，二硫键的数量接近 30,000，这在 3NF 中可能足够快，但我决定通过一个非规范化的表，所以这里没有画出来。）

所有共价键的预期总数至少是三级结构中原子数的两倍，在这种情况下，3NF 会非常慢，因此需要使用星型模式形式进行非规范化。 在该模式中，一些表有两个外键检查，这是因为在两个原子之间创建了一个键，所以每个原子需要有自己的primary_structure_id 、 atom_name_id和residue_id 。

填充d_atom_name维表有两种方法：从数据中，从另一个来源，即我前面提到的化学成分字典。 它的格式类似于 PDB 格式：只有RESIDUE和CONECT行有用。 这是因为RESIDUE的第一列包含一个残基三字母代码，而CONECT包含原子的名称及其连接，它们也是原子名称。 因此，从这个文件中，我们可以解析所有原子名称并将它们包含在我们的数据库中，尽管我建议您允许数据​​库包含未列出的原子名称的可能性。

 RESIDUE PRO 17 CONECT N 3 CA CD H CONECT CA 4 NC CB HA CONECT C 3 CA O OXT CONECT O 1 C CONECT CB 4 CA CG HB2 HB3 CONECT CG 4 CB CD HG2 HG3 CONECT CD 4 N CG HD2 HD3 CONECT OXT 2 C HXT CONECT H 1 N CONECT HA 1 CA CONECT HB2 1 CB CONECT HB3 1 CB CONECT HG2 1 CG CONECT HG3 1 CG CONECT HD2 1 CD CONECT HD3 1 CD CONECT HXT 1 OXT END HET PRO 17 HETNAM PRO PROLINE FORMUL PRO C5 H9 N1 O2

在这个项目中，编码速度比执行速度和存储消耗更重要。 我决定不规范化——毕竟，我们的目标是生成一个包含介绍中提到的列的表。

在这一部分中，将只解释最重要的表。

主要表格有：

• proteins ：带有实验名称和代码的表。
• ps ：包含sequence 、 chain_id和code的主结构表。
• ts ：包含从原始数据中提取并转换为ATOM记录格式的三级/四级结构的表。 这将用作临时表，并且可以在提取后删除。 不包括配体。
• sources ：从中得出实验数据的生物体列表。
• tax_names , taxonomy_path , taxonomy_paths ：来自 NCBI 分类数据库的林奈分类名称，用于从sources中列出的生物中获取分类路径。
• nr ：从 NR 集中提取的 NCBI 非冗余蛋白质列表。
• pdb_ssbond ：给定 PDB 文件中的二硫键列表。

过滤和处理数据

从 RCSB PDB 存储库的快照中检索数据。

每个文件都使用 Perl 脚本导入到我们的 PostgreSQL 数据库中的单个表raw_pdb中。 该脚本使用每块 10,000 次插入的事务。

raw_pdb的结构是这样的：

导入脚本如下所示：

 #!/usr/bin/perl -w use FindBin '$Bin'; use DBI; $dbName = 'bioinf'; $dbLogin = 'ezop'; $dbPass = 'XYZ'; $conn = DBI->connect("DBI:Pg:database=$dbName;host=localhost", "$dbLogin", "$dbPass", {'RaiseError' => 1, 'AutoCommit' => 0}); die "./pdb_lines_unos.pl <table> <file>" if not defined($ARGV[0]); $recordCount = 0; $table = $ARGV[0]; $fName = $ARGV[1]; open F, "zcat $fName|"; while (<F>) { chomp; $linija = $_; $recordCount += 1; $sql = "insert into $table (code, line_num, line_cont) values (?, ?, ?)"; $conn->do($sql, undef, $fName, $recordCount, $linija); $conn->commit() if ($recordCount%10000 == 0); } close F; $conn->commit(); 1;

导入行后，将使用我们将在下面定义的函数对其进行解析。

根据raw_pdb数据，我们通过解析相应的记录生成表ts 、 ps 、 proteins 、 sources 、 sources_organela和ss_bond 。

ps表具有三个重要的列： chain 、 length和sequence 。 使用 C-α 原子生成蛋白质序列，用于每个链和按残基序列排序的每个残基，仅采用第一个插入和第一个替代位置。 chain取自TS.chain列， length只是sequence字符串的预先计算的长度。 由于本文仅分析单链和链内连接，因此我们在此处的分析中跳过了多链蛋白质。

在SSBOND记录中，所有二硫键都存储在pdb_ssbond表中，该表继承自pdb_ssbond_extended表。 pdb_ssbond看起来像这样：

我还添加了这些索引：

 "pdb_ssbond_pkey" PRIMARY KEY, btree (id) "ndxcode1" btree (code, chain_id1, res_seq1) "ndxcode2" btree (code, chain_id2, res_seq2)

假设截止前二硫键的分布是高斯分布（未使用例如 KS-test 进行测试），因此计算同一蛋白质中半胱氨酸之间的每个距离的标准偏差，以发现允许的键长范围并进行比较到截止。 截止值与计算的平均值加减三个标准差相同。 我们扩展了范围以引入更多可能的二硫键，这些二硫键未在SSBOND行的 PDB 文件中列出，但我们通过计算ATOM记录之间的距离自行插入。 ssbonds的选择范围在 1.6175344456264 和 2.48801947642267 Å 之间，即平均值 (2.05) 加减四个标准差：

 select count(1) as cnt , stddev(dist) as std_dev , avg(dist) as avg_val , -stddev(dist) + avg(dist) as left1 , stddev(dist) + avg(dist) as right1 , -2 * stddev(dist) + avg(dist) as left2 , 2 * stddev(dist) + avg(dist) as right2 , -3 * stddev(dist) + avg(dist) as left3 , 3 * stddev(dist) + avg(dist) as right3 , -4 * stddev(dist) + avg(dist) as left4 , 4 * stddev(dist) + avg(dist) as right4 , min(dist) , max(dist) from pdb_ssbond_tmp where dist > 0 and dist < 20;

TS表包含所有原子的坐标，但只使用半胱氨酸，它们的硫命名为" SG " 。 因此，创建了另一个仅包含" SG "硫原子的暂存表，以通过减少要搜索的记录数来加快处理速度。 仅选择硫时，组合的数量远远少于所有原子的情况——前者为 194,574，而后者为 122,761,100。 在此连接到自身的表中，使用欧几里德距离计算距离，并将结果导入到pdb_ssbond表中，但仅在距离介于之前计算的定义长度之间的情况下。 进行这种加速的原因是为了减少再次运行整个过程的时间——出于检查目的——将其保持在一天之内。

为了清除二硫键，我们使用以下算法：

• 当连接两边指向同一个氨基酸时删除
• 删除长度不在 1.6175344456264 和 2.48801947642267 之间的键
• 删除插入
• 去除由交替原子位置引起的键，但保留第一个位置

代码（以pdb_ssbond表作为第一个参数）是：

 CREATE OR REPLACE FUNCTION pdb_ssbond_clean2( clean_icodes boolean, clean_altloc boolean, mark_reactive boolean, mark_consecutive boolean) RETURNS void AS $$ declare _res integer; BEGIN delete from pdb_ssbond b where exists ( select a.id from pdb_ssbond a where a.code = b.code and a.id > b.id and ( ( a.chain_id1 = b.chain_id1 and a.res_seq1 = b.res_seq1 and a.chain_id2 = b.chain_id2 and a.res_seq2 = b.res_seq2) or ( a.chain_id1 = b.chain_id2 and a.res_seq1 = b.res_seq2 and a.chain_id2 = b.chain_id1 and a.res_seq2 = b.res_seq1) ) ) ; delete from pdb_ssbond where chain_id1 = chain_id2 and res_seq1 = res_seq2 and i_code1 = i_code2; update pdb_ssbond set is_consecutive = true , is_reactive = false; delete from pdb_ssbond where not dist between 1.6175344456264 and 2.48801947642267; if clean_icodes then delete from pdb_ssbond where i_code1 not in ('', ' ', 'A') or i_code2 not in ('', ' ', 'A') ; end if; if clean_altloc then delete from pdb_ssbond a where exists ( select 1 from pdb_ssbond b where b.code = a.code and b.chain_id1 = a.chain_id1 and b.res_seq1 = a.res_seq1 and b.i_code1 = a.i_code1 and b.ser_num < a.ser_num ) ; delete from pdb_ssbond a where exists ( select 1 from pdb_ssbond b where b.code = a.code and b.chain_id2 = a.chain_id2 and b.res_seq2 = a.res_seq2 and b.i_code2 = a.i_code2 and b.ser_num < a.ser_num ) ; end if; if mark_reactive then update pdb_ssbond set is_reactive = false ; update pdb_ssbond set is_reactive = true where exists ( select 1 from pdb_ssbond b where b.code = pdb_ssbond.code and b.chain_id1 = pdb_ssbond.chain_id1 and b.res_seq1 = pdb_ssbond.res_seq1 and b.i_code1 = pdb_ssbond.i_code1 and b.ser_num < pdb_ssbond.ser_num ) ; update pdb_ssbond set is_reactive = true where exists ( select 1 from pdb_ssbond b where b.code = pdb_ssbond.code and b.chain_id2 = pdb_ssbond.chain_id2 and b.res_seq2 = pdb_ssbond.res_seq2 and b.i_code2 = pdb_ssbond.i_code2 and b.ser_num < pdb_ssbond.ser_num ) ; update pdb_ssbond set is_reactive = true where (code, chain_id1, res_seq1, i_code1) in ( select code, chain_id1 as c, res_seq1 as r, i_code1 as i from pdb_ssbond intersect select code, chain_id2 as c, res_seq2 as r, i_code2 as i from pdb_ssbond ) ; update pdb_ssbond set is_reactive = true where (code, chain_id2, res_seq2, i_code2) in ( select code, chain_id1 as c, res_seq1 as r, i_code1 as i from pdb_ssbond intersect select code, chain_id2 as c, res_seq2 as r, i_code2 as i from pdb_ssbond ); end if; if mark_consecutive then update pdb_ssbond set is_consecutive = false; update pdb_ssbond set is_consecutive = true where not exists ( select 1 from pdb_ssbond a where a.code = pdb_ssbond.code and ( (a.chain_id1 = pdb_ssbond.chain_id1 and a.chain_id2 = pdb_ssbond.chain_id2) or (a.chain_id1 = pdb_ssbond.chain_id2 and a.chain_id2 = pdb_ssbond.chain_id1) ) and ( (a.res_seq1 between pdb_ssbond.res_seq1 and pdb_ssbond.res_seq2) or (a.res_seq2 between pdb_ssbond.res_seq1 and pdb_ssbond.res_seq2) or (pdb_ssbond.res_seq1 between a.res_seq1 and a.res_seq2) or (pdb_ssbond.res_seq2 between a.res_seq1 and a.res_seq2) ) and not ( a.code = pdb_ssbond.code and a.chain_id1 = pdb_ssbond.chain_id1 and a.chain_id2 = pdb_ssbond.chain_id2 and a.res_seq1 = pdb_ssbond.res_seq1 and a.res_seq2 = pdb_ssbond.res_seq2 ) ); end if; END; $$ LANGUAGE plpgsql;

这样，非冗余的蛋白质集被导入到连接到ps和proteins表的nr表中，并标记集（ set7 、 set40和set80 ）。 最后，根据蛋白质数量只分析一组。 从分析中删除 PDB 和 NR 之间的不匹配链。

没有二硫键的蛋白质与不属于任何一组的蛋白质一起被排除在研究之外。 数据用 DSSP 处理，这些文件有分辨率问题或需要处理的原子过多的文件也被排除在外。 由于没有保留链间连接，因此仅使用具有单链的蛋白质作为分析结果，尽管它们很容易通过计算具有不同链的连接数从ssbond表中计算出来。

对于剩余的蛋白质，对键的总数和重叠键的数量进行更新，并且对每一组进行更新。

源生物是从SOURCE记录中提取的。 如果它是未知的、合成的、设计的、人工的或混合的，它就会被从研究中丢弃。 只有当它们的侧链不可见时，才会排除低分辨率蛋白质。

SOURCE记录存储在sources表中，其中包含分类行。 在某些情况下，分类法缺失或不正确。 在这些情况下，需要专家的手动更正。

从从 NCBI 下载的源代码和分类法中，该类被分配给每个主要结构。 如果无法分配类别，则将蛋白质从分析列表中删除。 知道正在使用生物数据库，建议由生物学家对所有源记录和 NCBI 分类类别进行额外检查，否则不同域之间的分类可能会出现问题。 为了将源蜂窝位置与分类 ID 匹配，源表中的数据被提取到表sources_organela中，其中所有数据都以代码、标签和值的形式写入。 其格式如下图：

 select * from sources_organela where code = '1rav';

我们要使用的分类档案是一个包含七个转储文件的 ZIP 文件。 在这些文件中，最重要的两个是names.dmp和merged.dmp 。 这两个文件都是 CSV 制表符分隔的文件，详见文档：

• 文件merged.dmp包含以前分类ID 的列表，以及每个分类ID 合并到的当前分类ID。
• names.dmp按以下顺序包含这些重要列：
  • tax_id : 分类 ID
  • name_txt ：物种的名称，如果适用，唯一的名称（可以找到具有多个名称的物种，这有助于消除歧义）
• division.dmp包含我们将用作类的顶级域的名称。
• nodes.dmp是包含使用分类 ID 的有机体层次结构的表。
  • 它包含一个父分类 ID，一个可以在names.dmp中找到的外键。
  • 它还包含一个部门 ID，这对于我们正确存储相关的顶级域数据很重要。

通过所有这些数据和手动更正（设置正确的生活领域），我们能够创建taxonomy_path表的结构。 数据样本如下所示：

 select * from taxonomy_path order by length(path) limit 20 offset 2000;

在进行任何分析之前，为了避免偏差，必须检查序列的同一性水平。 Although the NR set contains sequences which are already compared between each other, an extra check is always recommended.

For each disulfide bond's prior statistical analysis, data is marked if it is reactive or overlapping. After marking overlaps, that information automatically reveals how many consecutive and non-consecutive bonds are inside each protein, and that data is stored in the proteins table from which all protein complexes are excluded in final result.

Each disulfide bond is marked also for its association to sets, by checking both bond sides to see if they belong to the same NR set. Where that is not the case, the disulfide bond is skipped for that set analysis.

To analyze the quantity of bonds by their variance, each length has to be put in a specific class. In this case, only five classes are chosen as written in the function below.

 CREATE OR REPLACE FUNCTION len2class(len integer) RETURNS integer AS $BODY$ BEGIN return case when len <= 100 then 1 when len > 100 and len <= 200 then 2 when len > 200 and len <= 300 then 3 when len > 300 and len <= 400 then 4 else 5 end; END; $BODY$ LANGUAGE plpgsql VOLATILE COST 100;

Most of the protein sizes are less than 400 amino acids, so length classification is done by splitting lengths into five ranges, each taking 100 amino acids, except the last one which takes the rest. Below you can see how to use this function to extract data for analysis:

 SELECT p.code, p.title, p.ss_bonds, p.ssbonds_overlap, p.intra_count, p.sci_name_src, p.len, p.tax_path, p.pfam_families, p.src_class, ( SELECT s.id FROM src_classes s WHERE s.src_class::text = p.src_class::text) AS src_class_id, p.len_class7, ( SELECT s.val FROM sources_organela s WHERE s.code = p.code::bpchar AND s.tag::text = 'EXPRESSION_SYSTEM_CELLULAR_LOCATION'::text) AS expression_system_cellular_location, ( SELECT s.val FROM sources_organela s WHERE s.code = p.code::bpchar AND s.tag::text = 'CELLULAR_LOCATION'::text) AS cellular_location, ps.sequence, ps.uniprot_code, ps.accession_code, ps.cc, ps.seq_uniprot, ps.chain_id FROM proteins p JOIN nr n ON n.code::text = p.code::text AND n.rep7 = true JOIN ps ps ON ps.code::text = n.code::text AND ps.chain_id = n.chain_id AND ps.het = false WHERE p.is_excluded = false AND p.chain_cnt = 1 AND p.is_set7 = true AND p.reactive_cnt7 = 0 ORDER BY p.code;

An example result with corrected titles and some manually added columns is shown below.

PostgreSQL 作为处理中介

在这项工作中，我们展示了如何处理数据，从获取到分析。 在处理科学数据时，有时需要标准化，有时则不需要。 当处理不会重复用于其他分析的少量数据时，在处理速度足够快的情况下将其非规范化就足够了。

这一切都在一个生物信息学数据库中完成的原因是 PostgreSQL 能够集成多种语言。 这包括 R，它可以直接在数据库中进行统计分析——这是未来关于生物信息学工具的文章的主题。

特别感谢 Toptal 的同事 Stefan Fuchs 和 Aldo Zelen 的宝贵咨询。