LoginSignup
35

More than 5 years have passed since last update.

Python+matplotlibでDendrogram付きHeatmap

Last updated at Posted at 2015-03-09

Rではheatmap(x)と打つだけで、x, y各軸に対してクラスタリングされたデンドログラム付きヒートマップが描画でき、遺伝子の発現量が・・・など日々呟いているbioinformaticiansにはお馴染みの機能なのですが、Python+matplotlib環境で同じことをやろうとした際、あまり情報が無かったので、試行錯誤してみました。

備忘録ついでに公開しときます。

データの準備

まずは、適当なデータとして、細胞組織ごとの遺伝子発現量の増減を表したようなデータを作ります。
ここでの組織のチョイスも、遺伝子のチョイスも、全て適当です。
後でクラスタリングしたときにそれっぽく見えるように、乱数で生成したデータを若干小細工してpandas.DataFrameオブジェクトにします。

#!/usr/bin/env python3

genes = [
    'HIST1H4H', 'SPRN', 'DNLZ', 'PYCARD', 'PRDM11', 'DKKL1', 'CYBRD1', 'DMAP1',
    'MT1E', 'PDGFRL', 'SERTM1', 'PIFO', 'FAM109A', 'ST5', 'ABCA7', 'FAM160A1',
    'SAMD15', 'NUAK1', 'GLTP', 'HIST3H2A', 'SCN5A', 'PNPLA5', 'SFRP5', 'CCBE1',
    'PTCD1', 'RFTN1', 'SYTL2', 'FAM65B', 'NFKBIZ', 'RHOG', 'KIF3A', 'DYRK1B',
    'NXPH2', 'APLN', 'ZNF526', 'NRIP3', 'KCNMA1', 'MTSS1', 'ZNF566', 'TNC',
    'GPX2', 'AQP3', 'TSACC', 'SNX15', 'TRIM22', 'THAP6', 'GRIP1', 'DLGAP3',
]
tissues = [
    'brain', 'kidney', 'lung', 'heart',
    'liver', 'pancreas', 'testis', 'placenta',
]


def draw_heatmap(a):
    print(a.head())

def _main():
    from random import choice
    import numpy as np
    from pandas import DataFrame
    v = np.random.random([4, len(tissues)]) * 1.2 - 0.6
    w = np.random.random([3, len(genes)]) * 0.8 - 0.4
    v = np.vstack([choice(v) for x in genes])
    w = np.vstack([choice(w) for x in tissues]).T
    a = DataFrame(v + w, index=genes, columns=tissues)
    draw_heatmap(a)

if __name__ == '__main__':
    _main()
実行結果
             brain    kidney      lung     heart     liver  pancreas  \
HIST1H4H -0.085630  0.074054  0.058026 -0.142751 -0.767515 -0.348885   
SPRN     -0.424203  0.251821 -0.012052 -0.037645 -0.000477  0.714727   
DNLZ      0.372402  0.532086  0.097971 -0.102806 -0.727570  0.109148   
PYCARD   -0.561378  0.114647  0.706732  0.681139  0.718306  0.577552   
PRDM11   -0.698969 -0.022945  0.240133  0.214540  0.251708  0.439960   

            testis  placenta  
HIST1H4H  0.324709 -0.319531  
SPRN     -0.815679 -0.010529  
DNLZ      0.402956 -0.241284  
PYCARD   -0.728135  0.077015  
PRDM11   -0.773637  0.031513  

これ以降、draw_heatmap(a)関数のみをいじっていきます。

手っ取り早く表示する場合

資料作成などではなく、とりあえずデータだけ確認したいときは、次のようにmatplotlibimshow()にデータを突っ込めば表示してくれます。

def draw_heatmap(a):
    from matplotlib import pyplot as plt
    plt.imshow(a, aspect='auto', interpolation='none')
    plt.colorbar()
    plt.xticks(range(a.shape[1]), a.columns)
    plt.yticks(range(a.shape[0]), a.index)
    plt.show()

fig01.png

見栄えを良くする

