PyData入門 ######################################## Pythonでデータ解析をする時の、全体的な話です。今後適宜、資料を改変していきます。 `SPARQLthonでの発表 `_ での発表がもとになっています。 準備 **************************************** 以下のモジュールを使えるようにしておきます。 これから環境を構築するのであれば、これらすべてを一度にインストールできる、`Anaconda `_ の利用が便利です。 - `IPython notebook `_ Pythonインタラクティブシェルの機能強化版。コードと結果を1つにまとめて保存できるので便利。 - `matplotlib `_ グラフ描画のためのライブラリ - `pandas `_ 表やベクトル形式のデータを扱うための、多彩で高性能な機能が満載のライブラリ - `scikit-learn `_ 機械学習アルゴリズムを手軽に使えるライブラリ Pandasを使ってみる **************************************** データの読み込みと整形 ======================================== .. figure:: num_data.png CSVデータをエクセルで表示 このデータは、肝臓ガンの患者さんの遺伝子発現データです。行方向に遺伝子が並んでいます。列には、患者さんが並んでいますが、背景が灰色になっている列は、患者さんに関係なく、遺伝子の付加的な情報です。タブ区切りのテキストデータになっているこのファイルを、pnadasで読み込むには、次のようにします。 .. code-block:: python import pandas as pd _exp_data = pd.read_csv('pivot.txt',sep='\t',index_col=2) ファイル名の後に、区切り文字を指定します。index_colで、読み込んだあと、行の名前になる列を指定できます。関数の戻り値は、pandas.DataFrame型で、ちょうどエクセルの1つのシートと同じような表形式のデータになります。 .. code-block:: python _exp_data IPython notebookでは、これだけで簡単に中身を確認することができます。 列名を取得する。 .. code-block:: python _exp_data.columns Index([u'ucscid', u'chrom_pos', u'chrom', u'start', u'end', u'Genesymbol', u'UCSCid', u'category', u'mappedBases', u'exonLength', ... u'M_stR_538N', u'M_stR_548N', u'M_stR_553N', u'M_stR_576N', u'M_stR_581N', u'M_stR_612N', u'M_stR_627N', u'M_E_54N_DR', u'M_E_712N_DR', u'M_E_742N_DR'], dtype='object', length=203) 行の名前を取得する。 .. code-block:: python _exp_data.index Index([u'DDX11L1', u'DDX11L9', u'DKFZp434K1323', u'FLJ00075', u'F379', u'OR4F5', u'DQ595736', u'AK125248', u'DQ597235', u'LOC100132287', ... u'XLOC_014495', u'XLOC_014497', u'XLOC_014499', u'XLOC_014500', u'XLOC_014502', u'XLOC_014507', u'LINC00281', u'XLOC_014503', u'XLOC_014504', u'XLOC_014505'], dtype='object', name=u'Gene', length=38368) どちらも、pandasのIndex型です。 特定の列にアクセスする。 .. code-block:: python _exp_data['chrom'] または _exp_data.chrom 行を指定するときは、ix .. code-block:: python _exp_data.ix['DDX11L1'] 列名で全体をソートすることも可能。 .. code-block:: python _exp_data.sort('chrom') 列名のリストを渡せば、複数の列でソートすることも可能。 .. code-block:: python _exp_data.sort(['chrom','Genesymbol']) 数字だけのデータにしたいので、何列目からデータが入っているのか確認。 .. code-block:: python for i,v in enumerate(_exp_data.columns): print('{}\t{}'.format(i,v)) .. code-block:: python # Rでもお馴染みのスライス sample_ids = _exp_data.columns[12:] 12列目からの列名をsample_idsに格納。 .. code-block:: python # 列名で絞り込んで新たなDataFrameをexp_dataと命名 exp_data = _exp_data[sample_ids] 列の名前が、M_stR_75TやM_stR_135Tなどとなっていますが、最後の１文字がTの時は肝臓がんのデータ、Nのときは正常細胞のデータなので、この列名を、(75,T)や（135,T)などのタプルに変換。データ特有の変換なので、あまり気にしなくてもOKです。 .. code-block:: python def rename(c, _info, maps): for v in maps: if c in _info[v[0]].values: return (_info[(_info[v[0]] == c)].index[0],v[1]) raise ValueError('{}は指定の列の範囲に見付かりません。'.format(c)) from functools import partial exp_rename = partial(rename, _info=info, maps=[('dRNA_T','T'),('dRNA_N','N')]) exp_data.columns = list(map(exp_rename, exp_data.columns)) このように、列名や行名は、上書きできます。 ちょっとした統計処理 ======================================== .. code-block:: python # サンプル(列）ごとの平均 exp_data.mean() (75, T) 11.351534 (135, T) 12.739927 (146, T) 11.572280 (151, T) 12.490965 (153, T) 12.663369 ... .. code-block:: python # 遺伝子（行）ごとの平均 exp_data.mean(1) Gene DDX11L1 0.040685 DDX11L9 0.040093 DKFZp434K1323 11.384697 FLJ00075 10.763521 OR4F5 0.012345 ... 引き数を省略すると0の意味。行列なので、引き数が0の時は、行方向に計算が進む。なので、最終的に列の平均が得られる。列方向に計算が進めば、行の平均が得られるという訳です。他の問う計量も計算できます。stdで標準偏差など。 .. code-block:: python #グラフを描画するための準備 %matplotlib inline import pylab .. code-block:: python # 遺伝子ごとの平均のヒストグラム pylab.hist(exp_data.mean(1)) (array([ 3.83460000e+04, 1.30000000e+01, 2.00000000e+00, 1.00000000e+00, 4.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 2.00000000e+00]), array([ 0. , 2580.05383807, 5160.10767614, 7740.1615142 , 10320.21535227, 12900.26919034, 15480.32302841, 18060.37686648, 20640.43070454, 23220.48454261, 25800.53838068]), .. figure:: hist01.png .. code-block:: python # 出力を変数で受け取って # binを細かく _val = pylab.hist(exp_data.mean(1),bins=50) .. figure:: hist02.png .. code-block:: python # 平均がすごいデカイ奴がいる。 # 誰？ exp_data.mean(1) > 5000 Gene DDX11L1 False DDX11L9 False DKFZp434K1323 False FLJ00075 False F379 False .. code-block:: python # 条件にあう行（遺伝子）だけが抽出できる exp_data[exp_data.mean(1) > 5000] .. code-block:: python #すべてのデータをlogに変換 exp_log = exp_data.applymap(pylab.log10) # 自分用の関数を作る my_log = lambda x:pylab.log10(x) if x > 1e-6 else pylab.log10(1e-6) exp_log = exp_data.applymap(my_log) 0がマイナス無限大になってしまうので、適当なところで、打ち止めするために関数を自作しています。 .. code-block:: python # 遺伝子ごとの平均のヒストグラム _val = pylab.hist(exp_log.mean(1),bins=50) .. figure:: hist03.png .. code-block:: python #行の平均が-4より大きいものだけを残す exp_log_m4 = exp_log[exp_log.mean(1)>-4] _val = pylab.hist(exp_log_m4.mean(1),bins=50) .. figure:: hist04.png ヒストグラムからデータがフィルタリングされているのが分かります。 .. code-block:: python # 数が多いとデモデータとしては面倒なので、800行に削減 import random idx = random.sample(exp_data_m4_s05.index,800) data = exp_data_m4_s05.ix[idx] データの可視化 **************************************** .. code-block:: python # ヒートマップ表示 pylab.matshow(data) .. figure:: heat01.png .. code-block:: python # データの標準化 # 遺伝子ごとに、平均値を引いて標準偏差で割る _mean = data.mean(1) _std = data.std(1) ndata = data.subtract(_mean,0).div(_std,0) # 色付けの範囲をvmin、vmaxで指定できる pylab.matshow(ndata,vmin=-2,vmax=2) .. figure:: heat02.png .. code-block:: python # ボックスプロット(最初の5サンプル) _val = pylab.boxplot(ndata.values[:,0:5]) .. figure:: box.png .. code-block:: python # 散布図（0番目と1番目、標準化前） pylab.plot(data.iloc[:,0], data.iloc[:,1], 'go') .. figure:: plot01.png 3つ目の引き数goは、greenで○の意味。詳しくは、matplotlibのドキュメントを参照してください。 データ解析 **************************************** 統計処理、PCA、K-meansクラスタリングなどをやってみます。 がんと正常組織で発現差のある遺伝子を見つけるために、まずはサンプルを分類して、別々のリストを作る。 .. code-block:: python tumor = [] normal = [] for v in ndata.columns: if v[1] == 'T': tumor.append(v) elif v[1] == 'N': normal.append(v) else: print(v) 何も出力されなければ、うまく分類できた。 .. code-block:: python # t検定 import scipy.stats as ss ss.ttest_ind(data[tumor].ix['SYT3'], data[normal].ix['SYT3']) (-2.0378779730286061, 0.042955723848096199) 出力はタプル。t統計量、P値の順。 主成分分析 ======================================== .. code-block:: python # Scikit-learnを使ったPCA from sklearn.decomposition import PCA # 高次元データを、2次元に縮約 pca = PCA(n_components=2) # 行方向にサンプルが並ぶように転置 pca_data = pca.fit(ndata.T).transform(ndata.T) # 結果のプロット pylab.plot(pca_data[:,0],pca_data[:,1],'go') .. figure:: pca01.png .. code-block:: python # がんのサンプルを赤く、正常組織を青く for i,v in enumerate(ndata.columns): if v[1] == 'T': pylab.plot(pca_data[i][0], pca_data[i][1],'ro') else: # Normal pylab.plot(pca_data[i][0], pca_data[i][1],'bo') .. figure:: pca02.png K-meansクラスタリング ======================================== .. code-block:: python # K-meansクラスタリング from sklearn.cluster import KMeans # 先ほどのPCAと似たコードで実行できる。2クラスに分類 kmeans = KMeans(init='k-means++',n_clusters=2) pred = kmeans.fit_predict(ndata.T) print(pred) [1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 1 1 0 1 1 0 1 1 1 1 0 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 1 0 1 1 0 0 0 1 1 1 1 0 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] .. code-block:: python # がんと正常がちゃんと分類できているか？（ちょっと高度な書き方） from collections import Counter Counter(zip(pred, [v[1] for v in ndata.columns])) Counter({(1, 'T'): 89, (0, 'T'): 71, (0, 'N'): 31}) クラスタリングで1のクラスに入ったサンプルは、全部がんですが、0のクラスにサンプルが混在。2クラスにわけるのは失敗です。 階層的クラスタリング ======================================== 階層的クラスタリングを実行できるモジュールはありますが、樹形図とヒートマップを組み合わせる部分は自作が必要です。このコードはちょっと作りかけな感じがあるので、再利用するときは、ご注意ください。（近日中に直したいと思ってはいますが・・・） .. literalinclude:: cluster.py :linenos: .. code-block:: python _d = hc(ndata.T,'temp.png') .. figure:: hcluster.png まとめ **************************************** PyDataの紹介として、データの読み込みから、簡単なデータ解析までを概観しました。ちょっと突貫なので、暇を見つけてリニューアルします（たぶん） `お気軽にお問い合わせください `_