この記事は
今年読んだ一番好きな論文2017 Advent Calendar 2017
の2日目が空いているということにこの記事を書いてから気づいて、1000円相当の参加賞があることに目がくらんで加筆修正した記事です。2017年で一番好きかというとそうでm(ここで文章が途絶えている けれどもこういう状況でいい検定法がなくて困る人がいるかもと思って書いておく。
昨年度は猫じゃら賞を受賞し、ちまたで噂のがくしん本をいただきました(白目
今年読んだ一番好きな論文2016 - 驚異のアニヲタ社会復帰への道
読んだ。
Cell. 2016 Dec 1;167(6):1495-1510.e12.
腸内フローラの細菌数はoscillation していて、そのoscillation は宿主の遺伝子発現、メチル化状態に影響していて、腸管を超えて肝臓の発現oscillation にも、肝障害応答にも影響している、という話。メタゲノムと概日リズムすげぇ、2017年のノーベル生理学・医学賞はこれだったしね、という話。
このなかで、oscillation という現象が起きているのか起きていないのかを検定している。例えば、図3B を引用すると、消化管にいる微生物数のoscillation の有無が宿主の応答(トランスクリプトームやエピゲノム)に対してどのように影響しているかを確認するため、広域抗菌薬(とにかくバイキンは殺すやつ、具体的にはバンコマイシン、アンピシリン、カナマイシン、メトロニダゾールで、人に使ってもほとんどの菌が死ぬ)で腸内フローラを破壊している。
青は対照で、赤が抗菌薬である。対照は普通にoscillation している様子は見て分かるし、抗菌薬で細菌が常に死滅させられているとoscillation していない様子も見て分かる。
これが統計学的にoscillation しているかしていないのか(帰無仮説検定だけど)でいうと、対照は有意差があるのでoscillation していそうで、抗菌薬は棄却されないのでoscillation していなさそう、ということになる。
この論文で、oscillation の検定に、JTK_cycle という検定を使っている。oscillation という現象は時系列データであるし、時系列解析はかなり難しいので「各時刻でt検定しました()」とか「とりあえずANOVA しました()」みたいなことになるのだが、これらの検定は各データ点が独立という前提の検定なので、過去のデータに影響を受ける時系列データでは不適当なのであるが、かと言ってじゃあいい検定手法あるの? と聞かれるとなるほど、わからん、ということになる。
ということでoscillation に使える帰無仮説検定としてこの論文を紹介する。
JTK_CYCLE: An Efficient Nonparametric Algorithm for Detecting Rhythmic Components in Genome-Scale Data Sets
J Biol Rhythms. 2010 Oct;25(5):372-80.
HughesLab:JTK Cycle - OpenWetWare
R パッケージについてはMetaCycle
Bioinformatics. 2016 Nov 1;32(21):3351-3353.
oscillation の検定にはいろいろあるが、5種類試してみて、時系列のデータ間隔と全体の長さが大事だが、JTK_cycle
について言えば時系列のデータ間隔がかなり重要、と言っている。
cycSimu4h2d のデータセットはここからシミュレーションで作っている。
J Biol Rhythms. 2014 Aug;29(4):231-42.
ちょっと詳しい話はこちら。
PLoS Comput Biol. 2015 Mar 20;11(3):e1004094.
JTK は、ヨンクヒール・タピストラ・ケンダール検定の頭文字である。
上昇トレンドの検定 ヨンクヒール・タプストラ検定 - 驚異のアニヲタ社会復帰への道
ケンダール検定は、データのペアワイズの差分符号を調べて検定するが、これにヨンクヒール・タピストラ検定という、データに順序があるときの制約をいれて統計量を算出するノンパラメトリックな検定の合わせ技、ということになっている。
いま、時系列データ,
に対して
という統計量を計算する。これは組み合わせになるのでめちゃめちゃ計算量が増えるが、がんばれば正規分布に近似できるらしく、なんかいろいろ計算量を減らす工夫があるらしい。
また、JTK_cycle 自体には、最適なoscillation 周期を(正確確率検定的に) 最小のp値を与える周期として推定する。また、周期性があるとしてベースラインからの振幅も推定するようになっているらしい。
基本的には、手持ちのデータから計算できる符号のみに興味があるノンパラメトリックな手法なので、ほとんどの状況で使える。ノンパラメトリックなので検出力的には緩い場合もあるかもしれないが、その他の性能評価論文ではデータ間隔を密にするようにすればどうにかなるようなので、oscillation を検定する手法を知らないよりかははるかにマシ。
というわけでやってみる。cycSimu4h2d は、適当なcos 関数にピークを入れてみたり、矩形にしてみたり、定常性を保ったり失わせたり(傾き trend がつく) したりしてシミュレーションしている。
ここで、negt はたぶんoscillation がないシミュレーションで、そのほかは三角関数でoscillation があるようにしている。
oscillation がない緑は、JTK_cycle 検定で有意にならなかった(破線)ものが多く、oscillation があるものとしてシミュレーションしたデータは多くが有意になっている。
ここで、検定で帰ってくるp値は、遺伝子発現など多くの検定を期待しているので、(おそらく、論文中には書いてあるがソースコードはたどっていない)ボンフェローニ補正と、BH法のq値が返ってくる。最初に挙げたCell 論文ではトランスクリプトームの発現パターンが、対照と抗菌薬でどちらでもoscillation している遺伝子群(shared)、対照ではoscillation しているが抗菌薬で消失した遺伝子群(lost)、対照ではoscillation していないが抗菌薬でoscillation するようになった遺伝子群(gained) をそれぞれp値 0.05 q値 0.1 で取り出していて、それぞれ
shared: circadian
lost: DNA replication, nucleotide
gained: focal adhesion, pyruvate matabolism
というようなパスウェイがKEGG で出てくる、というような解析をしている。
curvID 0 4 8 12 16 20 24 28 32 36 40 44 1 cos2_1756 2.11 -1.24 -1.15 -2.04 0.31 0.24 1.40 0.77 0.80 -1.42 0.51 0.41 2 cos2_1537 0.64 -0.52 0.50 -0.08 -0.90 -0.57 0.93 2.09 0.02 -1.31 -2.75 -0.31 3 negt_558 -1.01 0.50 0.95 -0.70 -0.33 -0.59 0.57 0.32 1.45 -0.12 1.05 0.77 4 negt_1057 -0.05 -1.15 -0.19 0.52 0.44 -0.23 -2.56 0.77 0.41 0.45 1.58 -1.25 5 cosp_3926 2.25 1.93 0.97 -1.87 -2.72 -0.27 1.36 2.51 0.28 -1.36 -0.58 1.35 6 cosp_3960 -1.86 -0.76 0.24 3.58 1.07 -2.45 -3.50 -1.46 0.00 1.54 0.79 -0.93 7 cosp_2183 1.06 1.95 -0.18 -0.83 -1.04 -0.46 -0.13 2.39 -0.05 -1.35 -1.74 -1.28 8 negt_2142 -0.56 -0.34 -0.12 -0.11 2.32 -1.20 0.48 -1.57 0.53 -0.81 -0.18 -1.36 9 pekp_4656 -1.02 -0.81 -2.84 -1.42 1.39 2.11 -1.50 -2.08 -1.95 -1.99 0.45 1.37 10 pekp_4740 -2.67 -1.14 -2.18 -0.99 0.65 0.05 -2.36 -1.92 -2.58 -1.22 3.82 0.09 11 negt_5430 0.63 0.25 0.09 -0.90 0.44 -0.67 -0.21 -0.02 -0.14 0.55 1.14 0.33 12 negt_5514 0.18 -0.93 0.26 -1.30 -1.06 1.77 -0.82 0.90 -0.03 0.79 0.84 1.17 13 negt_6004 0.65 1.12 1.30 -0.86 0.38 1.25 -0.08 -0.29 1.23 -0.11 -0.75 0.22 14 squa_7485 2.05 -2.71 -2.48 -1.38 0.76 1.64 1.12 -1.77 -2.71 -4.13 3.46 2.21 15 negt_6385 -0.78 1.80 0.05 0.93 0.83 1.35 2.28 -0.89 -0.20 -1.03 -1.13 0.08 16 squa_6903 1.49 -4.39 -3.02 -1.05 2.65 1.01 0.94 -1.87 -4.12 -1.76 2.89 2.31 17 negt_8645 0.29 0.36 -0.39 0.44 1.12 0.80 0.78 0.59 0.63 -0.24 1.89 -0.50 18 negt_8219 0.09 -0.70 0.02 1.95 1.52 0.68 -1.67 -1.59 -0.96 0.18 0.89 -1.56 19 negt_8578 -1.68 -0.45 0.91 -0.87 1.12 -0.09 0.21 -0.44 0.86 0.52 0.12 -0.85 20 negt_9637 1.10 -0.06 -0.34 0.15 1.56 0.27 -0.58 1.01 0.04 0.15 -1.80 -0.97
library(MetaCycle) library(stringr) alpha <- 0.05 write.csv(cycSimu4h2d, file="cycSimu4h2d.csv", row.names=FALSE) # output integration file with analysis results from each method meta2d(infile="cycSimu4h2d.csv", filestyle="csv", outdir="example", timepoints="Line1",cycMethod="JTK") leg <- str_extract(cycSimu4h2d[,1], "^\\w{4}") p <- read.csv("./example/JTKresult_cycSimu4h2d.csv") cols <- 1:5 matplot(as.numeric(colnames(cycSimu4h2d)[-1]), t(cycSimu4h2d[,-1]), las=1, col=cols, type="o", pch=15, xlab="Time", ylab="Amplitude", lwd=ifelse(p$ADJ.P<alpha, 3, 2), lty=ifelse(p$ADJ.P<alpha, 1, 3)) legend("bottomright", legend=levels(factor(leg)), col=cols, lty=1, lwd=3) legend("topleft", legend=sprintf("p %s %.2f", c("<", ">"), alpha), lty=c(1,3), lwd=3)
CycID BH.Q ADJ.P PER LAG AMP 1 cos2_1756 0.49574766 0.247873828 28 26 1.1766793 2 cos2_1537 0.09273990 0.032458965 24 4 0.9898153 3 negt_558 1.00000000 1.000000000 20 8 0.7502920 4 negt_1057 0.68662969 0.377646329 24 14 0.6292402 5 cosp_3926 0.02880291 0.004320437 24 4 1.7893915 6 cosp_3960 0.02880291 0.001933622 28 12 2.4629848 7 cosp_2183 0.02880291 0.004320437 28 2 0.9186064 8 negt_2142 1.00000000 1.000000000 24 12 0.7223294 9 pekp_4656 0.08770330 0.017540660 24 20 1.7087341 10 pekp_4740 0.09273990 0.032458965 28 16 1.5126607 11 negt_5430 1.00000000 1.000000000 20 0 0.2726404 12 negt_5514 1.00000000 1.000000000 28 24 0.8294560 13 negt_6004 1.00000000 1.000000000 28 4 0.6008885 14 squa_7485 0.24211385 0.096845539 24 22 3.0772840 15 negt_6385 1.00000000 1.000000000 28 18 0.6149837 16 squa_6903 0.09273990 0.032458965 24 20 3.2208690 17 negt_8645 1.00000000 1.000000000 28 20 0.3630139 18 negt_8219 0.49574766 0.247873828 24 16 1.1203698 19 negt_8578 1.00000000 1.000000000 24 12 0.6831608 20 negt_9637 0.93150729 0.558904371 20 16 0.3930070
