#はじめに
シリーズ第4回です。

このシリーズは、PDBデータをもとに、SketchUpでタンパク質の立体モデルを作成するシリーズです。
前回までに、ポリペプチドの主鎖と側鎖をedgeとして描画できるようになりました。

今回は、主鎖をfaceを持った立体(円柱が連なったもの)にして色を塗りたいと思います。コメント欄で教えていただいたカッコイイhashの使い方など(@scivola さん、ありがとうございます。)も取り入れていきたいと思います。

タンパク質の立体構造をSketchUpでモデル化その1
タンパク質の立体構造をSketchUpでモデル化その2【アミノ酸の側鎖を描画】
タンパク質の立体構造をSketchUpでモデル化その3【アミノ酸の側鎖も含めて描画】
タンパク質の立体構造をSketchUpでモデル化その4【主鎖を立体に】
タンパク質の立体構造をSketchUpでモデル化その5【玉の描画】
タンパク質の立体構造をSketchUpでモデル化その6【ついでにDNAをモデル化】
タンパク質の立体構造をSketchUpでモデル化その7【最終回】

#方針
エッジの開始端にcircleを描いて、これに対してfollowmeツールを使います。この時、描画済みのedgeを指定します。

絵に描けばこんなです。circleを書く時は、中心座標と半径に加えて、circleの3次元上の向き(つまり3次元ベクトル)を指定します。向きは要するに、円を含む平面の垂線ベクトルを指定するとのことです。ここでは、1つ目のpointから2つ目のpointへと進むベクトルがちょうど垂線ベクトルになります。

Vector3dオブジェクトを得る

Point3dオブジェクトとPoint3dオブジェクトの間で引き算をすると、ベクターが得られます。

p1 = Geom::Point3d.new(1,2,3)
p2 = Geom::Point3d.new(2,2,5)

vector = p1-p2
p vector
# Vector3d(-1, 0, -2)

動作確認

円を描いて、faceにし、さらにフォローミーツールを適用します。

#モデルを取得
model = Sketchup.active_model
entities = model.active_entities

#エッジを追加
p1 = Geom::Point3d.new(1,2,3)
p2 = Geom::Point3d.new(2,2,5)
p3 = Geom::Point3d.new(5,5,5)
edges = entities.add_edges([p1,p2,p3]) #followme用に取っておく

#垂線ベクトル
vector = p1-p2

#circleを追加する。半径は適当に、0.5とする。
r = 0.5
circleEdges = entities.add_circle(p1, vector , r)

#円をfaceに変換。
face = entities.add_face(circleEdges)

#faceにフォローミーツールを適用
face.followme(edges)

動作確認です。

ちゃんと描けました。なお、垂線ベクトルをp2-p1にすると、Faceの裏地が表になった立体になってしまいます。

##色を塗る
これが難航しました。faceが1つでもあれば、それにつながっているface(間接的につながっているものも含めて)を一括して得ることのできるメソッド(all_connected)があるのですが、followmeツールを適用したfaceは、適用時に消されてしまうらしく、followmeで作った立体に対する参照を得られないのです(followmeツールの返り値はbool値)

他に良い方法があるのかもしれませんが、entitiesから一つ一つentityを取り出して、それがFaceであるかをtypenameで調べ、そうだったら先に作った円の中心座標を含むかどうかをclassify_pointで調べ、そうであったらつながっているFaceをall_connectedで取得して、それぞれに色を塗ることにしました。色を塗るには、Faceのmaterialプロパティに値("red"など)をセットします。

for x in entities
    if x.typename == "Face"
        result = x.classify_point(p1)
        if result == Sketchup::Face::PointInside
           connected_entities = x.all_connected
           for y in connected_entities
                if y.typename == "Face"
                   y.material= "red"
                end
            end
        end
    end
end

やっと目処が立ちました。

##関数の定義
Point3dオブジェクトの配列と、円柱の半径、外側の色の3つが与えられれば、立体モデルが作れることになりますので、関数 addStrand()を定義します。

# 立体的なストランドを追加するメソッド、addStrand。
#第1引数は、Point3Dの配列。ストランドの半径と色を指定
#第2引数は、円柱の半径
#第4引数は、色

