Pythonで行/列名付きのmatrixを出力するいくつかの方法

正攻法, 変わり種など4つ。

• 2014/11/12　追記

　方法3に問題が合ったので修正いたしました。
　それと方法4を追加。

今回出力したいのはこのようなマトリクス

	Ile	Val	Leu	Phe	Cys	Met	Ala	Gly	Thr	Trp	Ser	Tyr	Pro	His	Asp	Glu	Asn	Gln	Lys	Arg
Ile	IleIle	IleVal	IleLeu	IlePhe	IleCys	IleMet	IleAla	IleGly	IleThr	IleTrp	IleSer	IleTyr	IlePro	IleHis	IleAsp	IleGlu	IleAsn	IleGln	IleLys	IleArg
Val	ValIle	ValVal	ValLeu	ValPhe	ValCys	ValMet	ValAla	ValGly	ValThr	ValTrp	ValSer	ValTyr	ValPro	ValHis	ValAsp	ValGlu	ValAsn	ValGln	ValLys	ValArg
Leu	LeuIle	LeuVal	LeuLeu	LeuPhe	LeuCys	LeuMet	LeuAla	LeuGly	LeuThr	LeuTrp	LeuSer	LeuTyr	LeuPro	LeuHis	LeuAsp	LeuGlu	LeuAsn	LeuGln	LeuLys	LeuArg
Phe	PheIle	PheVal	PheLeu	PhePhe	PheCys	PheMet	PheAla	PheGly	PheThr	PheTrp	PheSer	PheTyr	PhePro	PheHis	PheAsp	PheGlu	PheAsn	PheGln	PheLys	PheArg
Cys	CysIle	CysVal	CysLeu	CysPhe	CysCys	CysMet	CysAla	CysGly	CysThr	CysTrp	CysSer	CysTyr	CysPro	CysHis	CysAsp	CysGlu	CysAsn	CysGln	CysLys	CysArg
Met	MetIle	MetVal	MetLeu	MetPhe	MetCys	MetMet	MetAla	MetGly	MetThr	MetTrp	MetSer	MetTyr	MetPro	MetHis	MetAsp	MetGlu	MetAsn	MetGln	MetLys	MetArg
Ala	AlaIle	AlaVal	AlaLeu	AlaPhe	AlaCys	AlaMet	AlaAla	AlaGly	AlaThr	AlaTrp	AlaSer	AlaTyr	AlaPro	AlaHis	AlaAsp	AlaGlu	AlaAsn	AlaGln	AlaLys	AlaArg
Gly	GlyIle	GlyVal	GlyLeu	GlyPhe	GlyCys	GlyMet	GlyAla	GlyGly	GlyThr	GlyTrp	GlySer	GlyTyr	GlyPro	GlyHis	GlyAsp	GlyGlu	GlyAsn	GlyGln	GlyLys	GlyArg
Thr	ThrIle	ThrVal	ThrLeu	ThrPhe	ThrCys	ThrMet	ThrAla	ThrGly	ThrThr	ThrTrp	ThrSer	ThrTyr	ThrPro	ThrHis	ThrAsp	ThrGlu	ThrAsn	ThrGln	ThrLys	ThrArg
Trp	TrpIle	TrpVal	TrpLeu	TrpPhe	TrpCys	TrpMet	TrpAla	TrpGly	TrpThr	TrpTrp	TrpSer	TrpTyr	TrpPro	TrpHis	TrpAsp	TrpGlu	TrpAsn	TrpGln	TrpLys	TrpArg
Ser	SerIle	SerVal	SerLeu	SerPhe	SerCys	SerMet	SerAla	SerGly	SerThr	SerTrp	SerSer	SerTyr	SerPro	SerHis	SerAsp	SerGlu	SerAsn	SerGln	SerLys	SerArg
Tyr	TyrIle	TyrVal	TyrLeu	TyrPhe	TyrCys	TyrMet	TyrAla	TyrGly	TyrThr	TyrTrp	TyrSer	TyrTyr	TyrPro	TyrHis	TyrAsp	TyrGlu	TyrAsn	TyrGln	TyrLys	TyrArg
Pro	ProIle	ProVal	ProLeu	ProPhe	ProCys	ProMet	ProAla	ProGly	ProThr	ProTrp	ProSer	ProTyr	ProPro	ProHis	ProAsp	ProGlu	ProAsn	ProGln	ProLys	ProArg
His	HisIle	HisVal	HisLeu	HisPhe	HisCys	HisMet	HisAla	HisGly	HisThr	HisTrp	HisSer	HisTyr	HisPro	HisHis	HisAsp	HisGlu	HisAsn	HisGln	HisLys	HisArg
Asp	AspIle	AspVal	AspLeu	AspPhe	AspCys	AspMet	AspAla	AspGly	AspThr	AspTrp	AspSer	AspTyr	AspPro	AspHis	AspAsp	AspGlu	AspAsn	AspGln	AspLys	AspArg
Glu	GluIle	GluVal	GluLeu	GluPhe	GluCys	GluMet	GluAla	GluGly	GluThr	GluTrp	GluSer	GluTyr	GluPro	GluHis	GluAsp	GluGlu	GluAsn	GluGln	GluLys	GluArg
Asn	AsnIle	AsnVal	AsnLeu	AsnPhe	AsnCys	AsnMet	AsnAla	AsnGly	AsnThr	AsnTrp	AsnSer	AsnTyr	AsnPro	AsnHis	AsnAsp	AsnGlu	AsnAsn	AsnGln	AsnLys	AsnArg
Gln	GlnIle	GlnVal	GlnLeu	GlnPhe	GlnCys	GlnMet	GlnAla	GlnGly	GlnThr	GlnTrp	GlnSer	GlnTyr	GlnPro	GlnHis	GlnAsp	GlnGlu	GlnAsn	GlnGln	GlnLys	GlnArg
Lys	LysIle	LysVal	LysLeu	LysPhe	LysCys	LysMet	LysAla	LysGly	LysThr	LysTrp	LysSer	LysTyr	LysPro	LysHis	LysAsp	LysGlu	LysAsn	LysGln	LysLys	LysArg
Arg	ArgIle	ArgVal	ArgLeu	ArgPhe	ArgCys	ArgMet	ArgAla	ArgGly	ArgThr	ArgTrp	ArgSer	ArgTyr	ArgPro	ArgHis	ArgAsp	ArgGlu	ArgAsn	ArgGln	ArgLys	ArgArg

