2016-07-01

BioPythonで系統樹を書いてNewick出力

タイトルの通り。
MegaよりもSeaViewよりも軽いし落ちない。

#!/usr/bin/env python
#coding: UTF-8

#よく使うモジュールは何も考えずに書く

import sys 
import re
import os
import glob
import copy
import commands
from bisect import bisect
from Bio import SeqIO
import urllib2
import datetime
import matplotlib
from Bio import Phylo
from Bio.Phylo import NewickIO
from Bio import Entrez
argvs = sys.argv  # コマンドライン引数を格納したリストの取得
argc = len(argvs) # 引数の個数

#name = argvs[1]

from Bio.Phylo.TreeConstruction import DistanceCalculator
from Bio import AlignIO
aln = AlignIO.read('hogehoge.fasta', 'fasta')
#print aln
#SingleLetterAlphabet() alignment with 5 rows and 13 columns
calculator = DistanceCalculator('identity')
dm = calculator.get_distance(aln)
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
constructor = DistanceTreeConstructor(calculator, 'nj')
tree = constructor.build_tree(aln)
print tree
Phylo.write(tree, 'hugahuga.tre', 'newick')

2016-03-14

nt、env_ntデータベースから相同性ベースで配列を取り出す

ある遺伝子ホモログについて、代表配列20個程度がある状態がスタート。

この配列をnt、env_ntデータベースに投げてヒットした領域をmulti fastaとして取り出したい。そんなのしょっちゅうやってるよと思ったが意外とめんどくさかった。

なんでかというとファイルのデカさ。ntより小さめのenv_nt でさえ100GB以上ある。このサイズだと大好きなSeqIO.parseで処理できない。どうしよう。分割して保存したらサーバー上の僕の領域が確実にパンクする。

アプローチを変えてみる。遺伝研スパコンの力でblastはいけたので、ヒットした配列のGenBank IDを取得。http://www.ncbi.nlm.nih.gov/nuccore/**のアドレスからソースを取得し無理矢理配列をゲットしようと思ったがあまり美しくないなと思い直し、biopythonでできないかしらと調べるとコードまで出てきてハッピー。

Fetching Genbank Entries For List Of Accession Numbers.

リンク先のコードのrettype="gb"をrettype="fasta"にすれば万事解決。

2016-01-13

統計の勉強-共分散

統計の問題集をちまちま解いているのですが共分散の所でちょっと詰まった。

問題としては独立な二つの正規分布(norm1, norm2)があり、それらの分散は何れも１。

このとき、norm1と(norm1+norm2)/√2の共分散を求めよ、という感じ。

回答みても直感的に理解できなかったのでrで書いてみる。

#norm1 <- rnorm(10000, mean = 100, sd = 1)　
#norm2 <- rnorm(10000, mean = 1000, sd = 1)

#散布図を作成し平均を描画
plot(norm1,norm2)
abline(h = mean(norm2),v = mean(norm1),col = "red")

まず分散1で適当な正規分布を2つ作成。

散布図はこうなる。赤線はそれぞれの平均。

f:id:kmooog:20160113105141p:plain

ここで相関係数は

"散布図に平均点で縦横に軸を描き込んで、プロットされている各人の点からそれぞれの軸に下ろした垂線と軸とで囲まれた長方形の（正負つきの）面積の平均"

相関係数のわかりやすい解釈（続き） - たけみたの脱社会学日記

なので適当に外れ値選んで長方形を書く。

# 外れ値において長方形作成
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean(norm2),mean(norm2),norm2[which.max(norm1)],norm2[which.max(norm1)]), col="gray")

f:id:kmooog:20160113105350p:plain

これらの分布に相関はないので共分散は0.

次にnorm2を(norm1+norm2)に変更し同様の解析。

# norm1,norm2+norm1で同様の処理
plot(norm1,(norm2+norm1))
abline(h = mean(norm2+norm1),v = mean(norm1),col="red")
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean(norm2+norm1),mean(norm2+norm1),norm2[which.max(norm1)]+norm1[which.max(norm1)],norm2[which.max(norm1)]+norm1[which.max(norm1)]), col="gray")

