xlsの日付はfloat型

最近はxlsでダウンロードしてきたデータをxlrdで読みだしてSQLAlchemy経由でoracleにimportするようなものを作っている。

Excelで日付が文字列として入っている場合は strptime, fromtimestamp, mktimeを組み合わせてOracleに突っ込めばいい。

普通にExcelのDate型の場合、実態はfloatらしいので読み出すときにちょっと困るが、Date型ではいってる場合はxldate_as_tupleっていうxlrdのユーティリティ関数を使って

publication_date = datetime(*xldate_as_tuple(sheet.cell(r, 21).value, wb.datemode))

みたいにできる。

でもExcelなのでというか、parse出来ないよエラー(floatの値が何故か小さすぎる?)とかValueErrorとか普通に出るというか、元データがダーティ過ぎてそっちの対応に四苦八苦してたw

ggplotっぽい配色

ggplotっぽい配色が欲しかったので探したらsofで見つけたがこれはRのコードなのでPythonで

colorsysというモジュールを使う。

でもこれは0-1の範囲の数字がタプルで返ってきて使いにくいので255をかけた後16進数表記にして返す関数を書いた。

from colorsys import hls_to_rgb

def color_list(n, l=0.65, s=1.0):
    return ['#%02x%02x%02x' % tuple(map(lambda x: int(x * 255), \
        hls_to_rgb((i + 3 - n)/ float(n), l, s))) for i in range(n)]

10個ほど生成させてみるといい感じ

Bioinformatics版のProject Euler

RosalindというProject Eulerっぽいサイトがあったので少し遊んでみた。

楽しくbioinformaticsの問題を解きながらアルゴリズムを覚えていく感じなので、BioPythonを覚えるのにちょうどいいかなと思った。

GC Contentを計算するメソッドくらいあるだろうとは思うんだけど知らないので素朴に計算。リスト内包表記でmaxで取り出すのがいいかなと思ったけど、面倒くさいのでこれまたループで素朴に。

import sys
from Bio import SeqIO

def count_gc(seq):
    return float((seq_record.seq.count("G") + seq_record.seq.count("C"))) * 100 / len(seq_record.seq)

if __name__ == '__main__':
    fasta = sys.argv[1]
    highest_id = ""
    highest_content = 0
    for seq_record in SeqIO.parse(fasta, "fasta"):
        c = count_gc(seq_record)
        if c > highest_content:
            highest_id = seq_record.id
            highest_content = c

    print highest_id
    print "{0:.6f}".format(highest_content)

僕はBioPerlでバイオインフォマティクス関連のプログラミングをしてたので、BioPythonで配列解析をしたことは殆ど無いんですよね。

あらかじめクラスタ数を決めないでクラスタリングする方法(Affinity Propagation)

K-meansのように予めクラスタ数を指定すると、「そのクラスタ数は正しいの?」っていう疑問が浮かぶと思う。

「なんらかの統計値に基づいて適切なクラスタに分割して欲しい」そんな願いを叶えるのがAffinity Propagationというクラスタリングアルゴリズムである

exemplara(セントロイドとかクラスタ中心)になるべきパラメータ(responsibility)とクラスタメンバに属しやすさ(availability)を交互に更新していって収束させる手法なので、K-meansのような初期値依存性がないらしい。

クラスタ数は類似度行列の対角要素(自分との類似度)に依存する(デフォルトはmedian)のでここを変更するとクラスタ数も変わるんだけどね。

Scikit-learnではAffinity Propagationが実装されているのでsykのケミストリースペースを作ってクラスタリングしてみた。ちなみにスライドのPCAの説明は間違っていた(pca.fit(fps).transform(fps)としなければいけなかった)。

from rdkit import Chem
from rdkit.Chem import AllChem, DataStructs
from sklearn.decomposition import PCA
from sklearn.cluster import AffinityPropagation
from ggplot import *
import numpy as np
import pandas as pd

suppl = Chem.SDMolSupplier('syk.sdf')