方法1: 大人しくnumpy使う

正攻法 (numpy入ってる人にとっては)。

import numpy as np

Aminos = ["Ile", "Val", "Leu", "Phe", "Cys", "Met", "Ala", "Gly", "Thr", "Trp", "Ser", "Tyr", "Pro", "His", "Asp", "Glu", "Asn", "Gln", "Lys", "Arg"]

data = [[a1+a2 for a2 in Aminos] for a1 in Aminos]
with open("matrix.txt", "w") as f:
        f.write("\t" + "\t".join(Aminos) + "\n")
        np.savetxt(f, np.hstack([zip(Aminos), data]), fmt='%s', delimiter="\t")

方法2: とりあえずprintしてawkに投げる

awkに馴染みがあれば、まぁ楽といえば楽か

Aminos = ["Ile", "Val", "Leu", "Phe", "Cys", "Met", "Ala", "Gly", "Thr", "Trp", "Ser", "Tyr", "Pro", "His", "Asp", "Glu", "Asn", "Gln", "Lys", "Arg"]
print
for a1 in Aminos:
    print a1
for a1 in Aminos:
    print a1
    for a2 in Aminos:
        print a1+a2

% python matrix_test.py | awk 'BEGIN{c=1}{if(c%21 == 0){print $0}else{printf $0"\t"};c++}' > matrix.txt

方法3: printによる出力

見た目上は出来てるが実際は出来てない方法、とその解決方法。

Aminos = ["Ile", "Val", "Leu", "Phe", "Cys", "Met", "Ala", "Gly", "Thr", "Trp", "Ser", "Tyr", "Pro", "His", "Asp", "Glu", "Asn", "Gln", "Lys", "Arg"]

print "\t",
for a1 in Aminos:
    if Aminos[-1] == a1:
        print "%s\n" % a1,
    else:
        print "%s\t" % a1,
for a1 in Aminos:
    print "%s\t" % a1,
    for a2 in Aminos:
        if Aminos[-1] == a2:
            print "%s\n" % (a1+a2),
        else:
            print "%s\t" % (a1+a2),

問題点1

なぜ問題かは次のようにすれば分かる

for i in range(10):
    print "#",

% python print_test.py
# # # # # # # # # #