norm1が足された分、すべての点はY軸方向正の方向へ移動し、相関が出てくる。

f:id:kmooog:20160113105728p:plain

ここで今回の図と最初の図の長方形の面積の差を考える。

最初の図の長方形は

(norm1-norm1の平均)*(norm2-norm2の平均)

二番目の図の長方形は

(norm1-norm1の平均)*{(norm1+norm2)-(norm1の平均+norm2の平均)}

この差分を取るとnorm1^2-norm1の平均^2(=norm1の分散)となる。

よって今回の共分散はnorm1の分散と等しく1。

ここまではOK。次はnorm1と(norm1+norm2)/√2でやってみる

f:id:kmooog:20160113111535p:plain

おっとそうか図が変わらないのか。元々norm1とnorm2に相関無いから共分散に影響を与えるのはこの場合もnorm1の分散で、面積としてはy軸が1/√2になってるので今回の分散はnorm1の分散/√2 = 1/√2。

なんかグラフィカルに理解しないと頭に入ってこない。

確率変数の計算得意になりたい。

#適当な正規分布作成
norm1 <- rnorm(10000, mean = 100, sd = 1)　
norm2 <- rnorm(10000, mean = 1000, sd = 1)


#散布図を作成し平均を描画
plot(norm1,norm2)
abline(h = mean(norm2),v = mean(norm1),col = "red")
which.max(norm1)

# 外れ値において長方形作成
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean(norm2),mean(norm2),norm2[which.max(norm1)],norm2[which.max(norm1)]), col="gray")


var(norm1,norm2)

# norm1,norm2+norm1で同様の処理
plot(norm1,(norm2+norm1))
abline(h = mean(norm2+norm1),v = mean(norm1),col="red")
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean(norm2+norm1),mean(norm2+norm1),norm2[which.max(norm1)]+norm1[which.max(norm1)],norm2[which.max(norm1)]+norm1[which.max(norm1)]), col="gray")



#norm1が足された分、すべての点はY軸方向正の方向へ移動し、相関が出てくる

var(norm1,(norm1+norm2))
sum((norm1-mean(norm1))*( (norm1+norm2)-mean( (norm1+norm2))))/10000



plot(norm1,(norm2+norm1)/1.4142)
abline(h = mean(norm2+norm1)/1.4142,v = mean(norm1),col = "red")
sum((norm1-mean(norm1))*( (norm1+norm2)/1.4142 -mean( (norm1+norm2)/1.4142)))/10000
polygon( c(norm1[which.max(norm1)],mean(norm1),mean(norm1),norm1[which.max(norm1)]) ,c(mean((norm1+norm2)/1.4142),mean(norm1+norm2)/1.4142,(norm2[which.max(norm1)]+norm1[which.max(norm1)])/1.4142,(norm2[which.max(norm1)]+norm1[which.max(norm1)])/1.4142), col="gray")

2015-10-28

今日一番の気付き

今日まではpython スクリプトをqsubで投げる時わざわざ.shファイルを作るようなつまらない人生を送ってきましたが、以下のコマンドで済むことを知ってQOL(クオリティ・オブ・ライフ)が一気に上がりました。

qsub -cwd -S $(which python) hogehoge.py

2015-10-24

blasrのインストールに再挑戦

環境は遺伝研スパコン。

LPM: Local Package Managerにパッケージが上がっていたので試してみた。

hdf5は一発で通るのにblasrは通らない。

---

iblasr/BlasrHeaders.h:22:23: 致命的エラー: libconfig.h: そのようなファイルやディレクトリはありません

---

致命的エラー。そうですか。

仕方ないのでblasrに必要なpbcoreを自力で入れようと試みる。

cd pbcore

python setup.py install

---

The installation directory you specified (via --install-dir, --prefix, or

the distutils default setting) was:

/usr/local/lib/python2.7/site-packages/

---

デフォルトではpythonのsite-packagesに入れることになってるのか？やり直し。

