#はじめに
初回は、タンパク質の主鎖を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の文法が頭に入らず、コードにするのに時間がかかります。頑張って覚えてしまわないとダメですね。