上の例に対して、少し見栄えを良くしていきます。

  • 目盛りは不要なので消します。
  • ラベルが見にくいので、サイズや位置を調整します。
  • 二色法マイクロアレイの名残か、Down regulationが赤、Neutralが黒、Up regulationが緑・・・という色分けが比較的多いようなので、これに倣います。
  • スケールが正負で対称となっていないので、上限と下限を指定します。
from matplotlib.colors import LinearSegmentedColormap
microarray_cmap = LinearSegmentedColormap('microarray', {
    'red': [(0.0, 1.0, 1.0), (0.5, 0.2, 0.2), (1.0, 0.0, 0.0)],
    'green': [(0.0, 0.0, 0.0), (0.5, 0.2, 0.2), (1.0, 1.0, 1.0)],
    'blue': [(0.0, 0.0, 0.0), (0.5, 0.2, 0.2), (1.0, 0.0, 0.0)],
})

def draw_heatmap(a, cmap=microarray_cmap):
    from matplotlib import pyplot as plt
    plt.imshow(a, aspect='auto', interpolation='none',
               cmap=cmap, vmin=-1.0, vmax=1.0)
    plt.colorbar()
    plt.xticks(range(a.shape[1]), a.columns, rotation=15)
    plt.yticks(range(a.shape[0]), a.index, size='small')
    plt.gca().get_xaxis().set_ticks_position('none')
    plt.gca().get_yaxis().set_ticks_position('none')
    plt.show()

fig02.png

赤→黒→緑、と変わるColormapはmatplotlibのプリセットに無いので、LinearSegmentedColormapで自前で定義しています。
0を完全な黒とすると、0近傍の差が明度が低すぎて分かりにくいので、敢えて中間は濃い灰色にしています。

クラスタリングをする

似た発現パターンを持つ遺伝子をクラスタリングすることで、見やすくします。

def draw_heatmap(a, cmap=microarray_cmap):
    from matplotlib import pyplot as plt
    from scipy.spatial.distance import pdist
    from scipy.cluster.hierarchy import linkage, dendrogram

    metric = 'euclidean'
    method = 'average'

    plt.subplot(1, 2, 1)
    ylinkage =linkage(pdist(a, metric=metric), method=method, metric=metric)
    ydendro = dendrogram(ylinkage, orientation='right', no_labels=True,
                         distance_sort='descending')
    a = a.ix[[a.index[i] for i in ydendro['leaves']]]

    plt.subplot(1, 2, 2)
    plt.imshow(a, aspect='auto', interpolation='none',
               cmap=cmap, vmin=-1.0, vmax=1.0)
    plt.colorbar()
    plt.xticks(range(a.shape[1]), a.columns, rotation=15)
    plt.yticks(range(a.shape[0]), a.index, size='small')
    plt.gca().xaxis.set_ticks_position('none')
    plt.gca().yaxis.set_ticks_position('none')
    plt.gca().invert_yaxis()
    plt.show()

Scipyの機能を用いて、階層クラスタリングを行い、デンドログラムと並べて描画します。

fig03.png

  • 距離の定義は、euclideanの他にも様々なものを使用できる。詳しくはScipyのドキュメント参照。
  • 階層クラスタリング方法も、averageの他にも様々なものを使用できる。詳しくはScipyのドキュメント参照。
  • デンドログラム描画後、葉の順番に従いDataFrameの行を入れ替える。
  • クラスタリングはしたいけれど、デンドログラムは要らない場合、dendrogram(..., no_plot=True)を指定すれば良い。
  • y軸に対し、imshow()は左上を原点とするのに対し、dendrogram()は左下を原点とするので、invert_yaxis()を使って軸の方向を入れ替える。

見栄えを良くする

上の例の見栄えを良くします。

  • デンドログラムとヒートマップをつなげて、y軸ラベルは右側に移動します。
  • デンドログラムの枠は要らないので消します。
  • デンドログラムの色も要らないので、黒のみにします。
  • デンドログラムの幅を、もうすこし詰めます。
