はじめに
RCSB Protein Data Bank (RCSB PDB) には非常に多くのタンパク質の構造が登録されていますが、稀に残基番号がずれている場合があります。そのようなPDBファイルの場合、これに紐付いている論文内では残基番号に矛盾が無いのですが、他の論文に書かれている残基番号とは異なっていて面倒です。特にタンパク質に変異を入れてリガンドとの活性を色々と測定している論文などを読むときは、同時にタンパク質の立体構造も確認するので、番号がずれているとイライラします。なぜ残基番号がずれるのか理由はよく分かりませんが、上記のような場合に修正したいと思うのではないでしょうか。今回はPDBの残基番号修正するプログラムをOpenBabelを利用して作成したいと思います。
計算対象など
今回はPDB codeが1LKEのものを例にします。今まで私が扱って来た残基番号のズレたPDBがどれなのか探したのですが、残基番号がズレているのは稀なので、これはネットで色々と検索してやっと見つけたタンパク質構造(PDB)です。下図はPyMOLで1LKEのPDBファイルを読み込んだときの画面です。画面の上部がアミノ酸配列が一文字略号で表示されてますが、アスパラギン酸が5番目から始まってます。
しかしUniProtでアミノ酸配列を調べてみると、5番目はイソロイシンであることがわかります。
実はこのアスパラギン酸はUniProtのアミノ酸配列的には20番目に相当します。
計算の方針
このようなPDBの残基番号のズレを修正するプログラムを作成するにあたり、以下のようなステップを考えます。
- PDBファイル中のアミノ酸配列をFASTAファイルに出力し、UniProtから基準となるアミノ酸配列のFASTAファイルを取得する。
- 2つのアミノ酸配列をClustal Omegaを利用してアライメントする。
- アライメント済みのFASTAファイルとPDBファイルを読み込む。
- FASTAファイルに従ってPDBファイル中の残基番号を修正する。
- 修正後のPDBファイルを出力する。
上記手順のうち1と2はプログラムで自動化せずにマニュアルで実行することとします。したがってプログラムを作成するのは3以降となります。PDBファイルはOpenBabelのライブラリを利用すれば簡単に読み込むプログラムを作成できます。さらにOpenBabelはFASTAファイルにも対応しています。しかしながらアライメント済みのFASTAファイル中にはギャップを示す”-“記号が含まれますが、OpenBabelはギャップ記号には対応していないようなので、FASTAファイルを取り扱う部分については自作することとします。前編ではFASTAファイルを取り扱う部分について見ていくこととしましょう。
アミノ酸配列のアライメント
先ずはPDBファイル中のアミノ酸配列とUniProtからのアミノ酸配列をアライメントします。以下のようなOpenBabelのコマンドでPDBファイルからアミノ酸配列を取り出します。
$ obabel 1lke.pdb -O 1lke.fasta出力されたFASTAファイルは次のようになります。
>1lke.pdb 216 bp; generated with OpenBabel 3.1.0
DGACPEVKPVDNFDWSQYHGKWWQVAAYPDHITKYGKCGWAEYTPEGKSVKVSRYSVIHG
KEYFSEGTAYPVGDSKIGKIYHSYTIGGVTQEGVFNVLSTDNKNYIIGYFCSYGHMDLVW
VLSRSMVLTGEAKTAVENYLIGSPVVDSQKLVYSDFSXそしてUniProtからFASTAファイルをダウンロードします。このFASTAファイル（ファイル名”P09464.fasta”）の内容は次の通りです。
>sp|P09464|BBP_PIEBR Bilin-binding protein OS=Pieris brassicae OX=7116 PE=1 SV=2
MQYLIVLALVAAASANVYHDGACPEVKPVDNFDWSNYHGKWWEVAKYPNSVEKYGKCGWA
EYTPEGKSVKVSNYHVIHGKEYFIEGTAYPVGDSKIGKIYHKLTYGGVTKENVFNVLSTD
NKNYIIGYYCKYDEDKKGHQDFVWVLSRSKVLTGEAKTAVENYLIGSPVVDSQKLVYSDF
SEAACKVNN2つのFASTAファイルを用意できたら結合します。
$ cat 1lke.fasta P09464.fasta > All.fasta2つの配列を結合したこのファイルをClustal Omegaを利用してアライメントします。
$ clustalo -i All.fasta -o aligned.fastaちなみにClustal OmegaはAnacondaを利用すると、次のコマンドにより簡単にインストールできます。
$ conda install -c bioconda clustalo最後にClustal OmegaによってアライメントされたFASTAファイルは次のような無いようになります。上の方がPDBファイルのアミノ酸配列で、所々にギャップが入ることによってUniProtのアミノ酸配列を基準とした残基番号に修正できます。
$ cat aligned.fasta
>1lke.pdb 216 bp; generated with OpenBabel 3.1.0
-------------------DGACPEVKPVDNFDWSQYHGKWWQVAAYPDHITKYGKCGWA
EYTPEGKSVKVSRYSVIHGKEYFSEGTAYPVGDSKIGKIYHSYTIGGVTQEGVFNVLSTD
NKNYIIGYFCSY-----GHMDLVWVLSRSMVLTGEAKTAVENYLIGSPVVDSQKLVYSDF
SX-------
>sp|P09464|BBP_PIEBR Bilin-binding protein OS=Pieris brassicae OX=7116 PE=1 SV=2
MQYLIVLALVAAASANVYHDGACPEVKPVDNFDWSNYHGKWWEVAKYPNSVEKYGKCGWA
EYTPEGKSVKVSNYHVIHGKEYFIEGTAYPVGDSKIGKIYHKLTYGGVTKENVFNVLSTD
NKNYIIGYYCKYDEDKKGHQDFVWVLSRSKVLTGEAKTAVENYLIGSPVVDSQKLVYSDF
SEAACKVNNプログラム解説
今回のブログで紹介するプログラム（プログラム名をread_fasta.cppとします）はclustaloコマンドから出力されたFASTAファイル（aligned.fasta）から最初の配列を読み込んで、そのままFASTA形式で出力するというものです。このプログラムのコンパイルと実行方法は以下の通りです。
$ g++ read_fasta.cpp
$ ./a.out
>1lke.pdb 216 bp; generated with OpenBabel 3.1.0
-------------------DGACPEVKPVDNFDWSQYHGKWWQVAAYPDHITKYGKCGWA
EYTPEGKSVKVSRYSVIHGKEYFSEGTAYPVGDSKIGKIYHSYTIGGVTQEGVFNVLSTD
NKNYIIGYFCSY-----GHMDLVWVLSRSMVLTGEAKTAVENYLIGSPVVDSQKLVYSDF
SX-------それではプログラムを見ていきましょう。
C++のソースコードは以下の通りです。
1 #include <iostream>
2 #include <fstream>
3 #include <sstream>
4 #include <string>
5 #include <vector>
6 #include <map>
7
8 using namespace std;
9
10 // 1文字を3文字に変換するための変数設定
11 map<char, string> aa3;
12
13 void init(void) {
14 aa3['A'] = "ALA";
15 aa3['R'] = "ARG";
16 aa3['N'] = "ASN";
17 aa3['D'] = "ASP";
18 aa3['C'] = "CYS";
19 aa3['E'] = "GLU";
20 aa3['Q'] = "GLN";
21 aa3['G'] = "GLY";
22 aa3['H'] = "HIS";
23 aa3['I'] = "ILE";
24 aa3['L'] = "LEU";
25 aa3['K'] = "LYS";
26 aa3['M'] = "MET";
27 aa3['F'] = "PHE";
28 aa3['P'] = "PRO";
29 aa3['S'] = "SER";
30 aa3['T'] = "THR";
31 aa3['W'] = "TRP";
32 aa3['Y'] = "TYR";
33 aa3['V'] = "VAL";
34 aa3['-'] = "GAP";
35 aa3['X'] = "ANY";
36 }
37
38 class SEQ {
39 string name;
40 string seq;
41 vector<pair<string, int> > seq3l;
42 int pos;
43
44 public:
45 SEQ() : pos(0) {} // コンストラクタ
46
47 SEQ(const SEQ& Obj) { // コピーコンストラクタ
48 name = Obj.name;
49 seq = Obj.seq;
50 seq3l = Obj.seq3l;
51 pos = Obj.pos;
52 }
53
54 ~SEQ() {} // デストラクタ
55
56 bool ReadFASTA(istream* ist) {
57 bool _f = false;
58 string line;
59 pos = 0;
60
61 name.clear();
62 while (getline(*ist, line)) {
63 if (line[0] == '>') {
64 name = line.substr(1);
65 _f = true;
66 break;
67 }
68 }
69 if (!_f) {
70 return false;
71 }
72
73 seq.clear();
74 streampos spos;
75 while (getline(*ist, line)) {
76 if (line[0] == '>') {
77 break;
78 }
79 seq = seq + line;
80 spos = ist->tellg();
81 }
82 ist->seekg(spos);
83
84 this->ToTLC();
85
86 return (true);
87 }
88
89 void ToTLC(void) {
90 seq3l.clear();
91 for (int i = 0; i < seq.size(); i++) {
92 seq3l.push_back(pair<string, int>(aa3[seq[i]], i + 1));
93 }
94 }
95
96 string GetName(void) { return name; }
97
98 void SetName(string s) { name = s; }
99
100 int length(void) { return seq.size(); }
101
102 void Reset(void) { pos = 0; }
103
104 char GetChar(int p) {
105 if (p > seq.size()) return '\0';
106 return seq[p - 1];
107 }
108
109 string GetTLC(int p) {
110 if (p > seq.size()) return "END";
111 return seq3l[p - 1].first;
112 }
113
114 pair<string, int> Next(void) {
115 if (pos >= seq.size()) return pair<string, int>("END", pos);
116 return seq3l[pos++];
117 }
118
119 void print_seq(ostream* ost) {
120 string x;
121 *ost << ">" << name << endl;
122 int i = 0;
123 for (; i < seq.size() / 60; i++) {
124 x = seq.substr(i * 60, 60);
125 *ost << x << endl;
126 }
127 x = seq.substr(i * 60);
128 *ost << x << endl;
129 }
130 };
131
132 int main(int argc, char* argv[]) {
133 init();
134 SEQ s;
135 ifstream ifs("aligned.fasta");
136 s.ReadFASTA(&ifs);
137 s.print_seq(&cout);
138
139 return 0;
140 }それでは簡単に説明しましょう。10〜36行目ではアミノ酸の1文字表記から3文字表記へ変換するためのデータをmapを利用して設定しています。38行目でアミノ酸配列を扱うクラスを定義が始まります。45〜54行目でコンストラクタ、コピーコンストラクタ、デストラクタの3つのメソッドを定義しています。C++でクラスを作成する際にこの3つのメソッドが定義されていることが大切です。無くても処理系の方で何とかしてくることが多いですが、稀に想定外の動作してバグることがあるので、面倒ですがちゃんと定義しましょう。とは言ってもこの例では大したことは記述していませんけどね。56〜87行目ではFASTAファイルから1つ分のアミノ配列を読み込むメソッドを記述しています。このメソッドではFASTAファイルを読み込んでから1文字表記を3文字表記に変換しています。変換する部分は直ぐ下の89〜94行目の”ToTLC”というメソッドに記述しています。96〜117行目ではSEQクラス中のプロパティにアクセスするための色々なメソッドを記述しています。そして119〜129行目ではFASTA形式で出力する部分を記述しています。
最後の方の132〜140行目にmain関数が記述されています。ただ単にFASTAを読み込んでから、そのままFASTA形式で標準出力するという内容です。出力先をファイルにしたい場合はofstreamで出力先を定義して、そのポインタをprint_seqメソッドの引数に入れてあげればOKです。
ざっくりとした説明でしたが、PDBファイルの残基番号を修正するための準備がこれで整いました。後編ではこのプログラムにさらにOpenBabelを使用する部分を付け足してプログラムを完成させたいと思います。
Category: OpenBabel, プログラミング関連