Python で Wheel Diagram を描く

タンパク質の二次構造である α-ヘリックスの概要を掴むために、Wheel Diagram というものがある。読んで字の如く、螺旋状のヘリックスを長軸方向から見て車輪のようにプロットするというもの。α-ヘリックスのどの面にどんなアミノ酸が存在するかがひと目で分かるので便利。これを描くためにいつも使っていたサーバーが繋がらなくなっていたので、自前で作ってみた。

まずは、疎水性パラメータを取得する関数を作成してみた。
使用した疎水性パラメータは Kyte and Doolittle, 1982 のもの。膜貫通領域の予測などでも使われている。
この関数では、[配列のインデックス、周期、残基、疎水性] と疎水性モーメント、疎水性モーメントの大きさを出力する。

import math

def wheel_diagram(seq):
    result=[]
    RESIDUES = ["A","L","I","V","F","Y","W","M","C","G","P","N","Q","S","T","H","R","K","D","E"]
    kyte_doolittle = [1.8,3.8,4.5,4.2,2.8,-1.3,-0.9,1.9,2.5,-0.4,-1.6,-3.5,-3.5,-0.8,-0.7,-3.2,-4.5,-3.9,-3.5,-3.5]
    cycle = 0
    moment = [0,0]
    for index, res in enumerate(seq):
        result.append([index, cycle, res, kyte_doolittle[RESIDUES.index(res)]])
        moment[0] = moment[0] + (kyte_doolittle[RESIDUES.index(res)] * math.cos(index*100/180*math.pi))
        moment[1] = moment[1] + (kyte_doolittle[RESIDUES.index(res)] * math.sin(index*100/180*math.pi))
        if (index+1)%18 == 0 and index is not 0:
            cycle += 1
    #print(result)
    #print(moment)
    #print(math.sqrt(moment[0]**2+moment[1]**2))
    return result, moment, math.sqrt(moment[0]**2+moment[1]**2)

wheel_diagram("LFTIDSLLDDLRNTMKVAVHIIQQHTVLI")

実行してみると以下の通り。

([[0, 0, 'L', 3.8], [1, 0, 'F', 2.8], [2, 0, 'T', -0.7], [3, 0, 'I', 4.5], [4, 0, 'D', -3.5], [5, 0, 'S', -0.8], [6, 0, 'L', 3.8], [7, 0, 'L', 3.8], [8, 0, 'D', -3.5], [9, 1, 'D', -3.5], [10, 1, 'L', 3.8], [11, 1, 'R', -4.5], [12, 1, 'N', -3.5], [13, 1, 'T', -0.7], [14, 1, 'M', 1.9], [15, 1, 'K', -3.9], [16, 1, 'V', 4.2], [17, 1, 'A', 1.8], [18, 2, 'V', 4.2], [19, 2, 'H', -3.2], [20, 2, 'I', 4.5], [21, 2, 'I', 4.5], [22, 2, 'Q', -3.5], [23, 2, 'Q', -3.5], [24, 2, 'H', -3.2], [25, 2, 'T', -0.7], [26, 2, 'V', 4.2], [27, 3, 'L', 3.8], [28, 3, 'I', 4.5]], [4.110011070569221, -34.87093799688161], 35.11231276607343)

疎水性モーメントが 35.1 と大きめ（20を超えると大きい印象がある）なので、両親媒性ヘリックスであることがわかる。（追記：この値は長さで割ってない…）
これを可視化するためにプロットしてみる。

import matplotlib.pyplot as plt

def plot_wheel_diagram(seq):
    initial_radius = 3
    marker_size = 1
    result, moment_vec, moment_val = wheel_diagram(seq)
    prev_point = [initial_radius,0]
    for res in result:
        if res[2] in ["R", "K", "H"]:
            color="blue"
        elif res[2] in ["D", "E"]:
            color="red"
        elif res[2] in ["N", "Q", "S", "T"]:
            color="green"
        else:
            color="black"
        current_point = [math.cos(res[0]*100/180*math.pi)*(initial_radius+res[1]*marker_size),math.sin(res[0]*100/180*math.pi)*(initial_radius+res[1]*marker_size)]
        plt.text(current_point[0], current_point[1], res[2], color=color, fontsize=16, zorder=999)
        
        plt.plot([prev_point[0],current_point[0]],[prev_point[1],current_point[1]], color="gray", linewidth=0.2, zorder=100)
        prev_point = current_point
    final_radius = (initial_radius+result[-1][1]*marker_size)*1.2
    #plt.plot([0,moment_vec[0]/10],[0,moment_vec[1]/10])
    plt.text(0, 0.2, "{:.3g}".format(moment_val),horizontalalignment='center', zorder=20)
    an = plt.annotate("", xy=[moment_vec[0]/10,moment_vec[1]/10], xytext=[0,0], arrowprops=dict(shrink=0, width=1, headwidth=8,headlength=8, connectionstyle='arc3',facecolor='lightgray', edgecolor='lightgray'))
    an.set_zorder(10)
    plt.xlim(-final_radius, final_radius)
    plt.ylim(-final_radius, final_radius)
    plt.xticks([])
    plt.yticks([])
    plt.gca().invert_yaxis()
    plt.gca().set_aspect("equal")
    plt.show()

plot_wheel_diagram("LFTIDSLLDDLRNTMKVAVHIIQQHTVLI")

実行してみると以下の通りプロットされる。
今回は塩基性アミノ酸を青字で、酸性アミノ酸を赤字で、中性アミノ酸を緑字で、その他疎水性アミノ酸を黒字で表示した。疎水性の強い方向に矢印が向いて表示されている。

ブラウザ版: Plot wheel diagram

Python版: Github

めでたしめでたし。