タンパク質の立体構造をSketchUpでモデル化 その7【最終回】

#はじめに
PDB形式のタンパク質立体構造データをもとに、SketchUpを使ってタンパク質の立体モデルを作成するシリーズの最終回です。

[タンパク質の立体構造をSketchUpでモデル化 その1](https://qiita.com/yoho/items/1e15bbd4e3a7c4c7ff69)
[タンパク質の立体構造をSketchUpでモデル化 その2【アミノ酸の側鎖を描画】](https://qiita.com/yoho/items/05037dd8f128c77152d7)
[タンパク質の立体構造をSketchUpでモデル化 その3【アミノ酸の側鎖も含めて描画】](https://qiita.com/yoho/items/3c8e5eb48d68324cec60)
[タンパク質の立体構造をSketchUpでモデル化 その4【主鎖を立体に】](https://qiita.com/yoho/items/d5626464fd39903e0d51)
[タンパク質の立体構造をSketchUpでモデル化 その5【玉の描画】](https://qiita.com/yoho/items/c85a37e9937ddbbaab9c)
[タンパク質の立体構造をSketchUpでモデル化 その6【ついでにDNAをモデル化】](https://qiita.com/yoho/items/0c38e1bca3b1ce17fbb5)
[タンパク質の立体構造をSketchUpでモデル化 その7【最終回】](https://qiita.com/yoho/items/a8204082a216e35462d7)


前回までに、タンパク質の主鎖を細い円柱で、側鎖を細い線(edge)で描くこと、DNAについても、リン酸-デオキシリボースのバックボーンを細い円柱で、それ以外をedgeで描くことができていました。ところがコードがタンパク質用とDNA用に別れてしまい、冗長になってしまいました。

また、メインルーチンが長く、いかにも保守しにくそうであること、各行のデータに何度もsplit処理をするなど、無駄が多くなっていました。

今回は、これらをまとめてすっきりさせ、保守しやすくしたいと思います。なおRNAを含む立体構造データはほとんどないようでしたので、RNAの描画については考えないことにしました。RNAは分解しやすいので立体構造解析が難しいのでしょうね。

立体データはMMDB ID:71816のもので、逆転写酵素と2本鎖DNAの複合体です。

#コード
以下のメソッドを用意しました。
1. 主鎖を描画するメソッド
2. ポリマーの構成単位(つまりモノマー)を描画するメソッド
3. 球を描画するメソッド(ここでは不使用)
4. データタイプ(DNA or Protein)の判定メソッド
また、DNAとProteinの原子の結合状態データを1つにまとめ、定数である(変更しない)変数を、スクリプト上部に記述することにしました。


```rb
#タンパク質とDNAのPDBデータ形式をSketchUpの立体構造モデルに変換するRuby Scriptです。ご利用、改変は構いませんが、当スクリプトのこの行はそのままにするようお願いします。作者 @yoho at Qiita

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

#読み込むファイルへのパス
inputFilePath = "/Users/yoho/Downloads/mmdb_1D0E.pdb"

#描画する主鎖の原子
typesAndMainstrand = {"DNA"=>["P","O5*","C5*","C4*","C3*","O3*","P"],"Protein"=>["C" ,"CA", "N"]}

#描画用結合情報
$monomersAndBonds = {
    # DNAはindex 0,1が主鎖。塩基部分は、index 6以降
    "A" => [["P","O5*","C5*","C4*","C3*"],["C3*","O3*"],["P","O1P"],["P","O2P"],["C3*","C2*","C1*"],["C4*","O4*","C1*"],["C1*","N9","C8","N7","C5","C6","N1","C2","N3","C4","N9"],["C4","C5"],["C6","N6"]],#A
    "G" => [["P","O5*","C5*","C4*","C3*"],["C3*","O3*"],["P","O1P"],["P","O2P"],["C3*","C2*","C1*"],["C4*","O4*","C1*"],["C1*","N9","C8","N7","C5","C6","N1","C2","N3","C4","N9"],["C4","C5"],["C6","O6"],["C2","N2"]],#G
    "C" => [["P","O5*","C5*","C4*","C3*"],["C3*","O3*"],["P","O1P"],["P","O2P"],["C3*","C2*","C1*"],["C4*","O4*","C1*"],["C1*","N1","C2","N3","C4","C5","C6","N1"],["C2","O2"],["C4","N4"]],#C
    "T" => [["P","O5*","C5*","C4*","C3*"],["C3*","O3*"],["P","O1P"],["P","O2P"],["C3*","C2*","C1*"],["C4*","O4*","C1*"],["C1*","N1","C2","N3","C4","C5","C6","N1"],["C2","O2"],["C4","O4"],["C5","C5M"]],#T
    # アミノ酸は、index 0が主鎖。index 1がカルボニル基。それ以降は側鎖。
    "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"]

#関数
# 立体的なエッジ(円柱)を追加するメソッド。ストランドの半径と色を指定
def addStrand(strand,r,color)
    
    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 #監修 @scivola 先生
        next unless x.typename == "Face"
        result = x.classify_point(strand[0])
        next unless result == Sketchup::Face::PointInside
        connected_entities = x.all_connected
        for y in connected_entities
            y.material= color if y.typename == "Face"
        end
    end
end

def drawMonomers(hash,monomerType)
    #monomerTypeによって、主鎖以外の場所は描かないようにするなど変更したい
    if monomerType == "DNA"
        startIndex = 2
    elsif monomerType == "Protein"
        startIndex = 1
    else
        startIndex = 0
    end
    #monomerTypeはDNAかProteinのどちらか
    hash.each{|key,value|
        #最初のデータのindex 3からmonomerの種類を取得
        symbol = value[0][3]
        #それぞれの行にあるATOMの名前と、座標をHashにしまう。
        atomsAndCoordinates = Hash.new
        for data in value
            atomsAndCoordinates[data[2]] = Geom::Point3d.new(data[6].to_f, data[7].to_f, data[8].to_f)
        end
        #ヌクレオチドごとに、描画すべき結合を、edgeとして追加する。
        #原子情報が欠けている場合があるので、bondsから2つずつ取り出して、対応する座標データが両方得られる場合に限って、edgeを追加する。
        # DNAの場合はindex = 2から、Proteinの場合はindex 1から描画
        list = $monomersAndBonds[symbol][startIndex..-1]
        for bonds in list
            for index in 0..bonds.size - 2
                atomA = atomsAndCoordinates[bonds[index]]
                atomB = atomsAndCoordinates[bonds[index+1]]
                if  atomA != nil and atomB != nil
                    $entities.add_edges([atomA,atomB])
                else
                    puts "結合情報がありません" + key +" "+ symbol + " "+ bonds[index] + "-" + bonds[index + 1]
                end
            end
        end
    }
end

def makeSphere(center,radius)
    #第1引数:中心座標 Geom::Point3d あるいは[1,1,5] 3つの値の配列でも良い
    #第2引数:半径 数値

    #ベクトルを2つ使用。2つが直行していれば何でも良い。
    vector1 = Geom::Vector3d.new(1,0,0)
    vector2 = Geom::Vector3d.new(0,0,1)
    #円1を描く
    circleEdges = $entities.add_circle center, vector1, radius
    #円1に面をつける
    face = $entities.add_face(circleEdges)
    #円2を描く
    circleEdges2 = $entities.add_circle center, vector2, radius*2
    #円1のフェイスに対して、円2をエッジとしてフォローミーツールを適用
    face.followme(circleEdges2)
    #余計なエッジを消す
    for edge in circleEdges2
        edge.erase!
    end
end

def type(lines)
    #最初の1データで判断。
    if ["A","C","G","T"].include?(lines[0][3])
        return "DNA"
    elsif ["ALA" ,"CYS" ,"ASP" ,"GLU" ,"PHE" ,"GLY" ,"HIS" ,"ILE" ,"LYS" ,"LEU" ,"MET" ,"ASN" ,"PRO" ,"GLN" ,"ARG" ,"SER" ,"THR" ,"VAL" ,"TRP" ,"TYR"].include?(lines[0][3])
        return "Protein"
    else
        return "UnKnown Type"
    end
end
    
#ファイルを読み込む
lines = IO.readlines(inputFilePath)

#ストランドごとに行データを分ける。ATOMの書いてある行だけにする。
polymersAndLines = Hash.new{ |hash, key| hash[key] = [] }# 存在しないkeyを指定した時に、そのkeyに関連づけられた空の配列が返るカッコイイhash
for line in lines
    dataInLine = line.split
    if dataInLine[0] == "ATOM"
        key = dataInLine[4]
        polymersAndLines[key].push(dataInLine)
    end
end


colorIndex = 0
mainStrandThickness = 0.1 # 0なら描画しない

polymersAndLines.each{|key,lines|
    #data typeを判別
    type = type(lines)
    if !["DNA","Protein"].include?(type)
        put "unknown polymer type. Polymer: " + key
        next
    end
    mainStrandPoints = []
    monomersAndLines = Hash.new{ |hash, key| hash[key] = [] }# 存在しないkeyを指定した時に、そのkeyに関連づけられた空の配列が返るカッコイイhash
    for dataInLine in lines
        #主鎖情報をkeep
        if typesAndMainstrand[type].include?(dataInLine[2])
            mainStrandPoints.push(Geom::Point3d.new(dataInLine[6].to_f, dataInLine[7].to_f, dataInLine[8].to_f))
        end
        #各モノマーの情報をhashでkeep
        monomerIdentifier = dataInLine[3..5].join("_")
        monomersAndLines[monomerIdentifier].push(dataInLine)
    end
    #主鎖を描画
    if mainStrandThickness != 0
        addStrand(mainStrandPoints,mainStrandThickness,colors[colorIndex%14])#立体的にしないなら
    else
        $entities.add_edges(mainStrandPoints)
    end
    
    #各モノマーを描画
    drawMonomers(monomersAndLines,type)
    #主鎖の色を変更
    colorIndex += 1
}




```
#動作確認
逆転写酵素とDNAの複合体です。下の方にあるのがDNAです。タンパク質分子は2つですが、DNAが付いているのは片方だけのようです。これぐらいのサイズで私のパソコンでは描画まで2分程度でした。


<img width="483" alt="image.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/595771/0d49c3b9-b5ef-e2ea-596a-056c9e16fe9a.png">


#終わりに
作り始めた当初コロナウィルス のキャプシドタンパク質、スパイクタンパク質をモデリングしてリアルなコロナウィルス を作ろう！と思っていたので、まだゴールには達していません。また、タンパク質の表面を描画しようと思っていたのですが、これは原子のデータから表面データへと変換しなくてはならないので、また別のテーマということで今回は保留です(空間に格子状にたくさんの点を置いて、タンパク質までの距離が一定の範囲内にある点だけを残して、これをつないだら表面になるかしら？？)。

コロナ自粛で家にこもっていたら、楽しくなってきてはかどりました。Rubyはほぼ初めてでしたが @scivola さんからたくさんのアドバイスをいただき、素人ながら、それほど素人っぽくないコードになったと思います(少なくとも素人から見たら素人のコードには見えないかも！)。改めてお礼申し上げます。