タンパク質の立体構造を A-Frame で

はじめに

この記事は @yoho さんの以下の一連の記事に触発されて書いたものです。

タンパク質の立体構造をSketchUpでモデル化 その1 - Qiita
（その 7 まである）

公開されているタンパク質のデータを元に SketchUp というツールで立体構造を描くというもの。
SketchUp は 3 次元モデリングのためのツールで，有償版と無償版があります。プログラミング機能もあるのですが，その言語が Ruby なのだそうです。したがって @yoho さんの記事も Ruby で書かれています。
たいへん興味を惹かれました。

私は，同じデータを使って，やはり Ruby でプログラミングし，しかし別の手段で絵を描かせようと考えました。
別の手段？　んー，よし A-Frame にしよう。

A-Frame は私にはうまく説明できないけど，ウェブブラウザーが持つ高度で高速な 3 次元グラフィックスのためのデータを HTML ベースで記述し，表示するための JavaScript 製ソフトウエア，という理解でいいのだろうか？
もともとバーチャルリアリティー技術として Mozilla が開発したもの。
ナントカゴーゴルを着けて観たら両眼立体視になるのかな。知らんけど。

出来たもの

こんなんできました。（表示されるまで数秒かかります）

↓これはその一コマ

A-Frame の JavaScript を読み込んでいる以外は，自前の JavaScript も CSS も画像も一切使わず，HTML だけでできています。

たいがいのブラウザーで表示できると思います。
アニメーションを指定して，ゆっくり回転するようにしています。

Firefox の場合，カーソルキー ↑ で視点が近づき（分子の中に入ることもできます），カーソルキー ↓ で遠ざかります。

ファイル構成

以下の 6 ファイルです。

• Gemfile
• aframe-1.0.4.min.js（A-Frame の JavaScript）
• amino_acids.yaml（アミノ酸の原子の結合の情報）
• mmdb_5EPW.pdb（原子の位置などの情報）
• protein.rb（HTML 生成プログラム）
• template.html.slim（HTML のテンプレート）

次節からそれぞれについて説明します。

Gemfile

Gemfile はこんだけ：

Gemfile

source "https://rubygems.org"

gem "slim"

HTML を生成するのに Slim テンプレートエンジンを使っています。
Slim なんか使わなくても，ERB でもいいし，文字列リテラルに式展開を埋め込んで HTML を生成してもいいんですけどね。
Slim を使うとメンテナンス性が非常にいいと感じています。

aframe-1.0.4.min.js

A-Frame の JavaScript は下記でダウンロードできます。
https://aframe.io/docs/1.0.0/introduction/installation.html#include-the-js-build
ファイル名はバージョンが分かりやすいように適当に変えました。
なお，上の「出来たもの」の節のリンク先のコードでは，1 ファイルで済ませるためにこの JavaScript は以下の URL を見に行くようにしています。
https://aframe.io/releases/1.0.4/aframe.min.js

amino_acids.yaml

このファイルはアミノ酸の種類ごとに，どの原子とどの原子が繋がっているかを表したもので，@yoho さんの記事からパクって拝借して YAML 形式に直したもの。

amino_acids.yaml

"ALA":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB"]
"CYS":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "SG"]
"ASP":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG"]
  - ["CG", "OD1"]
  - ["CG", "OD2"]
"GLU":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG", "CD"]
  - ["CD", "OE1"]
  - ["CD", "OE2"]
"PHE":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG"]
  - ["CG", "CD1", "CE1", "CZ", "CE2", "CD2", "CG"]
"GLY":
  - ["N", "CA", "C"]
  - ["C", "O"]
"HIS":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG", "ND1", "CE1", "NE2", "CD2", "CG"]
"ILE":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB"]
  - ["CB", "CG1", "CD1"]
  - ["CB", "CG2"]
"LYS":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG", "CD", "CE", "NZ"]
"LEU":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG"]
  - ["CG", "CD1"]
  - ["CG", "CD2"]
"MET":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG", "SD", "CE"]
"ASN":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG"]
  - ["CG", "OD1"]
  - ["CG", "ND2"]
"PRO":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG", "CD", "N"]
"GLN":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG", "CD"]
  - ["CD", "OE1"]
  - ["CD", "NE2"]
"ARG":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG", "CD", "NE", "CZ"]
  - ["CZ", "NH1"]
  - ["CZ", "NH2"]
"SER":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "OG"]
"THR":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB"]
  - ["CB", "OG1"]
  - ["CB", "CG2"]
"VAL":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB"]
  - ["CB", "CG1"]
  - ["CB", "CG2"]
"TRP":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG"]
  - ["CG", "CD2", "CE2", "CZ2", "CH2", "CZ3", "CE3", "CD2"]
  - ["CG", "CD1", "NE1", "CE2"]