def draw_heatmap(a, cmap=microarray_cmap):
    from matplotlib import pyplot as plt
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    from scipy.spatial.distance import pdist
    from scipy.cluster.hierarchy import linkage, dendrogram

    metric = 'euclidean'
    method = 'average'

    main_axes = plt.gca()
    divider = make_axes_locatable(main_axes)

    plt.sca(divider.append_axes("left", 1.0, pad=0))
    ylinkage = linkage(pdist(a, metric=metric), method=method, metric=metric)
    ydendro = dendrogram(ylinkage, orientation='right', no_labels=True,
                         distance_sort='descending',
                         link_color_func=lambda x: 'black')
    plt.gca().set_axis_off()
    a = a.ix[[a.index[i] for i in ydendro['leaves']]]

    plt.sca(main_axes)
    plt.imshow(a, aspect='auto', interpolation='none',
               cmap=cmap, vmin=-1.0, vmax=1.0)
    plt.colorbar(pad=0.15)
    plt.gca().yaxis.tick_right()
    plt.xticks(range(a.shape[1]), a.columns, rotation=15, size='small')
    plt.yticks(range(a.shape[0]), a.index, size='small')
    plt.gca().xaxis.set_ticks_position('none')
    plt.gca().yaxis.set_ticks_position('none')
    plt.gca().invert_yaxis()
    plt.show()

fig04.png

  • subplot()の代わりに、make_axes_locatable()を使っています。
  • デンドログラムの枠をset_axis_off()で非表示にし、link_color_funcで線を黒のみにしています。
  • yaxis.tick_right()で、y軸ラベルを右側に持ってきます。このままではカラーバーと重なるので、colorbar(pad=...)として間隔を広げています。

x軸もクラスタリングする。

時系列データなどならば不要ですが、組織ごとのデータなどではx軸もクラスタリングした方がわかりやすい場合もあります。やることは、基本的に今までと同じです。

def draw_heatmap(a, cmap=microarray_cmap):
    from matplotlib import pyplot as plt
    from mpl_toolkits.axes_grid1 import make_axes_locatable
    from scipy.spatial.distance import pdist
    from scipy.cluster.hierarchy import linkage, dendrogram

    metric = 'euclidean'
    method = 'average'

    main_axes = plt.gca()
    divider = make_axes_locatable(main_axes)
    xdendro_axes = divider.append_axes("top", 0.5, pad=0)
    ydendro_axes = divider.append_axes("left", 1.0, pad=0)

    plt.sca(xdendro_axes)
    xlinkage = linkage(pdist(a.T, metric=metric), method=method, metric=metric)
    xdendro = dendrogram(xlinkage, orientation='top', no_labels=True,
                         distance_sort='descending',
                         link_color_func=lambda x: 'black')
    plt.gca().set_axis_off()
    a = a[[a.columns[i] for i in xdendro['leaves']]]

    plt.sca(ydendro_axes)
    ylinkage = linkage(pdist(a, metric=metric), method=method,
                       metric=metric)
    ydendro = dendrogram(ylinkage, orientation='right', no_labels=True,
                         distance_sort='descending',
                         link_color_func=lambda x: 'black')
    plt.gca().set_axis_off()
    a = a.ix[[a.index[i] for i in ydendro['leaves']]]

    plt.sca(main_axes)
    plt.imshow(a, aspect='auto', interpolation='none',
               cmap=cmap, vmin=-1.0, vmax=1.0)
    plt.colorbar(pad=0.15)
    plt.gca().yaxis.tick_right()
    plt.xticks(range(a.shape[1]), a.columns, rotation=15, size='small')
    plt.yticks(range(a.shape[0]), a.index, size='x-small')
    plt.gca().xaxis.set_ticks_position('none')
    plt.gca().yaxis.set_ticks_position('none')
    plt.gca().invert_yaxis()
    plt.show()

fig05.png

あとはsavefig('hoge.eps')などして論文に貼り付ければおしまい。

ところで、dendrogram(orientation=...)の意味が、横向きと縦向きで異なっているような気がするのですが・・・Scipyのバグ?

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
35