つまりprint,には改行コードの代わりに空白が入る。
これでは見かけは出来ていても再利用する際に空白とタブの混在でsplitが面倒になる。

これはPython3系であれば

for i in range(10):
    print("#", end="")

Python2系の場合は

import sys
for i in range(10):
    sys.stdout.write("#")

とすることで解消できる。

問題点2

今回の場合は出来ているが、listは重複を許すので

for a1 in Aminos:
    if Aminos[-1] == a1:
        print "%s\n" % a1,
    else:
        print "%s\t" % a1,

ではなく

for i in range(len(Aminos)):
    if len(Aminos) == i+1:
        print "%s\n" % Aminos[i],
    else:
        print "%s\t" % Aminos[i],

の方が確実かと。

方法4: 普通にwrite

Aminos = ["Ile", "Val", "Leu", "Phe", "Cys", "Met", "Ala", "Gly", "Thr", "Trp", "Ser", "Tyr", "Pro", "His", "Asp", "Glu", "Asn", "Gln", "Lys", "Arg"]

with open("matrix.txt", "w") as f:
    f.write("\t")
    for i in range(len(Aminos)):
        if len(Aminos) == i+1:
            f.write("%s\n" % Aminos[i])
        else:
            f.write("%s\t" % Aminos[i])
    for j in range(len(Aminos)):
        f.write("%s\t" % Aminos[j])
        for i in range(len(Aminos)):
            if len(Aminos) == i+1:
                f.write("%s\n" % (Aminos[j]+Aminos[i]))
            else:
                f.write("%s\t" % (Aminos[j]+Aminos[i]))

参考

How to print in Python without newline or space? - Stack Overflow

【Fedora 20】erlangのinstallに失敗する

% sudo yum install -y erlang

.
.
.

失敗:
  erlang-erts.x86_64 0:R16B-03.7.fc20                                                                   

完了しました!
% erl
zsh: command not found: erl

selinux-policyをupdateする必要があるらしい

sudo yum update -y selinux-policy

失敗したerlang-ertsだけ再インストール

sudo yum install -y erlang-erts

参考

Bug 1154327 – SELinux is preventing /usr/lib64/erlang/erts-5.10.4/bin/erlexec from 'execute' accesses on the file .

【JavaScript】地図上の任意の位置に円グラフを表示する

こんな感じで円グラフを表示する方法です。

http://www.cdc.gov/dhdsp/data_statistics/ems/fs_ems_mt.htm

作ってみて、とりあえず楽にカスタムできそうなやつが出来たのでコードをおいておきます。
結果としてはこんな感じ

jVectorMap and Raphaël - CodePen

Raphaëlは初めて使いましたがこれはいいですね、いじってて楽しいのでもう少し触ってみようかと思います。

　
コードなどは私のCodePenにあります。
今回は単純にこうできる、という紹介ですのであまりそれっぽくないですが...
外枠を太くしたり、といったことがg.Raphaelのpiechartに見た感じなさそうだったので大きさの同じcircleを作成してそれをアニメーションさせています。

あと、gridは座標が分かりやすいようつけているだけです。
　

使用ライブラリ, Version等

• jQuery JavaScript Library v1.10.2
• jVectorMap version 2.0.0
• Raphaël 2.1.2
• g.Raphael 0.51
• g.pie.js

CodePen

jVectorMap and Raphaël - CodePen

docker run ...のオプション指定が面倒になったのでコマンド作った

はじめに

user名とかを使いまわす人向けです。

最近dockerが楽しくて一日中弄ってます。

　
ただ、毎回思うのがオプション指定が多くて面倒

docker run -it -u cordea -w /home/cordea --name hoge hoge/hoge /usr/bin/zsh

　
/usr/bin/zsh に関してはDockerfileで指定すれば問題はない

ただ user と working directory はどうしたものやら
調べてもそもそも指定している人が少ない

というわけでpythonで簡単に書いてみました。

詳しいことはGitHub(下にリンクが貼ってあります)のREADME.mdを見て下さい。
　

使用方法

導入方法等はGitHubで。
ただ、調べればREADME.mdより詳しく書いてあるところがあるかと思います。

　
使用前

docker run -it -u cordea -w /home/cordea --name hoge hoge/hoge /usr/bin/zsh

使用後

dockerun hoge/hoge

もちろん以下のように一部のオプションだけ変更したり、もともとrunコマンドにあるオプションを使用するのも可能です。

