#はじめに
初回は、タンパク質の主鎖をSketchUpで立体モデル化しました。第二回では、アミノ酸の側鎖を描画するためのデータを整理しました。
今回は、主鎖と側鎖を合わせて立体モデル化します。
rubyは素人ですので、おかしな書き方があるかもしれません。自分用の忘備録という位置付けですが、コメントなど頂けますと幸いです。
@scivolaさんのコメントを参考に、コードを修正したところがあります。ありがとうございます。
タンパク質の立体構造をSketchUpでモデル化 その1
タンパク質の立体構造をSketchUpでモデル化 その2【アミノ酸の側鎖を描画】
タンパク質の立体構造をSketchUpでモデル化 その3【アミノ酸の側鎖も含めて描画】
タンパク質の立体構造をSketchUpでモデル化 その4【主鎖を立体に】
タンパク質の立体構造をSketchUpでモデル化 その5【玉の描画】
タンパク質の立体構造をSketchUpでモデル化 その6【ついでにDNAをモデル化】
タンパク質の立体構造をSketchUpでモデル化 その7【最終回】
#PDBデータファイルからアミノ酸ごとにデータを集める
hashを使って、タンパク質中の各アミノ酸ごとに、関連するデータ行をリストに入れます。
hashのkeyは、Phe_A_329、などとします。Aはストランドを区別するための識別子です。valueは、行のデータそのものです。
新しいATOMが先頭にある行データを調べるたびに、上記のkeyを生成し、これがhashに登録されているか調べます。登録されていなければ、行データをしまった配列を作って、これをkeyのvalueとして新規登録します。登録されていれば、keyに対応する配列に、行データを追加します。
aminoAcidAndLines = Hash.new
for line in lines do
dataInLine = line.split
if dataInLine[0] == "ATOM" then
key = dataInLine[3] + "_" + dataInLine[4] + "_" + dataInLine[5]
if aminoAcidAndLines.has_key?(key) then
aminoAcidAndLines[key].push(line)
else
aminoAcidAndLines[key] = [line]
end
end
end
#側鎖を描画する
上で作成したhashについて、登録されているデータを取り出して描画します。
まずは、何のアミノ酸についてのデータかを調べます。次に、ATOMデータから、原子の名前(NとかCAとかCBとか)を調べ、対応する座標(Point3dオブジェクト)を、紐づけておきます。
処理の後半では、第二回で準備したデータを参照して、アミノ酸ごとにどの原子とどの原子の間にedgeを作るか調べます。アミノ酸ごとに、例えば「アラニンだったらCAの原子からCBの原子に、edgeを書く」というようなデータが得られますが、実際のATOMデータには、必ずしもすべての原子の位置データが記されている話ではありません。
ですので、間に線を引くべき2つの原子それぞれがデータとしてあるかどうかを調べ、両方ある場合のみedgeを追加することにしました。結合情報がない場合は、その旨を出力することにしました。
なお、アミノ酸データには、側鎖の情報だけでなく、主鎖の情報も含まれていました。主鎖の部分を2度描いてしまうのは嫌なので、listをsliceして2番目のデータ(["C","O"])から描画することにしました。主鎖と重なる部分は、そもそもaminoacidAndBondsに含めなければよかったかな.....
# hashのkeyごとに、1つのアミノ酸に関する情報が含まれている。
aminoAcidAndLines.each{|key,value|
#最初の行を分割したアレイのindex 3からアミノ酸名を取得
aa = value[0].split[3]
#それぞれの行にあるATOMの名前と、座標をHashにしまう。
atomsAndCoordinates = Hash.new
for x in value do
data = x.split
atomsAndCoordinates[data[2]] = Geom::Point3d.new(data[6].to_f, data[7].to_f, data[8].to_f)
end
#アミノ酸ごとに、描画すべき結合を、edgeとして追加する。
#アミノ酸の原子情報が欠けている場合があるので、bondsから2つずつ取り出して、対応する座標データが両方得られる場合に限って、edgeを追加する。
list = aminoacidAndBonds[aa]
slicedlist = list.slice(1,list.size-1)# 先頭の["N","CA","C"]は描画しない
for bonds in slicedlist
for index in 0..bonds.size - 2 do
atomA = atomsAndCoordinates[bonds[index]]
atomB = atomsAndCoordinates[bonds[index+1]]
if atomA != nil and atomB != nil then
entities.add_edges([atomA,atomB])
else
puts "結合情報がありません" + key +" "+ aa + " "+ bonds[index] + "-" + bonds[index + 1]
end
end
end
}
#描画してみる
コードは全部で以下のようになりました。なおGeomモジュールをincludeする手もある、とコメントをいただきましたが、どこのPoint3dか常に明示した方がわかりやすいかなと思い、逐一指定することしました。
#モデルを取得
model = Sketchup.active_model
entities = model.active_entities
#ファイルを読み込む
lines = IO.readlines("/Users/yoho/Downloads/mmdb_5EPW.pdb")
#主鎖データをキープ
strands = []
strandIdentity = ""
for line in lines do
dataInLine = line.split
if dataInLine[0] == "ATOM" then
if strandIdentity != dataInLine[4] then
points = []
strands.push(points)
strandIdentity = dataInLine[4]
end
if ["C" ,"CA", "N"].include?(dataInLine[2]) then
points.push(Geom::Point3d.new(dataInLine[6].to_f, dataInLine[7].to_f, dataInLine[8].to_f))
end
end
end
#主鎖を描画edgeとして、追加
for strand in strands
entities.add_edges(strand)
end
#側鎖データをキープ
# hashを使用。strand_アミノ酸番号をkey、配列をvalueにする。配列は行全体を要素とする。
aminoAcidAndLines = Hash.new
for line in lines do
dataInLine = line.split
if dataInLine[0] == "ATOM" then
key = dataInLine[3] + "_" + dataInLine[4] + "_" + dataInLine[5]
if aminoAcidAndLines.has_key?(key) then
aminoAcidAndLines[key].push(line)
else
aminoAcidAndLines[key] = [line]
end
end
end
aminoacidAndBonds = {
"ALA" => [["N","CA","C"],["C","O"],["CA","CB"]],#A
"CYS" => [["N","CA","C"],["C","O"],["CA","CB","SG"]],#C
"ASP" => [["N","CA","C"],["C","O"],["CA","CB","CG"],["CG","OD1"],["CG","OD2"]],#D
"GLU" => [["N","CA","C"],["C","O"],["CA","CB","CG","CD"],["CD","OE1"],["CD","OE2"]],#E
"PHE" => [["N","CA","C"],["C","O"],["CA","CB","CG"],["CG","CD1","CE1","CZ","CE2","CD2","CG"]],#F
"GLY" => [["N","CA","C"],["C","O"] ],#G
"HIS" => [["N","CA","C"],["C","O"],["CA","CB","CG","ND1","CE1","NE2","CD2","CG"]],#H
"ILE" => [["N","CA","C"],["C","O"],["CA","CB"],["CB","CG1","CD1"],["CB","CG2"]],#I
"LYS" => [["N","CA","C"],["C","O"],["CA","CB","CG","CD","CE","NZ"]],#K
"LEU" => [["N","CA","C"],["C","O"],["CA","CB","CG"],["CG","CD1"],["CG","CD2"]],#L
"MET" => [["N","CA","C"],["C","O"],["CA","CB","CG","SD","CE"]],#M
"ASN" => [["N","CA","C"],["C","O"],["CA","CB","CG"],["CG","OD1"],["CG","ND2"]],#N
"PRO" => [["N","CA","C"],["C","O"],["CA","CB","CG","CD","N"]],#P
"GLN" => [["N","CA","C"],["C","O"],["CA","CB","CG","CD"],["CD","OE1"],["CD","NE2"]],#Q
"ARG" => [["N","CA","C"],["C","O"],["CA","CB","CG","CD","NE","CZ"],["CZ","NH1"],["CZ","NH2"]],#R
"SER" => [["N","CA","C"],["C","O"],["CA","CB","OG"]],#S
"THR" => [["N","CA","C"],["C","O"],["CA","CB"],["CB","OG1"],["CB","CG2"]],#T
"VAL" => [["N","CA","C"],["C","O"],["CA","CB"],["CB","CG1"],["CB","CG2"]],#V
"TRP" => [["N","CA","C"],["C","O"],["CA","CB","CG"],["CG","CD2","CE2","CZ2","CH2","CZ3","CE3","CD2"],["CG","CD1","NE1","CE2"]],#W
"TYR" => [["N","CA","C"],["C","O"],["CA","CB","CG"],["CG","CD1","CE1","CZ","CE2","CD2","CG"],["CZ","OH"]]#Y
}
# hashのkeyごとに、1つのアミノ酸に関する情報が含まれている。
aminoAcidAndLines.each{|key,value|
#最初の行を分割したアレイのindex 3からアミノ酸名を取得
aa = value[0].split[3]
#それぞれの行にあるATOMの名前と、座標をHashにしまう。
atomsAndCoordinates = Hash.new
for x in value do
data = x.split
atomsAndCoordinates[data[2]] = Geom::Point3d.new(data[6].to_f, data[7].to_f, data[8].to_f)
end
#アミノ酸ごとに、描画すべき結合を、edgeとして追加する。
#アミノ酸の原子情報が欠けている場合があるので、bondsから2つずつ取り出して、対応する座標データが両方得られる場合に限って、edgeを追加する。
list = aminoacidAndBonds[aa]
slicedlist = list.slice(1,list.size-1)# 先頭の["N","CA","C"]は描画しない
for bonds in slicedlist
for index in 0..bonds.size - 2 do
atomA = atomsAndCoordinates[bonds[index]]
atomB = atomsAndCoordinates[bonds[index+1]]
if atomA != nil and atomB != nil then
entities.add_edges([atomA,atomB])
else
puts "結合情報がありません" + key +" "+ aa + " "+ bonds[index] + "-" + bonds[index + 1]
end
end
end
}
#描画結果
こんなふうになりました。いわゆる針金モデルですね。
#今後
- 原子それぞれを球として描画
- 主鎖を立体構造に
といったことが課題でしょうか。タンパク質の表面を描画したいと思っていましたが、これはちょっと先になりそうです。
なかなかrubyの文法が頭に入らず、コードにするのに時間がかかります。頑張って覚えてしまわないとダメですね。