fps = []
for mol in suppl:
    fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2)
    arr = np.zeros((1,))
    DataStructs.ConvertToNumpyArray(fp, arr)
    fps.append(arr)

print len(fps)
pca = PCA(n_components=2)
x = pca.fit(fps).transform(fps)
af = AffinityPropagation().fit(x)

d = pd.DataFrame(x)
d.columns = ["PCA1", "PCA2"]
d["labels"] = pd.Series(af.labels_)
g = ggplot(aes(x="PCA1", y="PCA2", color="labels"), data=d) + geom_point() + xlab("PCA1") + ylab("PCA2")
ggsave("ap.png", g)

1404818907

RDKitでFatal Python error: PyThreadState_Get: no current thread

homebrewのRDKit入れたらこんなエラーに遭遇

$ python mkdesc.py 
Fatal Python error: PyThreadState_Get: no current thread
Abort trap: 6

調べたらboost関連らしい

brew rm boost
brew install --with-icu --build-from-source boost
brew uninstall rdkit
brew install rdkit

で動くようになったけど brew installで再度boostのインストールを始めたので2番目のboostのインストールは必要ないかもしれない

土曜日はpython-ggplotの話もします

ChEMBLからデータを取ってきて、活性が強いほど色が濃くなってドットのサイズも大きくなるようにしてみた。python-ggplotを使うと数行で描くことができるので調子いいですみたいな話をします。

あとはscikit-learnとrdkitで予測モデルを作ったりクラスタリングをしたりといった話もできればなーと考えているので、 興味があれば参加して下さい。

import pandas as pd
from ggplot import *
import numpy as np

d = pd.read_csv("syk.csv")
d["pIC50"] = 9 - np.log10(d["IC50"])
print d

p = ggplot(aes(x='MWT', y='ALOGP', color="pIC50", size="pIC50"), data=d) + geom_point()
#print p
ggsave("2dplot.png", p)

1404120520

scikit-learnで交差検定が地味に便利

train_test_splitっていうメソッドが便利。これでテストデータと訓練データを分けてくれる。

SVC with linear kernelでちょっと検定

>>> import numpy as np
>>> from sklearn import cross_validation
>>> from sklearn import datasets
>>> from sklearn import svm
>>> iris = datasets.load_iris()
>>> X_train, X_test, y_train, y_test = cross_validation.train_test_split(
...     iris.data, iris.target, test_size=0.4, random_state=0)
>>> X_train.shape, y_train.shape
((90, 4), (90,))
>>> X_test.shape, y_test.shape
((60, 4), (60,))
>>> clf = svm.SVC(kernel='linear').fit(X_train, y_train)
>>> clf.score(X_test, y_test)
0.96666666666666667

あるディレクトリにあるファイルのなかで一番新しく作られたものを調べる(Python)

globとosモジュールを使う

>>> import os, glob
>>> max(glob.iglob('*.txt'), key=os.path.getctime)
'list_sample.txt'

参考

ProductName みんなのPython 第3版
柴田 淳
SBクリエイティブ株式会社 / ?円 ( 2013-12-04 )


CheckiOはPython3の学習サイト

遊んでみたけど結構楽しい

checkio

pandasとggplotを学ぶための本

最近グラフはPandas+ggplotで描いていて、深夜のバッチ処理で毎日沢山のグラフを生成させている。

rpy2?

使ってないなーみたいな。

Python版ggplotはR版のAPIと同じものを提供することをゴールにしているのでggplotのクックブックを読んでおくといいことがあります。六角形の密度プロットとか早く実装されないかなー。

ProductName Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集
Winston Chang
オライリージャパン / 3570円 ( 2013-11-30 )


尚、下の本も参考になります。

ProductName グラフィックスのためのRプログラミング
H. ウィッカム
丸善出版 / 4200円 ( 2012-02-29 )


そして、データの操作はPandasで。

ProductName Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理
Wes McKinney
オライリージャパン / 3780円 ( 2013-12-26 )


機械学習はscikit-learnを使えばよい。