ランダムな塩基配列をつくる

長さは20-merで、ATGCの4文字を0-3の整数に割り当てておいて、あとは乱数を引いてdictionaryにアクセスする。

簡単で多分早い（はず）。random.randint(a, b)は、a以上b以下なのでbを含む。numpyとは異なるので注意されたし。

miRBaseのID変換

miRBase

　miRBase (http://www.mirbase.org/)は、microRNAのおそらく世界でもっとも使われているデータベース。さまざまな生物種のmiRNAについて、そのアノテーションや配列などの情報がのっている。

HTSeq

　HTSeqは、DESeqなど発現解析をするさいにraw read countをBAMファイルから計算するときに使われるhtseq-countというコマンドラインツールがあるが、このツールはHTSeqという名前空間をもったpython libraryとして提供されている。pip install htseqでインストールできる。

　SAM/BAM/GFF/GTF/Fasta/FastaqなどBioinformatics解析では、デファクトとなっているフォーマットのパーサーがあるので、非常に便利。Reference overview — HTSeq 0.6.1p2 documentation このURLには、各ファイル形式にたいするparser classがのっているので、これに任意のファイルパスを与えれば、適当にパースされた結果がイテレータとして返るつくりになっている。

　ちなみに、HTSeq.GFF_Readerクラスは、GenomicFeatureのイテレータを返す。GenomicFeatureクラスはパースしたGFFの情報がいろいろつめ込まれている。HTSeq.GenomicFeatureはiv (interval)というattributeをもち、ここにはstrandやchromosome、startやendなどの情報がはいっている。

miRBaseのID

　miRBaseにはいくつかIDがある。miRNAはその生合成過程において、二本鎖のpre-miRNAが切断されて5pと3pの二つのmature miRNAができる。3pと5pはお互いに相補的な配列を持ち、どちらか/両方がターゲットの遺伝子と塩基配列相補的に対合する。

　miRBaseでは、preはMIというプレフィックスが、matureはMIMATというプレフィックスがそれぞれユニークなIDに割り当てられている。すなわち、おなじpreからできた3pと5pは元は同じなので、共通のpreのIDへ遡ることができる。

mature miRNAからpre-miRNAへ

　これを踏まえて、任意のmature miRNA IDがどのpre-miRNA IDから生まれたのか？を知るためには、以下のようなコードになる。つまり、preとmatureのGTFファイル (なんでもよいが)を準備して、mature GTFに含まれている"Derives_from"属性をキーにして、pre.gtfから該当するカラムを見つける。これにより、もともとの (pre)のIDのほかにゲノム座標などを切り出すことができる。

ggplot2のプロットをできるだけ共通化したい

小さな問題

　ggplot2の場合であるがplot系の関数の使われ方として、おなじデータにたいして、異なるスケールで、ラベルを変えて、色を変えて、、と幾つかのバージョンのplotを作りたいことは多々ある。こういうのは、g.1 <- ggplot(...)、g.2 <- ggplot(...)などと、データやカラム名にべったりなコードを書いてしまって、後からも非常に読みにくい。とてもつらい。できるところは共通化して、ラベル名やカラム名は引数として渡しておけばすっきりするように思える。

　ただ、文字列としてカラム名を渡しても、エラー: geom_point requires the following missing aesthetics: x, yとかでエラーになる。これまでめんどくさくて調べずに放置していたが、ちゃんと方法があった。

ggplot2::aes_string()を使おう

　通常は、aes(x=y, y=y)であるところを、ggplot2::aes_stringでwrapしてあげることで、引数として与えたカラム名をよしなにやってくれる。コードに起こすと以下のようになる。

　こんな感じで、同じデータフレームで、カラム名が異なるplotを一つの関数としてwrapすることで、3つのplotを一つの関数へ共通化することができた。あと、可変長引数となっているので、ggplot::titleなどもどんどん渡していく。まあこれくらいのことなら、pairとかいうカラムにplotの種類を列挙して、ggplot::facet_wrap(~pair)一発でもよい気がするが、、、。

　それから、引数として受け取ったデータフレームdfのxカラム加工する場合、df$x <- log10(df$x+1)とかだと、'closure' 型のオブジェクトは部分代入可能ではありませんとか怒られるので、df[, x] <- log10(df[, x]+1)と書くこと。ハマった。R力高めたい。なんだかtipsというか、ただのベストプラクティスとなってしまった。

emacs+pyflakes+flymakeでpythonの構文チェックを勝手にしてもらう

リアルタイムに構文チェックしたい

flymakeについては、flymake (いままでこれ無しでどうやってプログラム書いてたんだろう) - にゃあさんの戯言日記 ここを参考にされたし。emacsで書いてる傍から、リアルタイムに構文チェックが走るというやつ。メジャーな言語なら多分あると思う。

pythonとか動的型付け言語だと、コンパイルフェーズがないから、変数宣言を忘れたり、構文エラーでも実行しないと、それがエラーなのかわからないので困る。特に、重い計算するときに、途中まで処理うまくいったのに、ある関数を呼び出したら変数が未定義でエラーとかつらい。なので、事前に潰しておきたいというモチベーション。ちなみに、emacs 24.4.1を対象としています。このくらい新しいバージョンのemacsだと、flymakeはもう入ってるからすぐ使える。

pyflakesのインストール

pip pyflakesすればよいです。which pyflakesでpathをメモっておく。pyenvでインストールしたpythonだと、$HOME/.pyenv/shims/pyflakesにあると思う。

emacsの設定

(add-hook 'find-file-hook 'flymake-find-file-hook)
(when (load "flymake" t)
  (defun flymake-pyflakes-init ()
    (let* ((temp-file (flymake-init-create-temp-buffer-copy
                       'flymake-create-temp-inplace))
           (local-file (file-relative-name
                        temp-file
                        (file-name-directory buffer-file-name))))
      (list "PATH_TO_YOUR_pyflakes"  (list local-file))))
  (add-to-list 'flymake-allowed-file-name-masks
               '("\\.py\\'" flymake-pyflakes-init)))
; show message on mini-buffer
(defun flymake-show-help ()
  (when (get-char-property (point) 'flymake-overlay)
    (let ((help (get-char-property (point) 'help-echo)))
      (if help (message "%s" help)))))
(add-hook 'post-command-hook 'flymake-show-help)

PATH_TO_YOUR_pyflackesを書き換えると、pyflakesがホックされてemacs上でほぼリアルタイムに入力に対して、チェッカーが動く。flymake-show-helpという関数が定義されているが、この関数がflymakeのメッセージをemacsのミニバッファに表示してくれるので便利。

様子

importしたのに使ってないとか、定義してない変数ですよとか親切。

networkxでグラフを描画するとmatplotlibのエラーでラベルが表示できない

まえがき

networkxというグラフ理論を扱うpythonのライブラリがある。networkxのグラフオブジェクトは、networkx.draw_networkx_nodes()やnetworkx.draw_networkx_edges()などのメソッドでグラフ構造を可視化することができる。

import matplotlib.pylab as plt
import networkx

def show(g):
    text_font = "sans-serif"
    node_color = 'blue'
    node_alpha = 0.4
    pos = networkx.pygraphviz_layout(g, prog="circo")
    networkx.draw_networkx_nodes(g, pos, node_color="pink", alpha=node_alpha+0.5)
    networkx.draw_networkx_edges(g, pos, edge_color="blue", alpha=node_alpha, arrows=False)
    networkx.draw_networkx_labels(g, pos, font_size=12, font_family='sans-serif')
    plt.show()

if __name__ == '__main__'
    G = networkx.Graph() # 無向グラフ
    G.add_nodes_from(['A', 'C', 'B', 'E', 'D', 'G', 'F', 'I', 'H']) # nodeを作る
    G.add_edges_from([('A', 'I'), ('A', 'C'), ('A', 'B'), ('C', 'F'), ('B', 'D'), ('E', 'F'), ('D', 'G'), ('G', 'H')]) # edgeを作る
    show(G) 

問題

nodesやedgesにint/float型ではない、str型を渡すと以下のようなエラーがでる。グラフは描画できるが、ノードに何も表示されないのでつらい。

  File "/Users/yukke/.pyenv/versions/2.7.8/lib/python2.7/site-packages/matplotlib/artist.py", line 59, in draw_wrapper
    draw(artist, renderer, *args, **kwargs)
  File "/Users/yukke/.pyenv/versions/2.7.8/lib/python2.7/site-packages/matplotlib/figure.py", line 1079, in draw
    func(*args)
  File "/Users/yukke/.pyenv/versions/2.7.8/lib/python2.7/site-packages/matplotlib/artist.py", line 59, in draw_wrapper
    draw(artist, renderer, *args, **kwargs)
  File "/Users/yukke/.pyenv/versions/2.7.8/lib/python2.7/site-packages/matplotlib/axes/_base.py", line 2092, in draw
    a.draw(renderer)
  File "/Users/yukke/.pyenv/versions/2.7.8/lib/python2.7/site-packages/matplotlib/artist.py", line 59, in draw_wrapper
    draw(artist, renderer, *args, **kwargs)
  File "/Users/yukke/.pyenv/versions/2.7.8/lib/python2.7/site-packages/matplotlib/text.py", line 538, in draw
    bbox, info, descent = self._get_layout(renderer)
  File "/Users/yukke/.pyenv/versions/2.7.8/lib/python2.7/site-packages/matplotlib/text.py", line 311, in _get_layout
    ismath=False)
  File "/Users/yukke/.pyenv/versions/2.7.8/lib/python2.7/site-packages/matplotlib/backends/backend_macosx.py", line 166, in get_text_width_height_descent
    six.text_type(s), family, size, weight, style)
ValueError: failed to convert font family name to ASCII

matplotlibのbackendを変えればよい

matplotlibは、GTKAggやTkAggなど、異なるグラフィックエンジンを使用できる。$HOME/.matplotlib/~/matplotlibrcにその設定を書くことでバックエンドを切り替えることができる。どうも、デフォルトのバックエンドが問題で、backend : TkAggを設定ファイルに追記したら、治った。

今度は、ちゃんとnodeにラベルが表示されている。めでたい。

Upgrading to emacs-24.4 from 21.X, threw a "Symbol's function definition is void" error

問題

emacs-21.Xを24.4 (current version)にhomebrew経由でアップグレードした時に、Symbol's function definition is void:package-desc-versというエラーがでて困った。package.elは、emacsの設定ファイルをマネージメントするやつで、確か24.4から標準で入っている。

調査

• grep -R package-desc-vers ~/emacs.d/すると、package.el、melpa.el関連で幾つかpackage-desc-vers()関数をcallしているファイルが分かった。
• Symbol's function definition is void: package-desc-vers · Issue #110 · purcell/emacs.d · GitHub 同じ状態のissueもある。

解決

package-desc-vers()はpackage--ac-desc-version()に置き換わっていた模様。emacs/package-x.el at master · emacs-mirror/emacs · GitHubを参照。

なので、古いやつをpackage--ac-desc-version()へと関数名をちまちま書き換えたら、動いた。editorの設定するのめんどくさいし、ずっと動いててほしいから、後方互換を捨てないでほしい。せめて、deprecated warningを出す、というおもてなしが欲しかった。

ノートパソコン.new

　Macbook Pro 15 inch (retina)が届いた。メリットは画面が綺麗で早い、デメリットは重いです。

　あと、クリーンな環境で必要なものを手作業で入れていったけど、2015年を目前にして、人類がzshやemacsの設定ファイルをちまちま移動させたり、直したりしているのはおかしい。でも、OS入れなおしてから2時間くらいでほぼ昔の環境が出来上がって安心感を得た。