"TYR":
  - ["N", "CA", "C"]
  - ["C", "O"]
  - ["CA", "CB", "CG"]
  - ["CG", "CD1", "CE1", "CZ", "CE2", "CD2", "CG"]
  - ["CZ", "OH"]

この，例えば CA というのは炭素原子なのですが，一つのアミノ酸に複数の炭素原子が含まれるので，区別のため，元素記号の C に A とか B とか G とかを付けているのですね（なんで C じゃなくて G かというと，Alpha，Beta，Gamma が由来だから）。

このデータは @yoho さんが作られたものですが，この原子の命名は次節で述べるタンパク質データで使われているものです。

mmdb_5EPW.pdb

あるタンパク質のデータ。
これがどこで公開されているのかよく分からないけど，@yoho さんの記事にまるごと掲載されていたので，パクってコピペして同じファイル名で保存しました。

拡張子が .pdb ですが，中身はただのテキストファイルです。

長いので本記事には載せませんが，冒頭に掲げた @yoho さんの以下の記事の末尾にあるので，そこからコピペしてください。
タンパク質の立体構造をSketchUpでモデル化 その1 - Qiita

protein.rb

これが Ruby スクリプトです。あとの節で解説します。

ruby protein.rb

とやると，result.html が生成します。それをブラウザーで開くだけ。

template.html.slim

HTML のテンプレートです。Slim というテンプレート言語で書いています。
以下のようにとても簡素です。

template.html.slim

doctype html
html lang="ja"
  meta charset="UTF-8"
  title protein
  script src="aframe-1.0.4.min.js"
  body
    main
      a-scene
        a-sky color="#334"
        a-camera position="0 0 80" near="3" fov="60"
        a-entity animation=animation
          - geometries[:spheres].each do |sphere|
            a-sphere *sphere
          - geometries[:cylinders].each do |cylinder|
            a-cylinder *cylinder

Slim を知らなくても，インデントが階層を表していると知れば，上半分はなんとなく分かるのではないでしょうか。

少しだけ説明します。

a-entity animation=animation

は a-entity 要素（つまり <a-entity> 云々 </a-entity>）を生成するのですが，

animation=animation

の左辺の animation は属性名です。右辺の animation はローカル変数を参照していて，この値は Ruby スクリプト（protein.rb）から渡されます。

また，

- geometries[:spheres].each do |sphere|
  a-sphere *sphere

の - のあとは Ruby コードです。
geometries も Ruby スクリプトから渡されるローカル変数です。

a-sphere *sphere

は a-sphere 要素（球を描く要素）を生成しますが，*sphere の sphere は a-sphere 要素に与える属性を保持しているハッシュです。
こんなふうに，ハッシュに * を付けて渡すとそれが要素の属性の指定になります。
また，

- geometries[:cylinders].each do |cylinder|
  a-cylinder *cylinder

も同様にして，円柱を描く a-cylinder 要素を生成しています。

A-Frame の座標系

次の節の説明で必要になる，A-Frame の座標系について簡単に説明します。

A-Frame は右手系の座標系を採用しています。
まず右手の親指，人差し指，中指を互いに垂直になるように曲げます。曲げ方は二通り考えられますが，どう考えても親指・人差し指をピンと立て，中指を内側に折ったほうが自然ですよね。
この状態で，親指，人差し指，中指をそれぞれ X 軸，Y 軸，Z 軸とします。指の向かう方向が正の方向。
これが右手系。

そして，A-Frame は画面に向かって右が X 軸，上が Y 軸，手前が Z 軸です。
なので，さきほどの右手でいうと，中指が自分の顔を指しているわけですね。

ところで，さきほどのテンプレートでは，カメラ（a-camera 要素）の位置が

position="0 0 80"

となっていました。
Z 座標が正になっているのは，手前が正だからです。
ちなみに単位は m（メートル）です。VR なので，ピクセルじゃなくて現実世界の単位が使われているのですね。
原点から 80 m 手前にカメラがあるわけです。

タンパク質のデータで，原子の位置は Å（オングストローム）単位で表されています。
このデータの数値を，スケールしないでそのまんま A-Frame に持ってきているので，例えば 1 Å 離れた原子は，1 m 離れているように描画されるわけですね。

この記事で使ったタンパク質では，すべての原子は半径 40 Å の球の中に収まっています。VR データとしては半径 40 m ですね。そういう原子たちを 80 m 離れたところから見ている，ということになります。

スクリプト

では Ruby スクリプトです。

protein.rb

require "matrix"
require "yaml"
require "bundler"
Bundler.require

AMINO_ACIDS = YAML.load_file("amino_acids.yaml")