dockerun -w / --name huge hoge/hoge

ただ、他のオプションを指定する場合はoptparseの仕様で"-- -a"のようにしないといけないみたいです。
まぁ、他のオプションを多用するならdockerunを書き換えるか、docker runを使うかするのが楽かと思いますが...

　
急ごしらえで作ったものでバグも多いかと思います。

バグを見つけたらご自身で直していただくか、ご連絡頂けますと幸いです。

　
探せば似たようなのはありそうだけどとりあえず満足。

GitHub

CORDEA/DockerFiles · GitHub

echo ${array[0]} ...それ大丈夫？

echo ${array[0]}

で出力される結果がshellによっては必ずしも望んだものにはならないという話

要するに配列のインデックスの開始がshellの種類によって0だったり1だったりするということです。
複数のshellを使う人には当たり前の知識かもしれませんが、私には新鮮だったのでメモ。

sh

sh-3.2$ array=(0 1 2 3)
sh-3.2$ echo ${array[*]}
0 1 2 3
sh-3.2$ echo ${array[0]}
0
sh-3.2$ echo ${array[1]}
1

bash

CORDEA@macrou:~$ array=(0 1 2 3)
CORDEA@macrou:~$ echo ${array[*]}
0 1 2 3
CORDEA@macrou:~$ echo ${array[0]}
0
CORDEA@macrou:~$ echo ${array[1]}
1

zsh

CORDEA@macrou% array=(0 1 2 3)
CORDEA@macrou% echo ${array[*]}
0 1 2 3
CORDEA@macrou% echo ${array[0]}

CORDEA@macrou% echo ${array[1]}
0

csh

[CORDEA@macrou ~]$ set array = ( 0 1 2 3 )
[CORDEA@macrou ~]$ echo ${array[*]}
0 1 2 3
[CORDEA@macrou ~]$ echo ${array[0]}

[CORDEA@macrou ~]$ echo ${array[1]}
0

tcsh

[CORDEA@macrou ~]$ set array = ( 0 1 2 3 )
[CORDEA@macrou ~]$ echo ${array[*]}
0 1 2 3
[CORDEA@macrou ~]$ echo ${array[0]}

[CORDEA@macrou ~]$ echo ${array[1]}
0

ksh

$ array=(0 1 2 3) 
$ echo ${array[*]}
0 1 2 3
$ echo ${array[0]}
0
$ echo ${array[1]}
1

GC含量を計算する

はじめに

multifasta形式のファイルを投げるとGC含量を計算します。
ATGC以外の文字(Nとか) が入った場合の動作は保証していません。

　
動作確認にはNucleic and/or Amino Acid contentsを使用させていただきました。

2018/7/16 追記

今でもたまに参照されているようなので、少し修正を加えたものを以下に置いておきました

github.com

Mathematica (Wolfram)