python setup.py install --prefix /home/***

---

running install

Checking .pth file support in /home/***/lib/python2.7/site-packages/

error: can't create or remove files in install directory

The following error occurred while trying to add or remove files in the

installation directory:

[Errno 2] No such file or directory: '/home/***/lib/python2.7/site-packages/test-easy-install-35953.pth'

The installation directory you specified (via --install-dir, --prefix, orthe distutils default setting) was:

/home/***/lib/python2.7/site-packages/

This directory does not currently exist. Please create it and try again, or choose a different installation directory (using the -d or --install-dir option).

---

なんかsite-packagesに入れたがってるっぽい。

ローカルにpython2.7を落としてsite-packageを指定してやり直し。

export PYTHONPATH=/home/***/lib/python2.7/site-packages/

python setup.py install --prefix /home/***

通りそうな雰囲気を出しつつも、

ImportError: No module named Cython.Build

で止まる。えぇー

とりあえずお手上げ。

2015-09-30

盗まれた自転車が戻ってきた話

一昨年の冬、自転車が盗まれた。

夕方にラボから戻って、家でちょっと準備してから夕飯の買い物に行こう、というその１０数分の間鍵をかけずアパートの駐輪場に置いておいたところを盗まれたのだ。

ママチャリではあったが買ったばっかりだったのでそれなりに悲しく、また今住んでいる地域の治安の悪さを実感して嫌な気持ちになる事件だった。

一応自転車屋さんの盗難保険に入っていたので、警察の人に来てもらい盗難届けを出し、保険でかなり割り引かれた状態で新しい自転車を手に入れられたのでそれっきり忘れていた。

しかし３日ほど前、突然電話がかかる。

「もしもし、柏警察署ですがこちら熊谷さんのお電話でよろしかったでしょうか？」

最近警察に迷惑をかけるようなことをした覚えもないし、ここ数年はややこしい人たちとの付き合いもない。即座に自転車が見つかったことだろうと察する。話を聞くと、ごく家の近所で盗まれた自転車が見つかったらしい。

東京に住んでいた頃は、自転車で走っているだけでしょっちゅう呼び止められて登録番号を確認されたものだが、柏に移ってからそのようなことが一切ないので、盗難された自転車が見つかることはまずないだろうと思っていた。柏の外れでも、社会秩序を保とうと努力している人はいるようだ。

「一応盗難事件の証拠品という扱いになりますので、本日現場に来れない場合は柏警察署まで来ていただくことになりますがよろしいでしょうか？」

よろしいも何もそうするしかないんでしょうそうしますよ、と昔さんざん職務質問を食らった経験から来る警察への卑屈さを滲ませながら、後日取りに行く約束をする（その日は学部時代にお世話になった先輩の結婚式の二次会があり、久々に高田馬場の汚い空気を吸いに行くところであった）。

「車ないので取りに行くとしたら乗って帰ることになると思うんですが、自転車は使える状態でしょうか？」

「見たところ特に問題ないですよ！全然乗って帰れると思います！」

そんな顛末でわざわざ柏警察署まで出向き、30分以上かけて乗って帰ってきた自転車がこちら。

f:id:kmooog:20150930152852j:plain

f:id:kmooog:20150930152906j:plain

まごうことなき粗大ゴミである。

もしも僕がピクサーアニメの監督で、自転車を擬人化した作品を作っていて、レイプされた後のママチャリを表現しなければいけないなら、この自転車は非常にいいモデルになるのではないか。

バイク風になるように、わざわざブレーキの位置を変えているあたり盗んだ犯人の人柄が透けて見えて好ましい。呼び鈴も外されてるけど、お巡りさんに問題なく乗れると言われたので、道路交通法や条例などの面でもセーフという判断でいいんですね？お巡りさん？これで注意されたらあなたの名前出しますよ？

とりあえず、ラボ内の移動用などに自転車欲しかったのでちょっと整備して使う（呼び鈴はつけます）。

ピエール瀧になりたい

備忘録として