生命科学の分野でも、遺伝子発現情報やDNAメチル化情報のクラスター解析に使われているNMF（non-negative matrix factorization）を応用した手法を実装したPython用ライブラリにに、
nimfaがあります。データマイニングソフトの
Orangeにも採用されているようなので、使って見ました。
まず、
こちらのGithubサイトから、レポジトリをzipファイルで取得して解凍。sudo python setup.py installで入れます。
Python Image Library（PIL)が必要なので、なければこちらはpipなどで手軽に入れておきます。
NMFはやはりいくつのクラスに分類するかがもっとも効いてくるパラメータですので、これを2から9でふるサンプルスクリプトを載せておきます。といっても、estimate_rankというメソッドがあるので便利です。データは、行方向に変数、列方向にサンプルが並んだテキストファイルを、ほとんど神といっても過言では無い
PandasのDataFrameとして読み込んでいます。一部、nimfaのサンプルコードを利用しています。
#!/usr/bin/env python
import pandas
import nimfa
import numpy as np
from scipy.cluster.hierarchy import linkage, leaves_list
from matplotlib.pyplot import savefig, imshow, set_cmap
def reorder(C):
"""
Reorder consensus matrix.
:param C: Consensus matrix.
:type C: `numpy.matrix`
http://nimfa.biolab.si/nimfa.examples.all_aml.html
"""
c_vec = np.array([C[i,j] for i in xrange(C.shape[0] - 1) for j in xrange(i + 1, C.shape[1])])
# convert similarities to distances
Y = 1 - c_vec
Z = linkage(Y, method = 'average')
# get node ids as they appear in the tree from left to right(corresponding to observation vector idx)
ivl = leaves_list(Z)
ivl = ivl[::-1]
return C[:, ivl][ivl, :],ivl
def plot(_C, file_name):
"""
Plot reordered consensus matrix.
:param C: Reordered consensus matrix.
:type C: `numpy.matrix`
:param rank: Factorization rank.
:type rank: `int`
http://nimfa.biolab.si/nimfa.examples.all_aml.html
"""
C,idx = reorder(_C)
imshow(np.array(C))
set_cmap("RdBu_r")
savefig(file_name)
def run_nmf(data):
V = np.matrix(data.as_matrix())
model = nimfa.mf(V, method = "lsnmf", max_iter = 200)
est_rank = model.estimate_rank(range=xrange(2,10),n_run=100)
for i in xrange(2,10):
print('rank = {}'.format(i))
print('cophenetic coef = {}'.format(est_rank[i]['cophenetic']))
print('sample prediction')
print(est_rank[i]['predict_samples'][0])
f_name = 'clustered_consensus_matrix_{}.png'.format(i)
plot(est_rank[i]['consensus'],f_name)
return est_rank
if __name__ == '__main__':
data = pandas.read_table('data.txt',index_col=0,na_values=(' '))
run_nmf(data)
Consensus行列を可視化するとこんな感じになります。