fasta := URLDownload["http://example.com/hoge.fasta"]
lines = Import[fasta, {"FASTA", "LabeledData"}]
gcContent = Map[Module[{map = Counts[Characters[Last[#]]]}, First[#] -> Interpreter["PercentFraction"][N[(map["G"] + map["C"]) / Total[map]]]] &, lines]
Do[Echo[Last[x], First[x] <> ": "], {x, prc}]

Python

使いかた

python gc.py sequence.ffn

コード

#!/bin/env python
# encoding:utf-8
#
# Author:   Yoshihiro Tanaka
# Created:  2014-10-09
#
import sys
infile = open(sys.argv[1], "r")
lines = infile.readlines()
infile.close()
atgc = ['A', 'T', 'G', 'C']
fflag = False
lst = []
for l in range(len(lines)):
    line = lines[l]
    if ">" in line:
        if fflag:
            sys.stdout.write(seq)
            print("gc%: " + str( (lst[2]+lst[3]) / float(sum(lst)) ))
            lst = []
        seq = line
        fflag = True
    else:
        for i in range(len(atgc)):
            try:
                lst[i] += line.count(atgc[i])
            except:
                lst.append(line.count(atgc[i]))
sys.stdout.write(seq)
print("gc%: " + str( (lst[2]+lst[3]) / float(sum(lst)) ))

awk

awk version 20070501
　
動作確認はMac OS Xのみ

使い方

awk -f gc.awk sequence.ffn

コード

#!/bin/awk -f
# encoding:utf-8
#
# Author:   Yoshihiro Tanaka
# Created:  2014-10-26
#

function calc_gc(atgc) {
    sum = 0;
    for (i in atgc) {
        sum += atgc[i]
    }
    gc = ((atgc["G"] + atgc["C"]) / sum)*100
    return gc
}
function plus_atgc(atgc) {
    if (substr($0, i, 1) in atgc) {
        return substr($0, i, 1)
    }
}
{
    if ($0~/^>/) {
        if (length(atgc) > 0) {
            printf("%s%0.2f\n","gc%: ",calc_gc(atgc))
        }
        print $0
        for (i=1; i<=4; ++i) {
            atgc[substr("ATGC", i, 1)] = 0
        }
    }
    else {
        for (i=1; i<=length($0); ++i) {
            ++atgc[plus_atgc(atgc)]
        }
    }
}
END {
    printf("%s%0.2f\n","gc%: ",calc_gc(atgc))
}

Bash

動作確認はDebianとCentOSしかしていません。
GNU sed しか動かないようです。BSD sed(Mac OSとか)の方はgsed使うか書き換えて下さい。
　

使い方

chmod u+x gc.sh
./gc.sh sequence.ffn

コード

#!/bin/bash
#
# Author: Yoshihiro Tanaka
# date: 2014-10-08

flag=false
lines=$(wc -l $1 | sed -e 's/^  *//g' | cut -d ' '  -f 1)
j=1
for line in `cat $1`;
    do
        if [ `echo $line | grep ">"` ]; then
            if $flag; then
                echo $seq
                for val in ${array[@]}; do
                    let sum+=$val
                done
                echo "gc%: $(echo "scale=4; ((${array[2]}+${array[3]}) / $sum) * 100" | bc | cut -c 1-5)"
                sum=0
                array=()
            fi
            seq=$line
            flag=true
        else
            tmp=($(echo $line | sed -e 's/\(.\)/\1\n/g' | grep "[ATGC]" | sort | uniq -c | sed -e 's/^  *//g' | cut -d ' ' -f 1))
            i=1
            for val in ${tmp[@]}; do
                let result=${array[$i]}+$val
                array[$i]=$result
                let i++
            done
            if [ $lines -eq $j ]; then
                for val in ${array[@]}; do
                    let sum+=$val
                done
                echo $seq
                echo "gc%: $(echo "scale=4; ((${array[2]}+${array[3]}) / $sum) * 100" | bc | cut -c 1-5)"
            fi
        fi
        let j++
    done

軽く説明

主に

tmp=($(echo $line | sed -e 's/\(.\)/\1\n/g' | grep "[ATGC]" | sort | uniq -c | sed -e 's/^  *//g' | cut -d ' ' -f 1))

このあたりについて。
　

ATGCで並んだ文字列をsedで一文字ずつ改行してgrepでATGCだけ取り出します。
これは放置していると改行コードまで入ってしまうためです。

次のsortとuniqによってATGCの並び順はA, C, G, Tになります。
その後のsedでは前半のスペースを全て消しています。
ちなみにawk '{print $1}' によって以下の部分を省略できます。

sed -e 's/^  *//g' | cut -d ' ' -f 1

今回はawkを使わずに書きたかったのであえてこの書き方にしてあります。

　
そのため、以下の部分では　array[2] (C)　と　array[3] (G)　を取り出して計算しています。

echo "gc%: $(echo "scale=4; ((${array[2]}+${array[3]}) / $sum) * 100" | bc | cut -c 1-5)"

それと、以下の部分でcutしているのは、bcを行った際にどうしても5.5500のように00が付いてしまうのが紛らわしいと感じたためです。

"scale=4; ((${array[2]}+${array[3]}) / $sum) * 100" | bc | cut -c 1-5

何かいい方法があればご教授ください。

【awk】ファイルを複数の区切り文字を使って出力する場合

最近小ネタばっかりですが

こんなファイルがあったとして

0,1,2,3,4,5,6,7,8,9
a,b,c,d,e,f,g,h,i,j

このような感じで出力したい場合

0	1	2	3	4:5:6:7:8:9
a	b	c	d	e:f:g:h:i:j

#!/bin/awk -f

BEGIN {
    FS=","
}
{
    for(i=1; i<=NF; i++) {
        if (i>5) {
            if (i==NF) {
                print $i
            } else {
                printf $i ":"
            }
        } else {
            printf $i "\t"
        }
    }
}

もう少しきれいに書けないものか...