Documents

QScriptでファイル入出力を行なう

QScriptを使用するとPerlライクな 簡単なテキストファイルの入出力が行なえます． このページではこの機能を分子構造の表示に応用したケースを紹介します．

↑

DSSPで計算した二次構造を読み込む

現時点でCueMol自体には二次構造を計算する機能がありません． ですので， 二次構造情報が入っていないPDBファイルで リボンモデルを表示する場合， DSSP等で計算した二次構造情報を読込んでやる必要があります．

Documents/QScriptのチュートリアル/Step3ではprocheckの出力から 手動で設定する例を紹介しましたが， ここではDSSPの出力ファイルを自動的に解析して割り当てる例を紹介します．

↑

ファイル

以下のファイルをダウンロードし，同じディレクトリに置いてください． そして，CueMolを起動してメニュー「File」→「Execute QScript...」から fileio1.qsを開いてください．

fileio1.png
↑

DSSPファイル

DSSPの出力ファイルは，はじめの24行はヘッダーで，25行目から残基ごとの情報になります．

  #  RESIDUE AA STRUCTURE BP1 BP2  ACC     N-H-->O    O-->H-N    N-H-->O    O-->H-N    TCO  KAPPA ALPHA  PHI   PSI    X-CA   Y-CA   Z-CA 
    1    1   L              0   0  105      0, 0.0     2,-0.2     0, 0.0   178,-0.0   0.000 360.0 360.0 360.0 148.0    1.9   27.4   55.8
    2    2   T     >  -     0   0   75      1,-0.1     4,-2.4   177,-0.0     5,-0.3  -0.500 360.0-113.9 -86.9 164.2    0.6   27.5   52.2
    3    3   A  H  > S+     0   0   69      1,-0.2     4,-2.2     2,-0.2     5,-0.2   0.923 118.8  49.9 -62.4 -41.6   -2.4   25.5   51.3
    4    4   N  H  > S+     0   0   89      1,-0.2     4,-2.7     2,-0.2    -1,-0.2   0.934 109.7  50.4 -63.2 -45.3   -0.2   23.3   49.1
    5    5   Q  H  > S+     0   0   23      1,-0.2     4,-2.2     2,-0.2    -1,-0.2   0.892 111.2  45.7 -67.0 -39.8    2.3   22.7   51.8

以上は24行目から抜粋したものです． RESIDUEの部分（６〜１０文字目）に残基番号，STRUCTUREの部分（２５文字目）に二次構造の情報が入っています． これらを取り出して処理すればよいわけです．

↑

スクリプトの解説

# スクリプトがあるディレクトリーのパス名を取得
$pwd = sys.getScriptPath();

# PDBファイルを開き，
#  "ecr"という名前の分子オブジェクトを作成
#  分子オブジェクトを$molに代入する．
$mol = readPDB("$pwd/ECR_MODEL.pdb","ecr");

# ローカル変数は以下の{}(スコープ)内のみで有効
{
    # ローカル変数の宣言．
    local $line;
    # テキストファイルdssp.txtを開いて入力ストリームを作成．
    # そして$finに代入している．
    local $fin = fistream("$pwd/dssp.txt");
    # 最初の24行はスキップ．（ログとして表示）
    foreach $i (1..24) {
        $line = $fin.readline();
        $line.chomp();
        sys.println($line);
    }

    # 25行目から行末までは，残基ごとの情報が入っている．
    # 各行を読み込み処理．
    while ($fin.ready()) {
        # 入力ストリーム$finから一行読込む
        $line = $fin.readline();

        # 行の6文字目から5文字は残基番号
        # 　substrは0始まりで文字の位置を数える（Perlと同じ）
        # 　integer()で文字列から整数に変換している
        $resid = integer($line.substr(5, 5));

        # 行の17文字目から1文字は二次構造→$codeに格納
        $code = $line.substr(16, 1);

        # $codeがHの場合はhelix, Eの場合はsheet, 
        # それ以外はcoilに変換
        if ($code == "H") {
            $code = "helix";
        }
        else if ($code == "E") {
            $code = "sheet";
        }
        else {
            $code = "coil";
        }

        # 現在の残基（$resid）を選択
        $mol.selectResid("_", $resid, $resid);
        # 現在の残基のsecondaryプロパティーを設定
        $mol.setProp("secondary", $code);
    }
}
# 以上のスコープから出るとローカル変数$finは破棄される
# →ファイルはクローズされる．

$mol.select(se/resn BLK/);
... (省略) ...
$r_p.setProp("sheet_tail.gamma", 4.0);

# 残基プロパティーsecondaryがhelixの部分を選択
$mol.select(se/rprop secondary=helix/);
# 赤色に着色
molvis.paint($r_p, color(1,0,0));
# 残基プロパティーsecondaryがsheetの部分を選択
$mol.select(se/rprop secondary=sheet/);
# 青色に着色
molvis.paint($r_p, color(0,0,1));

... (以下省略) ...
↑

Verify3Dのスコアで着色