# 描画に関わる元素の情報
# （半径はテキトー）
ELEMENTS = {
  "C" => {radius: 0.70, color: "#909090"},
  "N" => {radius: 0.65, color: "#3050F8"},
  "O" => {radius: 0.60, color: "#FF0D0D"},
  "S" => {radius: 1.00, color: "#FFFF30"},
}

class Atom
  # データファイルで各原子に付けられた名前（C，CA，N など）
  attr_reader :name

  # 元素名（C，N，O など）
  attr_reader :element

  # アミノ酸識別子
  attr_reader :amino_acid_id

  # アミノ酸名（LYS などの略称）
  attr_reader :amino_acid_name

  # ポリマー ID（？）
  attr_reader :polymer_id

  # 原子の位置
  attr_accessor :position

  def initialize(name, element:, polymer_id:, amino_acid_name:, amino_acid_id:, position:)
    @name = name
    @element = element
    @polymer_id = polymer_id
    @amino_acid_name = amino_acid_name
    @amino_acid_id = amino_acid_id
    @position = Vector[*position]
  end
end

atoms = []

# データファイルをパースして Atom オブジェクトの配列を得る
IO.foreach("mmdb_5EPW.pdb") do |line|
  type, _number, name, aa_name, p_id, aa_number, x, y, z, *, element = line.split
  next unless type == "ATOM"

  position = [x, y, z].map(&:to_f)
  # アミノ酸名，ポリマー ID，アミノ酸番号を使ってアミノ酸識別子とする
  amino_acid_id = [aa_name, p_id, aa_number].join(" ")

  atom = Atom.new name,
    element: element,
    polymer_id: p_id,
    amino_acid_name: aa_name,
    amino_acid_id: amino_acid_id,
    position: position
  atoms << atom
end

# 重心が原点になるよう全原子の位置をずらす
center_of_gravity = atoms.map(&:position).inject(:+) / atoms.length
atoms.each do |atom|
  atom.position -= center_of_gravity
end

# 結合
bonds = []

atoms.group_by(&:amino_acid_id).each do |amino_acid_id, a|
  amino_acid_name = a.first.amino_acid_name
  AMINO_ACIDS[amino_acid_name].each do |b|
    b.each_cons(2) do |x1, x2|
      a1 = a.find{ _1.name == x1 }
      puts "Not found: #{x1} in #{amino_acid_id}" unless a1
      a2 = a.find{ _1.name == x2 }
      puts "Not found: #{x2} in #{amino_acid_id}" unless a2

      if a1 && a2
        bonds << [a1, a2, amino_acid_name]
      end
    end
  end
end

# 原子数，結合数，原点から最も遠い原子までの距離を表示
puts "%d atoms." % atoms.count
puts "%d bonds." % bonds.count
puts "Max distance from the origin: %0.3f" % atoms.map(&:position).map(&:norm).max

# 二つの座標を結ぶ円柱を描画するための基本パラメーターを作る
def cylinder_param(v1, v2)
  v21 = v2 - v1
  n = v21 / v21.norm # 単位方向ベクトル
  x, y, z = n.to_a
  if x ==0 && z == 0
    rot_y = 0
  else
    rot_y = -Math.atan2(z, x) / Math::PI * 180
  end
  {
    # 円柱の中心座標
    center: (v1 + v2) * 0.5,
    # 円柱の高さ
    length: v21.norm,
    # Y 軸まわりの回転角（度）
    rot_y: rot_y,
    # Z 軸まわりの回転角（度）
    rot_z: -Math.acos(y) / Math::PI * 180
  }
end

# 数値の組や Vector から丸めて空白で区切った文字列を作る
def rounded_numericals(ary, digit = 5)
  ary.to_a.map{ |v| v.round(digit) }.join(" ")
end

# A-FRAME の図形情報
geometries = {spheres: [], cylinders: []}

# 原子の情報から球の描画情報を作る
atoms.each do |atom|
  elem_prop = ELEMENTS[atom.element]
  geometries[:spheres] << {
    position: rounded_numericals(atom.position),
    radius: (elem_prop[:radius] * 0.7).round(3),
    color: elem_prop[:color]
  }
end

# 結合の情報から円柱の描画情報を作る
bonds.each do |atom1, atom2, amino_acid_name|
  cy = cylinder_param(atom1.position, atom2.position)
  geometries[:cylinders] << {
    position: rounded_numericals(cy[:center]),
    color: "#DDD",
    radius: 0.15,
    height: cy[:length].round(3),
    rotation: rounded_numericals([0, cy[:rot_y], cy[:rot_z]], 3)
  }
end