def addStrand(strand,r,color)
    model = Sketchup.active_model
    entities = model.active_entities
    
    #エッジを加える
    strandEdge = entities.add_edges(strand)

    #circleを追加する。
    vector = strand[0] - strand[1]
    circleEdges = entities.add_circle(strand[0], vector , r)

    #円をfaceに変換。
    face = entities.add_face(circleEdges)
    #faceにフォローミーツールを適用
    face.followme(strandEdge)

    #ストランドに色を付ける。entitiesの中から、フォローミーに使ったFaceの中心座標を含むFaceを探す。
    for x in entities
       if x.typename == "Face"
           result = x.classify_point(strand[0])
           if result == Sketchup::Face::PointInside
               connected_entities = x.all_connected
               for y in connected_entities
                   if y.typename == "Face"
                       y.material = color
                   end
               end
           end
       end
    end
end

#コード全体
前回までのコードを多少変更して、以下のようにしました。色については、あらかじめいつくかの色の名前を配列にして保持しておき、アクセス用のindexを1ずつ増やしながらアクセスすることにしました。循環するように、indexを色の数で割り算したあまりをindexとして使います。

#モデルを取得
model = Sketchup.active_model
entities = model.active_entities

#ファイルを読み込む
lines = IO.readlines("/Users/yoho/Downloads/mmdb_5EPW.pdb")

#関数
# 立体的なストランドを追加するメソッド。ストランドの半径と色を指定
def addStrand(strand,r,color)
    model = Sketchup.active_model
    entities = model.active_entities
    
    #エッジを加える
    strandEdge = entities.add_edges(strand)
    
    #circleを追加する。
    vector = strand[0] - strand[1]
    circleEdges = entities.add_circle(strand[0], vector , r)
    
    #円をfaceに変換。
    face = entities.add_face(circleEdges)
    
    #faceにフォローミーツールを適用
    face.followme(strandEdge)
 
    #ストランドに色を付ける。entitiesの中から、フォローミーに使ったFaceの中心座標を含むFaceを探す。
    for x in entities
       if x.typename == "Face"
           result = x.classify_point(strand[0])
           if result == Sketchup::Face::PointInside
               connected_entities = x.all_connected
               for y in connected_entities
                   if y.typename == "Face"
                       y.material = color
                   end
               end
           end
       end
    end
end

#描画用情報
aminoacidsAndBonds = {
    "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
}
#使用する色のリスト
colors = ["SpringGreen","SlateBlue","MediumVioletRed","MediumBlue","Lime","LightSeaGreen","Indigo",
"Green","Gold","DodgerBlue","DarkOrchid","Crimson","Coral","RoyalBlue"]


#アミノ酸の主鎖データをキープ
strands = []
strandIdentity = ""
for line in lines
    dataInLine = line.split
    if dataInLine[0] == "ATOM"
        if strandIdentity != dataInLine[4]
            points = []
            strands.push(points)
            strandIdentity = dataInLine[4]
        end
        #アミノ酸バックボーンを集める。
        if ["C" ,"CA", "N"].include?(dataInLine[2])
            points.push(Geom::Point3d.new(dataInLine[6].to_f, dataInLine[7].to_f, dataInLine[8].to_f))
        end
    end
end

#主鎖を追加
colorIndex = 0
for strand in strands
    if strand.size >=2
        addStrand(strand,0.1,colors[colorIndex%14])#立体的にしないなら entities.add_edges(strand)
        colorIndex += 1
    end
end

#モノマー(アミノ酸)ごとにデータをキープ
# hashを使用。アミノ酸_strand_アミノ酸番号をkey、配列をvalueにする。配列は行全体を要素とする。
monomersAndLines = Hash.new{ |hash, key| hash[key] = [] }# 存在しないkeyを指定した時に、空の配列が返るカッコイイhash
for line in lines
    dataInLine = line.split
    if dataInLine[0] == "ATOM"
        key = dataInLine[3..5].join("_")
        monomersAndLines[key].push(line)
    end
end
    
# hashのkeyごとに、1つのアミノ酸に関する情報が含まれている。
# アミノ酸の側鎖を描画
monomersAndLines.each{|key,value|
    #最初の行を分割したアレイのindex 3からアミノ酸名を取得
    aa = value[0].split[3]
    #それぞれの行にあるATOMの名前と、座標をHashにしまう。
    atomsAndCoordinates = Hash.new
    for x in value
        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 = aminoacidsAndBonds[aa]
    slicedlist = list.slice(1,list.size-1)# 先頭の["N","CA","C"]は描画しない
    for bonds in slicedlist
        for index in 0..bonds.size - 2
            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
}

#描いた絵
以下のようになりました。当然、回したりできます。

#終わりに
あとは球を追加できるようにすること、ですね。
DNAについても立体構造データがPDBにあるのでそちらもついでに描画する予定です。