実は，上記で使用しているECR_MODEL.pdbは， ヒトのグルココルチコイドレセプターから 昆虫のエクジソンレセプターを比較モデリングした モデル構造です． 例えばモデル構造の評価にVerify3Dを使用した場合， そのスコアによって残基を着色すれば モデルの良し悪しが分かりやすくなります．

↑

ファイル

以下のファイルをダウンロードし，同じディレクトリに置いてください． そして，CueMolを起動してメニュー「File」→「Execute QScript...」から

  • PDBファイル: fileECR_MODEL.pdb
  • Verify3Dのスコア：fileverify3d.txt
    これは，Verify3DのウェブサイトにECR_MODEL.pdbをサブミットすると出てくる結果ページの，「Display the raw average data」ボタンをクリックすると出てくるページからコピペしたものです．
  • QScriptファイル: filefileio2.qs
fileio2.png
↑

スクリプトの解説

主なポイント

  • 複雑なフォーマットのテキストファイルを解析するには，正規表現を使えば簡単に実現できます．（このケースではsubstrでも出来ますが，例として正規表現を使っています．）QScriptではPerlとほぼ同じ正規表現が使用できます．
  • テキストから読み込んだ数値を原子に割り当てるには，温度因子あるいはoccupancyに代入しておく．
  • simple, traceレンダラーの場合はこの例のように簡単に温度因子(or occupancy)でグラジエント着色できる．
  • tubeやribbonレンダラーの場合は次のセクションで説明． 
    $pwd = sys.getScriptPath();
$mol = readPDB("$pwd/ECR_MODEL.pdb","ecr");

# スコアの最大値と最小値
$smin=1.0e10;
$smax=-1.0e10;
{
    local $line;
    # 正規表現オブジェクトを作成 →$reに代入
    #  何回も使用する場合は，このようにループ等の外で
    #  予め正規表現オブジェクトを作っておいた方が高速．
    local $re = re/avg\s+\w\s+(\d+)\s+([\d\.]+)/;
    # 入力ストリームを開く →$finに代入
    local $fin = fistream("$pwd/verify3d.txt");
    # read verify3d score
    while ($fin.ready()) {
        $line = $fin.readline();
        # 正規表現に入力行$lineがマッチしなかった場合は
        # エラーを出して継続
        if (!$re.match($line)) {
            sys.println("Invalid line: $line");
            continue;
        }
        # グルーピングを使用してマッチした部分を取り出す
        #   at(1)で，1番目の括弧中にマッチした部分
        #   (Perlで言うところの$1)を取り出せる
        $resid = integer($re.at(1));
        #   at(2)で，2番目の括弧中にマッチした部分
        #   (Perlで言うところの$2)を取り出せる
        $score = real($re.at(2));

        # 最大値と最小値の更新
        if ($score<=$smin) $smin = $score;
        if ($score>=$smax) $smax = $score;

        # 現在の残基を選択
        $mol.selectResid("_", $resid, $resid);
        # 選択された残基の温度因子を$scoreに書き換える
        $mol.forEachAtom() ($atom) {
            $atom.setBfac($score);
        };
    }
    # sys.printfで最大最小スコアを表示．
    #   printfはC言語のものと大体同じ仕様だが，
    #   QScriptでは関数に可変引数を渡せないので，
    #   このようにリストにして渡してやる必要がある．
    sys.printf("max score: %f, min score: %f\n", [$smax, $smin]);
}
... (途中省略) ...
$mol.select(se/resi 1:239/);
$r_p = $mol.createRend("p", "trace");
$r_p.setProp("linew", 5.0);
# 温度因子による着色モード
$r_p.setProp("coloring.colormode", 1);
# 低いほうは$sminが赤色, 高いほうは$smaxが白色になるように着色
#  Verify3Dのスコアは高いほうが良いので，
#  良い構造の部分が白に，悪い構造の部分が赤になる．
$r_p.setProp("coloring.lowpar", $smin);
$r_p.setProp("coloring.highpar", $smax);
$r_p.setProp("coloring.lowcol", color(1.0, 0.0, 0.0));
$r_p.setProp("coloring.highcol",  color(1.0, 1.0, 1.0));
... (以下省略) ...
↑

Tube, ribbon rendererの場合

tubeやribbonレンダラーの場合は「スクリプトによる着色」を使用して着色する必要があります． この場合は別に温度因子やoccupancyに値を代入しなくてもよく， 独自の配列（実際には連想配列；QScriptではdict型）を作ってそれを参照して着色するスクリプトを設定します．

なぜ配列に単純なリストlist型を使わずに連想配列dict型を使うかというと， list型はリンクトリストで実装されているためで， インデックスでアクセスするとO(N)の時間がかかり遅いからです． 連想配列は木構造で実装されているため， インデックスでアクセスしてもO(logN)の時間でアクセスできます． ただし，そのためには数値を文字列に変換する必要があります．

以下がスクリプトファイルになります．PDBやverify3d.txtは上記と同じものを使用してください．

fileio3.png
Last-modified: Sun, 07 Aug 2005 02:25:53 JST (7048d)