# アニメーション（Y 軸周りのゆっくり回転）の設定
animation = <<~EOT.gsub(/\n/, ";")
  property: rotation
  from: 0 0 0
  to: 0 360 0
  dur: 20000
  easing: linear
  loop: true
EOT

template = Slim::Template.new("template.html.slim")

IO.write "result.html", template.render(nil, geometries: geometries, animation: animation)

あまり詳しく説明しませんが，質問歓迎です。

元素の半径と色

ELEMENTS 定数で元素の原子半径と描画色を定義しています。

原子にはそもそも半径などというものは無いので，何らかの意味で半径的なものを考えるのですが，目的に応じてさまざまな半径の定義があるそうです。
今回使った原子半径の値はどこから持ってきたのか忘れました。
しかも，その値を直接は使わず，見た目でいい感じになるように 0.7 倍しちゃってます。この倍率も目的に応じていろいろ変えてみるとよさそうです。

次に元素の色ですが，Wikipedia CPK配色 に掲載されていた「立体分子模型画像作成ソフトJmolの配色」というものを使いました。
これは，よく見かけるプラスチックの分子模型に使われている CPK 配色（もしくはその仲間）に近くて，多くの元素に拡張したものになっているようです。
CPK 配色と大きく違うのは，炭素が黒でなく灰色であることです。黒だと背景が黒だったら見えませんしね。

原子の位置

Atom クラスは，個々の原子の情報を表すクラスです。

3 次元の位置や変位は Vector オブジェクトで表します。

重心が原点にくるよう，まず全原子の位置の平均

atoms.map(&:position).inject(:+) / atoms.length

として重心を求め，これを差っ引くことで位置を修正しています。

原子間の結合

原子間の結合は A-Frame の a-cylinder オブジェクトで表すのですが，円柱の中心座標を position="2.34 -31.93 1.34" のように与えます。
二つの原子の座標を表す Vector オブジェクトを v1，v2 とすると，(v1 + v2) * 0.5 でその中間の座標が得られますね。

円柱の高さは (v1 - v2).norm で得られます。Vector#norm はベクトルのノルム（大きさ）です。

問題は円柱の向きなんですよね。
何も指定しないと，Y 軸（垂直上向き）を軸とした円柱になります。
これを回転させます。rotation 属性で rotation="23 -90 45" のように与えるのですが，この三つの数値は，それぞれ，X 軸，Y 軸，Z 軸まわりの回転角（度）を表しています。

ここまではさまざまな解説記事にも書いてあるのですが，問題は回転の順序と向きです。rotation の説明をするなら必ず書いてほしい！
X 軸を回してから Y 軸を回すのと，Y 軸を回してから X 軸を回すのでは結果が違うので，順序が大事なんですよ。
試したところ，どうやら，Z 軸 → Y 軸 → X 軸の順に回すようです。

これでようやく回転の向きの話ができます。
X 軸まわりの回転は，X 軸の正の向きに向かってネジを進めるときの回転の向きを正とします。言い換えると，X 軸の正の向きに親指を突き立て，他の四本の指を握ったとき，この四本の指の向きが正の回転の向きです。

ここまで分かれば円柱の回転について考えることができます。
Y 軸方向を向いている円柱を，適当に回転させて v1 から v2 へ向かうベクトルと同じ方向にしてやらなければなりませんが，これは cylinder_param メソッドで計算しています。
説明は割愛しますが，[0, 1, 0] というベクトルを Z 軸周りに回転させ，それを Y 軸周りに回転させます。このときのそれぞれの回転角を，アークタンジェント（正接関数の逆関数），アークコサイン（余弦関数の逆関数）を使って求めています。
図にすればそれほど難しいものではないのですが，描く元気はありませんでした。

発展

今回作ったのはただ立体模型がくるくる回るだけ。ブラウザーの機能を使ってもカメラの位置を変えることくらいしかできません。

おそらく簡単な JavaScript を追加すれば，好きな向きに回転させたりできるのでしょう。
また，結合を非表示にして原子だけにしたり，特定の種類のアミノ酸だけを目立たせたり，といったことも考えられます。

筆者は生化学とかまったく分からないので，ただ「絵がデター！　ひゃほー！」としか思っていません。
どういう表示をさせたらどういう科学的意義があるのかも見当がつきません。

そういえば学生の頃，生物学の講義で一枚のプリントが配られ，「これが 20 種のアミノ酸だ。構造を全部覚えてくるように」と言われました。
こりゃ大変だと必死で（半分以上？）覚えましたが，試験には一問も出ず，たしか講義でもこのアミノ酸の構造を用いた話は一切無かったように思います。なんだったんだろう？
私にとってアミノ酸とは，この一件が全てです。
しかし，この講義は戦時中の体験談などの雑談が面白くて 30 年経った今でも覚えています。