oscillation (振動)を検定する

MikuHatsune2017-12-18

この記事は
今年読んだ一番好きな論文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 は、ヨンクヒール・タピストラ・ケンダール検定の頭文字である。
上昇トレンドの検定 ヨンクヒール・タプストラ検定 - 驚異のアニヲタ社会復帰への道
ケンダール検定は、データのペアワイズの差分符号を調べて検定するが、これにヨンクヒール・タピストラ検定という、データに順序があるときの制約をいれて統計量を算出するノンパラメトリックな検定の合わせ技、ということになっている。
いま、時系列データ\bf{x}=\{x_1,\dots,x_n\}, \bf{y}=\{y_1,\dots,y_n\} に対して
\tau(\bf{x},\bf{y})=\frac{\sum_{1\leq i \leq j \leq n}sgn(x_j-x_i)\cdot sgn(y_j-y_i)}{\frac{1}{2}n(n-1)}
という統計量を計算する。これは組み合わせになるのでめちゃめちゃ計算量が増えるが、がんばれば正規分布に近似できるらしく、なんかいろいろ計算量を減らす工夫があるらしい。
また、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