二項分布まわりのまとめ１「ポアソン分布の導出」

このまま負の二項分布、ベータ二項分布まではまとめる予定。

二項分布は以下の式
{ f(k | n, p) = {}_{n}C{}_{k} p^k (1-p)^{n-k} }

nが十分に大きく、pが十分に小さい場合を考える。
ここで np = \lambda とすることで
ポアソン分布
{ f(k | \lambda) = \frac{\lambda^k}{k!} e^{- \lambda}}
を導出する。

ネイピア数があるのでたぶん{ e^{- \lambda} = \lim_{n \to \infty} (1 + \frac{ - \lambda}{n})^{n} }こんな感じなんでしょう。

k!やλを早めに集めて
{ f(k | n, p) = \frac{n^{k} p^{k}}{k!} \frac{n!}{(n-k)! n^{k}} (1-p)^{n-k} }
{ right side = \frac{n^{k} p^{k}}{k!} \frac{n(n-1)..(n-k+1)}{n^{k}} (1-p)^{n-k} }
{ right side = \frac{n^{k} p^{k}}{k!} 1 (1 - \frac{1}{n}) (1- \frac{2}{n})...(1 - \frac{(n-k+1)}{n}) (1-p)^{n-k} }

np = λより
{ right side = \frac{\lambda^{k}}{k!} 1 (1 - \frac{1}{n}) (1- \frac{2}{n})...(1 - \frac{(n-k+1)}{n}) (1-\frac{\lambda}{n})^{n-k} }
{ right side = \frac{\lambda^{k}}{k!} 1 (1 - \frac{1}{n}) (1- \frac{2}{n})...(1 - \frac{(n-k+1)}{n}) \frac{(1-\frac{\lambda}{n})^{n}}{(1-\frac{\lambda}{n})^{k}} }

{ n \to \infty}の極限をとると{ \frac{\lambda^{k}}{k!} (1-\frac{\lambda}{n})^{n}}が残る。

あとは、上述の通り。本当によく忘れる。

Rで日付->数値の変換

時系列のデータを扱うとき、
"7/1/2018 3:33:00 AM"みたいなデータから、
"1530383580"のように一つの数字に変換したい時がある。

そのためには、Rならas.POSIXct()を使えば良い。
以下参照。
https://stackoverflow.com/questions/8215404/change-from-date-and-hour-format-to-numeric-format

>as.numeric(as.POSIXct("7/1/2018  3:33:00 AM", format="%m/%d/%Y  %H:%M:%S"))
[1] 1530383580

ちなみに、あちこちに情報があるが、
タイムゾーンに関する警告も出る。以下のようにすれば良い。

>Sys.setenv("TZ" = "Asia/Tokyo")

以下参照
http://uribo.hatenablog.com/entry/2017/11/30/194734

管理者権限なしでR-3.4.0 + autoconf-2.6.9をインストール

基本的には以下参照。

http://labo.utsubo.tokyo/2017/01/30/r%E3%82%92%E3%82%BD%E3%83%BC%E3%82%B9%E3%82%B3%E3%83%BC%E3%83%89%E3%81%8B%E3%82%89%E3%82%A4%E3%83%B3%E3%82%B9%E3%83%88%E3%83%BC%E3%83%AB/

これで大丈夫と見せかけて、Bioconductorなどを
利用しようとすると、htslibなどのインストール時に
今度は「autoconfが古過ぎる」と言われてしまう。

autoconfを入れるには、m4-1.4.18のインストールが必要。
こちらはもう少し簡単だった。以下参照。
http://tetsuyai.hatenablog.com/entry/20110216/1297939037

CX_reportからCGmapへの変換

みんな大好きかつ、引くほど遅いmethylation extractorが素敵なBismark。
https://www.bioinformatics.babraham.ac.uk/projects/bismark/

(遅いのがいやならmethylpyか自分で書けばよい
https://github.com/yupenghe/methylpy)

CX_report.txtと名前のついた出力を吐くので
以降の解析にはこのファイルを使えばよい。

しかし最近では、BSSeeker2の影響か、
https://github.com/BSSeeker/BSseeker2
CGmapという形式が若干広まってきている。
MethGoとかHOMEとか。
http://methgo.readthedocs.io/en/latest/
https://github.com/ListerLab/HOME

そのため、以下のようなワンライナーで書き換える。

awk -F "\t" -v OFS="\t" '($4+$5>0){print $1,($3=="+"?"C":"G"),$2,$6,substr($7,1,2),$4/($4+$5),$4,$4+$5}'

Bisulfighter用には、また後日書く予定。
https://epigenome.cbrc.jp/bisulfighter

SimpleChangepointCalculatorについては
http://itolab.med.kyushu-u.ac.jp/CPT/
BMapの出力がバイナリっぽいので、
ソースコードを読まないといけない。

JBrowseでトラックを追加する

JBrowse · A fast, embeddable genome browser built with HTML5 and JavaScriptはApacheベースで割と簡単にインストールでき、比較的サクサク動くので使いやすい。

が、インストールするときにデータを置く位置をデフォルトから変えると、結構面倒臭い。
さらに、頻繁にアップデートしているのはありがたいのだけど、データの追加方法を調べてみると、見たページによって違う。

結果としては、trackList.jsonにデータ情報を書き込めば良い。
tracks.confに書く方法もあるようだが(誤読かも)、どうにもうまくいかなかった。

trackList.jsonの場所はどこでもいい。
とりあえず、"urlTemplate"以外の部分をしっかりと入力するのが大事。
そうすると、入力したファイルへのパスについてのエラーが出てきて、"urlTemplate"をどう書いたらいいかが分かる。
何も反映されない、表示もされない間は、それ以外の部分が間違っている。
具体的には以下のような感じ。

echo ' {
         "autoscale": "local",
         "logScaleOption": true,
         "label" : "mylabel",
         "key" : "mykey",
         "storeClass" : "JBrowse/Store/SeqFeature/BigWig",
         "urlTemplate" : "../../mydata.bw",
         "type" : "JBrowse/View/Track/Wiggle/XYPlot",
     } ' | bin/add-track-json.pl data/json/MyGenome/trackList.json