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

本文概述

• 输出如下
• 材料和方法
• NR格式
• 数据库构造和解析数据
• 过滤和处理数据
• PostgreSQL作为处理中介

蛋白质数据库(PDB)生物信息学数据库是世界上最大的由实验确定的蛋白质, 核酸和复杂装配体结构的存储库。所有数据都是使用实验方法收集的, 例如X射线, 光谱学, 晶体学, NMR等。

本文介绍了如何从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个字符, 其中标签占用六个或更少的字符, 外加一个或多个空格, 总共八个字节。在某些情况下, 标签和数据之间通常也没有空格, 通常用于CONECT标签。

标题

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行中存储的数据是实验中每个原子的坐标数据。有时, 实验具有插入, 突变, 替代位置或多个模型。这导致多次重复同一原子。选择正确的原子将在后面说明。

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线。

模型和ENDMDL

MODEL线标记结构模型的起始位置, 并包含模型的序列号。

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

债券

这些线在半胱氨酸氨基酸之间包含二硫键。二硫键可以存在于其他残基类型中, 但在本文中仅分析氨基酸, 因此仅包含半胱氨酸。以下示例来自代码为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-alpha和C端。这张图片没有显示二硫键, 但是简化了该蛋白质的空间排列：

通过将肽键合的原子与C-alpha原子相连来创建管道。半胱氨酸的颜色与CPK着色方法中的硫磺颜色相同。当胱氨酸足够接近时, 它们的硫会形成二硫键, 从而增强结构。否则, 该蛋白质将结合过多, 并且其结构在较高温度下将不稳定。

CONECT

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

资源

这些记录包含有关从中提取分子的来源生物的信息。它们包含子记录, 以便在分类法中更轻松地定位, 并且具有与标题记录相同的多行结构。

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, 则为。提到的列代表接受p值分别为10e-7、10e-40和10e-80的集合。仅将p值截止值为10e-7的集合用于分析。

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

#---------------------------------------------------------------------------------------------------------------------------
# 1  2       3     4    5 6     7    8 9     A    B C     D    E F    G      H      I      J     K     L   M   N     O P Q
#---------------------------------------------------------------------------------------------------------------------------
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和残渣_id。

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

RESIDUE   PRO     17
CONECT      N      3 CA   CD   H   
CONECT      CA     4 N    C    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

在该项目中, 编码速度比执行速度和存储消耗更为重要。我决定不进行标准化-毕竟, 我们的目标是生成一个包含简介中提到的列的表。

在这一部分中, 将仅解释最重要的表。

主要表是：

• 蛋白质：带有实验名称和代码的表格。
• ps：主要结构表, 其中将包含序列, chain_id和代码。
• ts：包含从原始数据提取并转换为ATOM记录格式的三级/四级结构的表。该表将用作登台表, 提取后可以删除。配体不包括在内。
• 来源：从中获得实验数据的生物列表。
• tax_names, taxonomy_path, taxonomy_paths：来自NCBI分类数据库的Linnean分类名称, 用于从来源中列出的生物获取分类路径。
• 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, 蛋白质, 来源, sources_organela和ss_bond表。

ps表具有三个重要列：链, 长度和序列。使用C-alpha原子为每个链和按残基序列排序的每个残基生成蛋白质序列, 仅采用第一个插入和第一个备用位置。链是从TS.chain列中获取的, 而length只是序列字符串的预先计算的长度。由于本文仅用于分析单链和链内连接, 因此此处的分析将跳过多链蛋白。

在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 between之间, 这是平均值(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和蛋白质表连接的nr表中, 并标记组(set7, set40和set80)。最后, 根据蛋白质数量仅分析一组。从分析中删除了PDB和NR之间不匹配的链。

没有二硫键的蛋白质以及不属于任何集合的蛋白质都被排除在研究之外。数据是用DSSP处理的, 并且这些文件的分辨率有问题或要处理的原子太多也被排除在外。尽管没有链间连接, 但可以通过计算具有不同链的连接数从ssbond表中轻松计算出来, 因此仅使用单链蛋白作为分析结果。

对于剩余的蛋白质, 将对键的总数和重叠键的数目进行更新, 并对每个集合进行更新。

源生物是从源记录中提取的。在未知, 合成, 设计, 人工或混合的情况下, 将其从研究中丢弃。低分辨率蛋白质也只有在其侧链不可见时才被排除。

SOURCE记录存储在源表中, 该表包含分类法行。在某些情况下, 分类法缺失或不正确。在这种情况下, 需要对专家进行人工纠正。

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

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

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

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

使用所有这些数据和手动更正(设置正确的生活域), 我们能够创建taxonomy_path表的结构。数据采样如下所示：

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

在进行任何分析之前, 为了避免偏差, 必须检查序列的同一性水平。尽管NR集包含已相互比较的序列, 但始终建议进行额外检查。

对于每个二硫键先前的统计分析, 如果数据是反应性的或重叠的, 则标记该数据。标记重叠后, 该信息会自动显示每个蛋白质内部有多少个连续键和非连续键, 并且数据存储在蛋白质表中, 从中排除了所有蛋白质复合物的最终结果。

通过检查两个键侧是否属于同一NR集, 还标记了每个二硫键与组的关联。如果不是这种情况, 则跳过二硫键进行该组分析。

要通过键的方差分析键的数量, 必须将每个长度放在特定的类别中。在这种情况下, 仅选择五个类, 如下面的函数所示。

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;

蛋白质的大多数大小都少于400个氨基酸, 因此, 通过将长度分为五个范围来完成长度分类, 每个范围包含100个氨基酸, 最后一个除外。你可以在下面看到如何使用此功能提取数据进行分析：

 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;

带有正确标题和一些手动添加列的示例结果如下所示。

PostgreSQL作为处理中介

在这项工作中, 我们展示了如何从获取到分析来处理数据。处理科学数据时, 有时需要规范化, 有时则不需要。当处理少量数据时, 这些数据将不会再用于其他分析, 那么足以在处理速度足够快的情况下将其标准化即可。

之所以在一个生物信息学数据库中全部完成, 是因为PostgreSQL能够集成多种语言。其中包括R, 它可以直接在数据库中进行统计分析, 这是以后有关生物信息学工具的文章的主题。

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