開髮用於二硫鍵研究的生物信息學數據庫

蛋白質數據庫 (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;

Before any analysis, to avoid biases, sequences have to be checked for their level of identity. 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 的寶貴諮詢。