X分で学ぶJulia
- 2016/11/20 追記: この記事は内容が一部古くなっています。最新版はこちらを参照して下さい
Julia Advent Calendar 2014 1日目の記事です。
まだまだJuliaは一般に知られていない言語ということもありまして、Juliaの基本的な機能を最短で学ぶ記事を初日に書いた次第です。 X<30くらいだと思います。
バージョンはJuliaの最新リリース版であるv0.3系を基にしていますが、特に異なる点は次期バージョンであるv0.4に関しても触れています。
Juliaの公式サイトは The Julia Language で 、Juliaのバイナリは Julia Downloads から入手できます。
文法
百聞は一見にしかず。まずはJuliaのコードをざっと見てみましょう。
function mandel(z) c = z maxiter = 80 for n = 1:maxiter if abs(z) > 2 return n-1 end z = z^2 + c end return maxiter end function randmatstat(t) n = 5 v = zeros(t) w = zeros(t) for i = 1:t a = randn(n,n) b = randn(n,n) c = randn(n,n) d = randn(n,n) P = [a b c d] Q = [a b; c d] v[i] = trace((P.'*P)^4) w[i] = trace((Q.'*Q)^4) end std(v)/mean(v), std(w)/mean(w) end
From: http://julialang.org/
どうでしょう？PythonやRubyをやったことがる人なら初見でも大体の意味が分かるのではないでしょうか?
関数定義・分岐・反復などの構文はそれぞれfunction ... end, if ... end, for ... end, while ... endのようにキーワードで始まりendで終わります。
ちょうどRubyと同じような感じです。
インデントはPythonのように必要ではありませんが、4スペースでひとつのインデントを表すのが慣習です。
また、変数の宣言や型の指定は通常必要ありません。
こうしたコードはJuliaのLLVMベースのJITコンパイラによりコンパイルされ、C言語などで書いたコードとそれほど変わらない速度で実行できます。
変数
変数は特別に宣言せずとも初期化と同時に使用できます。
functionやforなどほとんどのブロックは変数のスコープを新たに作ります。これはPythonなどと動作が異なりますので注意が必要です。例外的にスコープを作らないのはif ... endとbegin ... endです。
for i in 1:10 x = i end println(x) # ERROR: x not defined
if true x = 10 end println(x) # OUTPUT: 10
トップレベルで定義された変数はグローバル変数になります。
グローバル変数はコンパイラの最適化を難しくするため、constをつけて定義すべきです。
数値型
Juliaの値は型を持ちます。Juliaでは動的に型がつき、様々な機能と密接に関わってきます。
整数型はsigned, unsignedと8bit, 16bit, 32bit, 64bit, 128bitの組み合わせとBool型とBigInt型で合計12種類あり、
それぞれsigned 64bitはInt64やunsigned 32bitはUint32など一貫性のある型名がつけられています(v0.4でUintはUIntに改名)。
浮動小数点の型も16bit, 32bit, 64bitとBigFloat型で合計4種類があります。
BigInt型とBigFloat型はそれぞれ任意精度の整数と浮動小数点数です。
他には複素数のComplex{T}型があり、Tという 型パラメータ(type parameter) で実部と虚部の数値の型を指定します。ちょうどHaskellの型変数(type variable)のようなものです。
科学計算のために作られているJuliaは、このように豊富な数値の方を持つ点が魅力のひとつです。
リテラル
大抵、何らかの値を作るリテラルは他の言語と同じです。
- 数値
- 文字(列)
- 文字
Char型:'a','樹' - 文字列
ASCIIString型:"deadbeef","""Triple Quote String"""
- 文字
- その他
- 真偽値
Bool型:true,false - シングルトン
Nothing型:nothing(v0.4ではVoid型)
- 真偽値
IntはInt32またはInt64のエイリアスです。
Juliaの対話実行環境(REPL)で確認してみましょう。
型はtypeof関数で確認できます。
~/s/JuliaAdvent $ julia
_
_ _ _(_)_ | A fresh approach to technical computing
(_) | (_) (_) | Documentation: http://docs.julialang.org
_ _ _| |_ __ _ | Type "help()" for help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 0.3.0 (2014-08-20 20:43 UTC)
_/ |\__'_|_|_|\__'_| |
|__/ | x86_64-apple-darwin13.3.0
julia> 42
42
julia> typeof(42)
Int64
julia> typeof(1.0)
Float64
julia> '漢'
'漢'
julia> typeof('漢')
Char
julia> typeof("deadbeef")
ASCIIString (constructor with 2 methods)
julia> typeof(true)
Bool
julia> typeof(nothing)
Nothing (constructor with 1 method)
ベクトルは[]で記述できます。
Juliaではインデックスは1始まりなので注意が必要です。
julia> [1,2,3]
3-element Array{Int64,1}:
1
2
3
julia> x = [1,2,3]
3-element Array{Int64,1}:
1
2
3
julia> x[1]
1
以降対話セッションは>で表現します。
型
Juliaの値はひとつの 具体的な 型を持ちます。一部例外を除いて、型が自動的に別の型にキャストされることはありません
(一部例外とは、1 + 1.5などの数値計算とコンストラクタです。 http://julia.readthedocs.org/en/latest/manual/conversion-and-promotion/)。
ここで、 具体的(concrete)な とわざわざ述べたのは、Juliaにはこれと対極をなす 抽象的(abstract)な 型があるためです。 適切な訳語がわからないため、ここでは原文通り具体的な型を「具体型」、抽象的な型を「抽象型」と表記します。
最も大きな違いは具体型はインスタンス化可能な一方で抽象型はインスタンス化ができない点です。
即ち、任意のxに対して、typeof(x)の結果は必ず具体型になります。
抽象型は具体型や他の抽象型のsupertypeとして機能し、型間のsubtype/supertypeの関係性は木構造のグラフをとります。
さらに、具体型は他の型のsupertypeにはなれませんので、必然的にグラフの葉が具体型に、内部ノードが抽象型となります。
このグラフの根にあるのがAnyという抽象型です。
Juliaでユーザーが型を定義するとデフォルトでこのAny型のsubtypeになります。
(厳密には、Noneという具体型を含む全ての型のsubtypeになっている抽象型があります。)
具体例で確認しましょう。Int64とInt32は共に具体型で、Integerという抽象型のsubtypeになっています。
具体型はisleaftype、subtype/supertypeの関係は<:という関数で確認できます。
> isleaftype(Int64) true > isleaftype(Int32) true > isleaftype(Integer) false > Int64 <: Integer true > Int32 <: Integer true
自分で型を作るのも簡単です。
以下のようにtypeまたはimmutableに続けて型名を書き、フィールドを定義します。
type Person name::UTF8String age::Int end immutable Point x::Float64 y::Float64 end
デフォルトコンストラクタがありますので、即座にインスタンス化できます。
インスタンスのフィールドへは.でアクセスできます。
person = Person("佐藤建太", 24) point = Point(1.0, 2.0) println("Hello, my name is $(person.name) and $(person.age) years old.")
関数/メソッド
Juliaの関数定義には二種類の書き方があります。
手続き型言語っぽい通常の関数定義と、関数型言語っぽい=を用いた記法です。
# 手続き型風 (複数行のときに好まれる) function add_one(x) x + 1 end # 関数型風 (一行で書けるときに好まれる) add_one(x) = x + 1
返り値を表すreturnはPerlやRubyのように省略可能です。
関数の引数は基本的に参照渡しのようになりますが、引数値の変数への束縛を変更できないため、immutableな型(例えばIntやASCIIString)の値は関数内で変更できません。
逆にmutableな型(例えばArrayなどの配列)の値は関数内で変更することができます。
以下の例では標準のsort!関数で破壊的に配列の中身をソートしています。
Juliaでは破壊的な関数の名前に!をつけるのが慣習です。
> x = [4, 1, 2, 3]
4-element Array{Int64,1}:
4
1
2
3
> sort!(x)
4-element Array{Int64,1}:
1
2
3
4
> x
4-element Array{Int64,1}:
1
2
3
4
Juliaの他のスクリプト言語と違う大きな特徴のひとつが関数の 多重ディスパッチ(multiple dispatch) です。 これは、引数の型の組み合わせによって呼び出されるメソッドの実体が異なるという振る舞いを実現できます。
簡単な例で動作を確認しましょう。次の関数fooは4つの異なる実装を持ちます。
function foo(x::Int) println("foo Int: $x") end function foo(x::Int, y::Int) println("foo Int Int: $x $y") end function foo(x::Float64, y::Float64) println("foo Float64 Float64: $x $y") end function foo(x::Int, y::Float64) println("foo Int Float64: $x $y") end
見ての通り、<仮引数>::<型>の形で引数の型を指定できます。
これらの型の異なるメソッドは、呼び出し側の引数の型にマッチするものが実行されます。
もちろん、ユーザが定義した型に対しても同様に振る舞います。
以下のように型にマッチするものがない場合、その関数呼び出しはエラーになります。
> foo(1) foo Int: 1 > foo(1, 2) foo Int Int: 1 2 > foo(1.0, 2.0) foo Float64 Float64: 1.0 2.0 > foo(1, 2.0) foo Int Float64: 1 2.0 > foo(1.0, 2) ERROR: `foo` has no method matching foo(::Float64, ::Int64)
このように、Juliaではメソッドは関数に従属します。これは、メソッドがクラスに従属するクラスベースのプログラミング言語とは大きく異なります。
キーワード引数やオプショナル引数、可変長引数なども使えます。
=の後に値を書くとオプショナルな引数として認識され、引数の区切りを,の代わりに;を使うとそれ以降の引数がキーワード引数として扱われます。
args...とすると、仮引数argsは定義内ではタプルとなります。
# オプショナル引数 head(seq, n::Int=10) = seq[1:n] # キーワード引数 function optimize(func; iter=100, rate=0.1, alg=GradientDescent) # ... end # 可変長引数 sumof(args...) = sum(args)
- http://julia.readthedocs.org/en/latest/manual/functions/
- http://julia.readthedocs.org/en/latest/manual/methods/
構文糖衣
ここではJuliaで見られる特殊な構文をざっと見て行きます。
係数
変数の前に数値をつけると、暗黙的に積だと解釈されます。
これは多項式の記述や12pxなど単位の記述に便利です。
> x = 2.1 2.1 > 2x 4.2 > 4x^2 + 3x + 2 25.94
範囲
配列から一部を取り出したり、ある範囲で反復すときなどに便利なのが範囲型です。
基本的な文法はstart:endもしくはstart:step:endのようにコロンで区切って記述します。
範囲自体もオブジェクトですので、変数に収めたり関数に渡したりもできます。
> 1:4
1:4
> for i in 1:4
println(i)
end
1
2
3
4
> 1:2:10
1:2:10
> for i in 1:2:10
println(i)
end
1
3
5
7
9
配列から一部を取り出すのにも範囲は使用されます。
終端はendで指定できます。
> x = [1,2,3,4,5]
5-element Array{Int64,1}:
1
2
3
4
5
> x[1:3]
3-element Array{Int64,1}:
1
2
3
> x[4:end]
2-element Array{Int64,1}:
4
5
> x[2:end-2]
2-element Array{Int64,1}:
2
3
内包表記
配列や辞書はPythonのような内包表記で記述できます。
> [x for x in 1:4]
4-element Array{Int64,1}:
1
2
3
4
> [c => i for (c, i) in zip('a':'f', 1:6)]
Dict{Char,Int64} with 6 entries:
'd' => 4
'f' => 6
'b' => 2
'e' => 5
'c' => 3
'a' => 1
ベクトル/行列
Juliaにおいてベクトルとは、一次元配列のことです。 ベクトルは以下のように記述します。
> [1, 2, 3]
3-element Array{Int64,1}:
1
2
3
,を抜かすと、2次元の配列(行列)として解釈されます。
> [1 2 3]
1x3 Array{Int64,2}:
1 2 3
複数行の行列の記述は、以下のように;で書けます。
> [1 2 3; 4 5 6]
2x3 Array{Int64,2}:
1 2 3
4 5 6
式の後に'をつけると行列の転置になります。
> [1 2 3; 4 5 6]'
3x2 Array{Int64,2}:
1 4
2 5
3 6
無名関数
->を使うことで、無名関数を記述できます。
> x -> 2x
(anonymous function)
> map(x -> 2x, [1,2,3])
3-element Array{Int64,1}:
2
4
6
ブロック引数
第一引数に関数をとる関数は、位置を逆転させて代わりにdo ... endブロックを取ることができます。
即ち、以下の２つが同じ意味になります。
# 関数引数が先 map(x -> 2x, [1,2,3]) # 関数引数の代わりにブロックをとる map([1,2,3]) do x 2x end
これは、関数の中身が大きい場合に特に便利です。
さらに、ファイルを開いて自動的に閉じる処理の記述にも使われています。 以下のコードはファイルを開いて一行づつ処理をするパターンです。
open("somefile.txt") do f for line in eachline(f) # ... end end
非標準文字列
Juliaでは文字列の前に識別子をおくと、通常の文字列ではなくマクロとして解釈されます。 これを利用して、正規表現を記述できます。
> r"\w-\d \w+"
r"\w-\d \w+"
> match(r"\w-\d \w+", "B-2 Spirit")
RegexMatch("B-2 Spirit")
他にはバージョン文字列の記述にも利用されています。
> v"1.2.3" v"1.2.3"
外部コマンド
バッククォートで囲うことでコマンドを記述できます。
作られたコマンドはrun関数で実行できます。
> `ls -la` `ls -la` > typeof(`ls -la`) Cmd > run(`ls -la`) total 40 drwxr-xr-x+ 4 kenta staff 136 11 30 20:57 . drwxr-xr-x+ 26 kenta staff 884 11 29 17:17 .. -rw-r--r--+ 1 kenta staff 12692 11 30 20:57 learnjulia.md -rw-r--r--+ 1 kenta staff 234 11 29 18:10 sample.jl
便利なマクロ
@マークで始まる式は、コンパイル時に別の式に展開されます。この機能を マクロ といいます。
別の言い方をすれば、通常の関数が実行時に値を取って値を返すのに対し、マクロとはコンパイル時に式を取って式を返す関数と言えます。
以下の例はアサーションに失敗すると例外を投げる@assertマクロと、変数の状態をプリントする@showマクロです。
これら２つはデバッグなどに大変役立ちます。
> x = 100 100 > @assert x < 10 ERROR: assertion failed: x < 10 in error at error.jl:21 > @show x x = 100 100
@printfマクロは書式を指定したプリントのためのマクロです。
> @printf "line %3d: %s\n" 42 "foo" line 42: foo
ベンチマークなどに使えるのが実行時間とアロケートされたメモリを測る@timeマクロです。
最初の呼び出し時にはJITコンパイル分のオーバーヘッドがかかりますので注意が必要です。
> fib(n) = if n < 2 1 else fib(n-1) + fib(n-2) end fib (generic function with 1 method) > @time fib(10) elapsed time: 0.001739938 seconds (47552 bytes allocated) 89 > @time fib(20) elapsed time: 9.9135e-5 seconds (96 bytes allocated) 10946 > @time fib(40) elapsed time: 1.160873359 seconds (96 bytes allocated) 165580141
@code_llvmと@code_nativeは関数呼び出しで実行されるLLVMとネイティブコードをプリントします。
パフォーマンスが上がらない時などに原因を探るのに便利です。
> @code_llvm fib(10)
define i64 @"julia_fib;42387"(i64) {
top:
%1 = icmp sgt i64 %0, 1, !dbg !1280
br i1 %1, label %L, label %if, !dbg !1280
if: ; preds = %top
ret i64 1, !dbg !1282
L: ; preds = %top
%2 = add i64 %0, -1, !dbg !1282
%3 = call i64 @"julia_fib;42374"(i64 %2), !dbg !1282
%4 = add i64 %0, -2, !dbg !1282
%5 = call i64 @"julia_fib;42374"(i64 %4), !dbg !1282
%6 = add i64 %5, %3, !dbg !1282
ret i64 %6, !dbg !1282
}
> @code_native fib(10)
.section __TEXT,__text,regular,pure_instructions
Filename: none
Source line: 1
push RBP
mov RBP, RSP
push R15
push R14
push RBX
push RAX
mov RBX, RDI
cmp RBX, 1
jle L67
Source line: 1
lea RDI, QWORD PTR [RBX - 1]
movabs R15, 4534292768
call R15
mov R14, RAX
add RBX, -2
mov RDI, RBX
call R15
add RAX, R14
L56: add RSP, 8
pop RBX
pop R14
pop R15
pop RBP
ret
L67: mov EAX, 1
jmpq L56
マクロの書き方はちょっと入門の範囲を超えるので省略します。
パッケージ
JuliaのパッケージはGitレポジトリです。Julia本体もパッケージもGitHubを中心に回っています。 現在、Mercurialなど他のバージョン管理システムには対応していません。 しかし、パッケージ開発者を除くユーザとしてはGit自体を学ぶ必要はありません。 パッケージのインストールはJuliaの対話セッションから行うことができます。
> Pkg.add("DocOpt")
INFO: Installing DocOpt v0.0.2
INFO: Package database updated
ほか、レポジトリのメタデータのアップデートはPkg.update()で可能です。
公式に利用可能なパッケージはhttp://pkg.julialang.org/から一覧できます。
それ以外にもPkg.cloneを使うことで、公式レポジトリにはないがGitHubなどのGitレポジトリから直接インストールすることもできます。
その他
この記事では深く扱わなかった重要な部分については、Juliaのマニュアルへのリンクを張っておきます。 Juliaのマニュアルはとても読みやすく、それだけでプログラミング言語の勉強になるほど質の良いものです。
- 例外 http://julia.readthedocs.org/en/latest/manual/control-flow/#exception-handling
- モジュール http://julia.readthedocs.org/en/latest/manual/modules/
- 並列計算 http://julia.readthedocs.org/en/latest/manual/parallel-computing/
- C/Fortranの呼び出し http://julia.readthedocs.org/en/latest/manual/calling-c-and-fortran-code/
- コンストラクタ http://julia.readthedocs.org/en/latest/manual/constructors/
- 標準ライブラリ http://julia.readthedocs.org/en/latest/stdlib/base/
Consensus CDSの領域をGRangesオブジェクトにする
Bioconductor - TxDb.Hsapiens.UCSC.hg19.knownGeneは便利なのですが、ある遺伝子に対応するtranscriptが多すぎてどうしたものかと困ってました。
で、UCSCのゲノムブラウザで見るとConsensus CDS(CCDS)のアノテーションもあるのですが、このBioconductorのパッケージからそれを取得する方法が分かりませんでした。
正確にはmakeTranscriptDbFromUCSC(genome="hg19", tablename="ccdsGene")などとしてCCDSの領域は取得はできるんですが、他のアノテーション(Gene IDなど)が一切消えるので全く使いものにならない感じです。
仕方がないのでCCDS ProjectのFTPサーバから元データを引っ張ってきてそれをパースしてアノテーションのあるCCDSのデータを作ったのでちょろっと共有しておきます。
まずはCCDS ProjectのFTPサーバからCCDSの領域とアノテーションを収めたファイルを取得します (ftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/)。 reference genomeのビルドによっていくつかありますので、自分の目的に合ったものを取得しましょう。 私はftp://ftp.ncbi.nlm.nih.gov/pub/CCDS/archive/Hs37.3/CCDS.current.txtを取りました。
次にこのデータをRがパースしやすいようにPerlで整形します。元データはこんな感じのでデータになってます。
CCDS.current.txt:
#chromosome nc_accession gene gene_id ccds_id ccds_status cds_strand cds_from cds_to cds_locations match_type 1 NC_000001.8 LINC00115 79854 CCDS1.1 Withdrawn - 801942 802433 [801942-802433] Identical 1 NC_000001.10 SAMD11 148398 CCDS2.2 Public + 861321 879532 [861321-861392, 865534-865715, 866418-866468, 871151-871275, 874419-874508, 874654-874839, 876523-876685, 877515-877630, 877789-877867, 877938-878437, 878632-878756, 879077-879187, 879287-879532] Identical 1 NC_000001.10 NOC2L 26155 CCDS3.1 Public - 880073 894619 [880073-880179, 880436-880525, 880897-881032, 881552-881665, 881781-881924, 883510-883611, 883869-883982, 886506-886617, 887379-887518, 887791-887979, 888554-888667, 889161-889271, 889383-889461, 891302-891392, 891474-891594, 892273-892404, 892478-892652, 894308-894460, 894594-894619] Identical 1 NC_000001.10 PLEKHN1 84069 CCDS4.1 Public + 901911 909954 [901911-901993, 902083-902182, 905656-905802, 905900-905980, 906065-906137, 906258-906385, 906492-906587, 906703-906783, 907454-907529, 907667-907803, 908240-908389, 908565-908705, 908879-909019, 909212-909430, 909695-909743, 909821-909954] Identical 1 NC_000001.10 HES4 57801 CCDS5.1 Public - 934438 935352 [934438-934811, 934905-934992, 935071-935166, 935245-935352] Identical 1 NC_000001.10 ISG15 9636 CCDS6.1 Public + 948953 949857 [948953-948955, 949363-949857] Identical 1 NC_000001.10 C1orf159 54991 CCDS7.2 Public - 1018272 1026922 [1018272-1018366, 1019732-1019762, 1019860-1019885, 1021257-1021391, 1022518-1022583, 1022881-1022976, 1025732-1025807, 1026851-1026922] Identical 1 NC_000001.10 TTLL10 254173 CCDS8.1 Public + 1115433 1120521 [1115433-1115719, 1115862-1115980, 1116110-1116239, 1117120-1117194, 1117740-1117825, 1118255-1118426, 1119299-1119470, 1120348-1120521] Identical 1 NC_000001.10 TNFRSF18 8784 CCDS9.1 Public - 1138970 1141950 [1138970-1139339, 1139778-1139865, 1140749-1140871, 1141764-1141950] Identical
CDS領域の範囲が一行に収められてます。 RマスターならRでもサクッといけるんでしょうが、私は以下のPerlスクリプトで事前にTSVファイルに変換しました。 (追記: CCDSの入力ファイルが0-based coordinateだったので1-basedに修正)
parse_ccds.pl:
#!/usr/bin/env perl use v5.10; use strict; use warnings; say "chromosome\tnc_accession\tgene\tgene_id\tccds_id\tccds_status\tcds_strand\tcds_from\tcds_to\tmatch_type\tstart\tstop"; while (<>) { next if /^#/; chomp; my ($chromosome, $nc_accession, $gene, $gene_id, $ccds_id, $ccds_status, $cds_strand, $cds_from, $cds_to, $cds_locations, $match_type) = split /\t/; my @cds_ranges = map { my ($start, $stop) = split /-/; {start => $start+1, stop => $stop+1} } (split /,\s/, (substr $cds_locations, 1, -1)); foreach (@cds_ranges) { say "$chromosome\t$nc_accession\t$gene\t$gene_id\t$ccds_id\t$ccds_status\t$cds_strand\t$cds_from\t$cds_to\t$match_type\t$_->{start}\t$_->{stop}"; } }
そして、tsvファイルを吐きます。
cat CCDS.current.txt | perl parse_ccds.pl > CCDS.current.tsv
あとはRで読み込むだけです。幸い、makeGRangesFromDataFrameを使うことで、読み込んだデータフレームを簡単にGRangesオブジェクトに変換できます。
> library(GenomicRanges)
> ccds = read.csv("CCDS.current.tsv", sep="\t")
> ccds = makeGRangesFromDataFrame(ccds, strand.field="cds_strand", keep.extra.columns=T)
> ccds
GRanges object with 287116 ranges and 8 metadata columns:
seqnames ranges strand | nc_accession gene gene_id ccds_id ccds_status cds_from cds_to
<Rle> <IRanges> <Rle> | <factor> <factor> <integer> <factor> <factor> <integer> <integer>
[1] 1 [801942, 802433] - | NC_000001.8 LINC00115 79854 CCDS1.1 Withdrawn 801942 802433
[2] 1 [861321, 861392] + | NC_000001.10 SAMD11 148398 CCDS2.2 Public 861321 879532
[3] 1 [865534, 865715] + | NC_000001.10 SAMD11 148398 CCDS2.2 Public 861321 879532
[4] 1 [866418, 866468] + | NC_000001.10 SAMD11 148398 CCDS2.2 Public 861321 879532
[5] 1 [871151, 871275] + | NC_000001.10 SAMD11 148398 CCDS2.2 Public 861321 879532
... ... ... ... ... ... ... ... ... ... ... ...
[287112] Y [16168463, 16168738] + | NC_000024.9 VCY1B 353513 CCDS56618.1 Public 16168169 16168738
[287113] Y [16835028, 16835148] + | NC_000024.9 NLGN4Y 22829 CCDS56619.1 Public 16835028 16953141
[287114] Y [16936067, 16936252] + | NC_000024.9 NLGN4Y 22829 CCDS56619.1 Public 16835028 16953141
[287115] Y [16941609, 16942398] + | NC_000024.9 NLGN4Y 22829 CCDS56619.1 Public 16835028 16953141
[287116] Y [16952292, 16953141] + | NC_000024.9 NLGN4Y 22829 CCDS56619.1 Public 16835028 16953141
match_type
<factor>
[1] Identical
[2] Identical
[3] Identical
[4] Identical
[5] Identical
... ...
[287112] Identical
[287113] Identical
[287114] Identical
[287115] Identical
[287116] Identical
-------
seqinfo: 24 sequences from an unspecified genome; no seqlengths
遺伝子毎のCDS領域を得たければ、Gene IDのアノテーションがccdsにありますのでsplitでGRangesListオブジェクトにできます。
> ccdsbygene = split(ccds.public, ccds.public$gene_id)
> ccdsbygene
GRangesList object of length 18383:
$1
GRanges object with 8 ranges and 8 metadata columns:
seqnames ranges strand | nc_accession gene gene_id ccds_id ccds_status cds_from cds_to
<Rle> <IRanges> <Rle> | <factor> <factor> <integer> <factor> <factor> <integer> <integer>
[1] 19 [58858387, 58858394] - | NC_000019.9 A1BG 1 CCDS12976.1 Public 58858387 58864802
[2] 19 [58858718, 58859005] - | NC_000019.9 A1BG 1 CCDS12976.1 Public 58858387 58864802
[3] 19 [58861735, 58862016] - | NC_000019.9 A1BG 1 CCDS12976.1 Public 58858387 58864802
[4] 19 [58862756, 58863052] - | NC_000019.9 A1BG 1 CCDS12976.1 Public 58858387 58864802
[5] 19 [58863648, 58863920] - | NC_000019.9 A1BG 1 CCDS12976.1 Public 58858387 58864802
[6] 19 [58864293, 58864562] - | NC_000019.9 A1BG 1 CCDS12976.1 Public 58858387 58864802
[7] 19 [58864657, 58864692] - | NC_000019.9 A1BG 1 CCDS12976.1 Public 58858387 58864802
[8] 19 [58864769, 58864802] - | NC_000019.9 A1BG 1 CCDS12976.1 Public 58858387 58864802
match_type
<factor>
[1] Identical
[2] Identical
[3] Identical
[4] Identical
[5] Identical
[6] Identical
[7] Identical
[8] Identical
...
<18382 more elements>
-------
seqinfo: 24 sequences from an unspecified genome; no seqlengths
RStudioから使えるR Markdownで日本語のPDFを作成する
RStudioではRの実行結果を埋め込んでドキュメントやレポートを作成するのにLaTeXベースのものとMarkdown(の亜種のR Markdown)ベースのものが選べるのですが、R Markdownから日本語を含んだレポートをPDF化するのに手間取ったのでその方法を記述します。 私の環境はMac OS X Mavericksですが、他の環境でも共通する点はあると思います。
準備
まずMarkdownをLaTeXにするPandocとLaTeXの処理系(XeLaTeX)が必要です。
恐らくPandocに関してはRStudioに同封されているので問題無いと思います。
LaTeXの処理系に関してはMacTeXを利用するのが最も手早いです。しかし、2GB以上もある巨大なファイルをダウンロードし、インストールする羽目になるので、コアだけを集めたBasicTeXでも十分です。
BasicTeXをインストールしたら、いくつか必要なパッケージをインストールします1:
sudo tlmgr update --self --all sudo tlmgr install framed sudo tlmgr install titling sudo tlmgr install collection-langjapanese
設定
続いてRStudioを起動し、必要な設定を行います。
“RStudio” ➠ “Preferences…” ➠ “Sweave”と選択し、取り敢えずデフォルトをの設定を以下のようにします。
- Weave Rnw files using: knitr
- Typeset LaTeX into PDF using: XeLaTeX
- Preview PDF after compile using: System Viewer
最後の設定はRStudioのPDFビューアーが気に入らなかっただけで、私の好みです(PDF Reader Xを使っています)。
また、この設定はデフォルトのもので、プロジェクトでは別個に設定する必要があります。
R Markdownファイルの作成
“File” ➠ “New File” ➠ “R Markdown…” と選択し、適当にタイトルを設定してR Mawkdownのひな形作成します。ここでは、タイトルは「日本語PDF作成」としました。
しかし、これだけではプロジェクトの設定が反映され、先ほどのデフォルト設定が無視されてしまう場合があります。 そこで、ファイルの先頭を以下のように変更します。
変更前:
--- title: “日本語PDF作成” output: pdf_document ---
変更後:
---
title: “日本語PDF作成"
output:
pdf_document:
latex_engine: xelatex
mainfont: Hiragino Kaku Gothic Pro
---
latex_engineとしてxelatexを設定することで、UTF8などのエンコーディングに対応できるようになり、mainfontを設定することで、出来上がったPDFファイルの日本語フォントが表示されなくなる問題を修正しています。
ここでは、Hiragino Kaku Gothic Proを使いましたが、MacのFont Bookから、日本語に対応した好きなフォントを選ぶこともできます。
最後に日本語を含むR Markdownを書き、”Knit PDF”をクリックするとPDFが表示されます。
参考
JuliaTokyo #2 でモンテカルロ法について話してきました
9月27日の JuliaTokyo #2 で『Juliaで学ぶ Hamiltonian Monte Carlo (NUTS 入り)』と題してHamiltonian Monte Carlo (HMC) とStanで使われているNUTSについてお話してきました。
内容
サンプルコード: bicycle1885/JuliaTokyo2HMC · GitHub
MCMC自体については30分ではほとんど説明できませんでしたが、会場のほとんどの方が既に知っているようでしたのでちょうど良かったようです。
全体の流れとしては、Metropolis-Hastings、Hamiltonian Monte Carlo、No-U-Turn Samplerの順に3つサンプラーを実装コードと正規分布からのサンプリングを例に紹介し、どのように前のサンプラーの弱点を克服していってるかを説明しました。
Unicodeプログラミング
サンプルコードはJuliaの愉しみのひとつであるUnicode文字をふんだんに使っていたところが面白いところです。
function build_tree(L::Function, ∇L::Function, θ::Vector{Float64}, r::Vector{Float64}, u::Float64, v::Int, j::Int, ϵ::Float64)
if j == 0
θ′, r′ = leapfrog(∇L, θ, r, v * ϵ)
C′ = u ≤ exp(L(θ′) - r′ ⋅ r′ / 2) ? Set([(θ′, r′)]) : Set([])
s′ = int(L(θ′) - r′ ⋅ r′ / 2 > log(u) - Δmax)
return θ′, r′, θ′, r′, C′, s′
else
θ⁻, r⁻, θ⁺, r⁺, C′, s′ = build_tree(L, ∇L, θ, r, u, v, j - 1, ϵ)
if v == -1
θ⁻, r⁻, _, _, C″, s″ = build_tree(L, ∇L, θ⁻, r⁻, u, v, j - 1, ϵ)
else
_, _, θ⁺, r⁺, C″, s″ = build_tree(L, ∇L, θ⁺, r⁺, u, v, j - 1, ϵ)
end
s′ = s′ * s″ * ((θ⁺ - θ⁻) ⋅ r⁻ ≥ 0) * ((θ⁺ - θ⁻) ⋅ r⁺ ≥ 0)
C′ = C′ ∪ C″
return θ⁻, r⁻, θ⁺, r⁺, C′, s′
end
end
ちなみにjulia-vimはlatex記法サポートしてます https://t.co/iT48OQLi6r #JuliaTokyo
— 佐藤ヾ(⌒(｢ 'ω')｢ｶﾞｵｰ建太 (@bicycle1885) 2014, 9月 27箴言
Knuth 「尚早な最適化は諸悪の根源」
M. Jackson 「最適化第一法則: 最適化するな, 第二法則: まだするな」
#JuliaTokyo
— 佐藤ヾ(⌒(｢ 'ω')｢ｶﾞｵｰ建太 (@bicycle1885) 2014, 9月 27+α
俺のRandomForests.jlが某巨大企業で使われているらしい！凄い！ #JuliaTokyo
— 佐藤ヾ(⌒(｢ 'ω')｢ｶﾞｵｰ建太 (@bicycle1885) 2014, 9月 27社内で改造しているようですが、手続き上のアレでプルリクは受けられないようです(´・ω・`)
コマンドを並列に実行するGNU parallelがとても便利
最近のコンピュータは複数のCPUコアを持っているので並列にコマンドを実行することができます。 たくさんの同じようなファイルに同じ処理を実行することは、私のやっているバイオインフォマティクスではよくあります。
しかし自分で並列に実行するスクリプトを書くことはそれほど簡単ではなく、ログや実行結果の確認など煩雑な処理を書かなければいけません。 この記事では、そうした処理を簡単にするGNU parallelというツールを紹介します。
GNU parallel
UNIX系のOSではインストールはとても簡単です。MacでしたらHomebrewを使って、Linuxでは各ディストリビューションのパッケージマネージャからインストールできます。 詳しくはGNU parallelのウェブページを参照して下さい(http://www.gnu.org/software/parallel/)。
Homebrew:
brew install parallel
パッケージマネージャが使えない環境でも、ソースコードからビルドすることができます。 ダウンロードページ(http://ftp.gnu.org/gnu/parallel/)から、最新のparallel(parallel-latest.tar.bz2)をダウンロードして展開し、ビルドしましょう。 必要ならばconfigureに--prefixも渡せます。
Source:
./configure && make && make install
なお、最初にparallelコマンドを実行する際に引用のお願いメッセージが出ますが、これはparallel --bibtexを一度実行すれば抑えられます。
Hello, GNU parallel !
インストールが終わったら、早速GNU parallelを実行してみましょう。
~/t/sample $ parallel echo ::: hello world ! hello ! world
GNU parallelのparallelコマンドは、実行するコマンドと、そのコマンドに適用するいくつかの引数を引数としてとります。
上の例では、echoが実行するコマンドで、:::の後に続く3つの引数hello, world, !がそれぞれechoコマンドに適用する引数です。
echoされたメッセージが元の引数の順番と一致していないことから分かる通り、echoコマンドは並列に実行されるため、各引数の実行順は不定になります。
:::は複数つなげることで、引数の組み合わせの積をつくることもできます。
~/t/sample $ parallel echo ::: hello, bye, ::: Alice Bob Charlie ::: ! hello, Alice ! hello, Bob ! hello, Charlie ! bye, Alice ! bye, Bob ! bye, Charlie !
ファイルや標準入出力から引数を与える
以下のようなファイル名のリストを収めたファイル(list.txt)があるとします。
~/t/sample $ cat list.txt foo.txt bar.txt baz.txt
-a <file>オプションで、ファイルの各行を引数としてコマンドを並列に走らせることができます。
list.txtの各ファイルを並列に圧縮するには以下のようにします。
~/t/sample $ ls bar.txt baz.txt foo.txt list.txt ~/t/sample $ parallel -a list.txt gzip ~/t/sample $ ls bar.txt.gz baz.txt.gz foo.txt.gz list.txt
パイプで引数を渡すこともできます。その時は-a -とオプションを渡します。
~/t/sample $ cat list.txt | parallel -a - gzip ~/t/sample $ ls bar.txt.gz baz.txt.gz foo.txt.gz list.txt
引数の区切り文字はデフォルトでは\n(newline)です。
-0(--null)オプションで区切り文字が\0(null)文字になります。
他にも、<や::::を使って引数の納められたファイルを渡すことができます。
~/t/sample $ parallel gzip < list.txt ~/t/sample $ parallel gzip :::: list.txt
引数の置き換え
コマンドに渡す引数の位置が最後でない場合や、特別な処理をしたい場合に引数の置換え位置を指定するをする必要があります。
その際には{}を使います。
~/t/sample $ parallel 'find {} -name "README*"' ::: ~/vendor/julia ~/vendor/vector
/Users/kenta/vendor/vector//old-testsuite/microsuite/README
/Users/kenta/vendor/vector//README.md
/Users/kenta/vendor/julia//contrib/mac/app/README
/Users/kenta/vendor/julia//contrib/README.ackrc.txt
...
さらに、{.}でファイル名の拡張子を除いた引数にしたり、{/}でファイル名だけ取り出したりできます。
~/t/sample $ parallel 'echo {.}' ::: tmp/foo.txt.gz tmp/bar.txt.gz
tmp/foo.txt
tmp/bar.txt
~/t/sample $ parallel 'echo {/}' ::: tmp/foo.txt.gz tmp/bar.txt.gz
foo.txt.gz
bar.txt.gz
実行されるコマンドの確認
--dry-runをつかえば、引数がどのようにコマンドに適用されてどのようなコマンドが実際に走るかを確認できます。
実際に実行する前に一度確認すると良いでしょう。
~/t/sample $ cat list.txt | parallel -a - --dry-run gzip gzip foo.txt gzip bar.txt gzip baz.txt
実行結果の取得
標準出力と標準エラーに吐かれた結果をファイルとして取得することも容易です。
--results <outputdir>とすることで出力を各引数毎にディレクトリに構造化して保存できます。
~/t/sample $ parallel --results results 'perl -E "say STDOUT \"stdout\"; say STDERR \"stderr\""' ::: A B C
stdout
stderr
stdout
stderr
stdout
stderr
~/t/sample $ tree results/
results/
└── 1
├── A
│ ├── stderr
│ └── stdout
├── B
│ ├── stderr
│ └── stdout
└── C
├── stderr
└── stdout
4 directories, 6 files
~/t/sample $ cat results/1/A/stderr
stderr
~/t/sample $ cat results/1/A/stdout
stdout
実行結果の確認
時間のかかるコマンドや引数が多い場合など、実行結果の確認が大変な場合があります。
そこで、--joblog <logfile>を使えばコマンドのexit statusが確認しやすくなります。
~/t/sample $ parallel --joblog joblog.txt 'sleep {}; exit {}' ::: 0 1 2 3
~/t/sample $ cat joblog.txt
Seq Host Starttime JobRuntime Send Receive Exitval Signal Command
1 : 1407646480.280 0.072 0 0 0 0 sleep 0; exit 0
2 : 1407646480.284 1.144 0 0 1 0 sleep 1; exit 1
3 : 1407646480.289 2.125 0 0 2 0 sleep 2; exit 2
4 : 1407646480.294 3.153 0 0 3 0 sleep 3; exit 3
一目見て分かる通り、実行時間や実際に実行されたコマンドなども記録されるため、問題が発生した時のトラブルシューティングに非常に役に立ちます。 また、終わったコマンドから順に追記されていくため、進捗状況の確認もできます。
並列ジョブの制御
ファイルのダウンロードなどを並列に行いたいが過剰に負荷を掛けたくない場合や、逆にCPU boundでない処理をCPUのコア数より多く並列させたいときなどに同時実行するジョブ数を指定できると便利です。
デフォルトでは、マシンのコア数分だけ並列させるため、共有サーバーなどではCPUを占拠してしまい迷惑になることも考えられます。
そのようなときは--jobs <N>(-j <N>)オプションを使って、
以下のように6つの引数を1, 2, 3, 6並列でそれぞれ実行して実行時間を見てみましょう。
~/t/sample $ time parallel --jobs 1 'sleep {}' ::: 1 1 1 1 1 1
7.06 real 0.26 user 0.10 sys
~/t/sample $ time parallel --jobs 2 'sleep {}' ::: 1 1 1 1 1 1
3.58 real 0.27 user 0.11 sys
~/t/sample $ time parallel --jobs 3 'sleep {}' ::: 1 1 1 1 1 1
2.44 real 0.28 user 0.11 sys
~/t/sample $ time parallel --jobs 6 'sleep {}' ::: 1 1 1 1 1 1
1.30 real 0.30 user 0.12 sys
--jobs 50%などと書くことで全体の半分のCPUコアを使う指定もできます。
~/t/sample $ time parallel --jobs 50% 'sleep {}' ::: 1 1 1 1 1 1
3.59 real 0.26 user 0.11 sys
また、メモリを沢山必要とする計算などでは、--noswapオプションを指定することで、メモリのスワップが発生している時には新しいジョブを実行しないようにすることもできます。
他の便利そうな機能
ここで紹介した機能はGNU parallelのごく一部で、他にも有用なオプションなどがたくさんあります。
一度、man parallel_tutorialやman parallelに目を通すことをオススメします。
私はまだ実際に実行して確認していませんが、GNU parallelにはSSH越しに複数のノードで並列にジョブを実行させる機能やパイプライン処理で中間のボトルネックになっている処理だけ並列化するなどの機能もあるようです。
これらの機能は確認し次第、ブログに追記していこうと思います。
statisticsにt検定を実装しました
ちょっとHaskellのstatisticsパッケージを見ていたら、統計的仮説検定に何故かt検定が無かったのでチョイチョイと実装しました。
https://github.com/bos/statistics/pull/66
masterの一歩手前のconfidenceブランチに取り込まれましたので、そのうち使えるようになると思います。
この実装で以下の3つの2標本検定が使えるようになります。
- Student's t-test
- Welch's t-test
- paired samples t-test
詳細については、Wikipediaを参照下さい。
statisticsパッケージの検定のAPIがに統一感がなくてイケてない感じがするので、なにか思いついたらそちらも修正したいです。あとp値や統計量が得られないのはさすがに良くない気がします。
JuliaTokyo #1 を開催しました
日本で(多分)最初のJuliaの勉強会 JuliaTokyo #1 を開催しました。イベントページはこちらです。 http://juliatokyo.connpass.com/event/6891/
当日は40人ほどの参加者が集まり、いつの間にJuliaがそんなに一般化したんだという感じでしたが、Julia歴1週間とか当日の朝から始めましたというように、何処からかJuliaが熱いらしいとの噂を聞きつけてJuliaを使ってみようという新しいもの好きの人たちが多い印象でした。
私はパッケージの作り方についてお話をしてきた次第です。
以下は、それぞれの発表と私のひとこと感想です。
メインセッション
Juliaのこれまでとこれから - @sorami
Juliaの背景的な情報やどんなエコシステムになっているかを紹介していました。 JuliaはPythonやRからユーザを取り込みつつ、共に進化していく道を辿るようです。
Julia100本ノック - @chezou
Juliaでnumpyの100 numpy exercisesを書いた話です。 Juliaのやり方とnumpyのやり方はベクトル化などの使い方がかなり違い、苦労があったようです。 あとJuliaのベクトル内包表記は速いのでガンガン使っていきましょう！
Juliaのパッケージを作ろう! - @bicycle1885
私の発表です。 Juliaのパッケージをどのようにつくるかを解説しました。 パッケージ名やディレクトリ構成などについて、ドキュメントに書いてない慣習などが大事だったりします。
ビジュアライズJulia - @nezuq
ビジュアライズJulia from nezuQ
"The Grammar of Graphics"のlayerやscaleの概念や、Gadfly.jlとの関係の話でした。 可視化関係のライブラリはまだそれほどなので、今後充実してくると良いですね。
メカ女子将棋 - @kimrin
Juliaで将棋AIを作った話でした。やっぱりJuliaはCなどに比べ開発がずっとしやすいですね。 盤面の表現にはJuliaの128bit整数を使っているようです。
A Brief JuliaCon Report / Natural Language Processing with Julia - Pontus Stenetorpさん
https://github.com/JuliaCon/presentations/blob/master/JuliaNLP/JuliaNLP.pdf
先月Chicagoで開催されたJuliaConのレポートと、そのときに発表されたJuliaでの自然言語処理のお話をしていただきました。 JuliaConでは、アカデミックの人がメインだったものの、実際にサービスにJuliaを投入している会社もあったようです。 Juliaで自然言語のパーサallenを作って、パフォーマンスも満足いくものになったようです。
LT
Juliaでちゃんとしたコードを書く - 関さん
当日の朝からJuliaを始めて、どうするとJuliaが遅くなるかのお話でした。 マニュアルのperformance tipsマジ重要です。読みましょう!!
Julia0.3でランダムフォレスト - @gepuro
DecisionTree.jlを使ってみたようです。 DecisionTree.jlはID3を決定木構築のアルゴリズムに使っているので、CARTが使いたければ私のRandomForests.jlを使いましょう！
Plotly Julia API - Yasunobu Igarashiさん
Plotlyというデータの可視化や共有のためのプラットフォームの紹介でした。 JuliaのAPIもあるんですね、すごい！
LightTableを使いましょう - @QuantixResearch
LightTable + Jewel でJuliaを開発しようという話でした。 インラインで変数の値が見れたり、補完もかなり強力なようです。あとで使ってみよう！
感想
新しいおもちゃを手にしたときの子供のような、Juliaに対する期待と興奮が垣間見れたように思います。 おぉ！Juliaすごい！いいじゃん！というような熱気がありました。 Twitterでしか見たことない人と実際にお話したり、情報交換をしたりも大変有益でした。 反省点としては、ちょっと初心者向けの話が少なすぎたかなと思います。次回は初心者向けセッションを設けて欲しいとの要望が多かったです。 会場を用意していただいた株式会社ブレインパッドの皆様、本当にありがとうございました。
あと、本家Juliaステッカーをもらえたのがすごい嬉しかった！
本場シカゴからやって来たステッカー貰った(((o(*ﾟ▽ﾟ*)o))) pic.twitter.com/1AJeQCp0ml
— 佐藤ヾ(⌒(｢ 'ω')｢ｶﾞｵｰ建太 (@bicycle1885) 2014, 7月 5