PPBLOG

　いくつかのグループ間（の平均値）に一定の順序性・傾向性がみられるとき、これらの群間において有意な差があるか検定したい。こんなときは、ヨンキー・タプストラのトレンドテスト(Jonckheere-Terpstra trend test, 単にJonckheere検定とも)が便利だ。こういう場面は、結構あると思う。薬剤の量に応じたある因子の反応データだとか、SNPのタイピング(例えば、AA型/AG型/GG型に傾向性がある場合とか)だとか、病期の重症度に応じた予後だとか。

　一般的に、多群間の検定には、One-way ANOVAであったり、Kruskal-Wallis testがよく使われ、（post-hocとしての）多重比較には、Tuker-CramerやBonferroni/Dunn、あるいはSteel-Dwass検定が使われる。多群の検定というのは、検定に対する思想に左右される部分があり(個人的にはいきなり多重比較検定をしたい人）、適応に頭を悩ますこともある。

　その点、トレンドテストというのは、一定の傾向性があるか否かを判定する、非常に分かりやすい検定法だ。薬剤量に応じた傾向性だとか、臨床病期に応じた傾向に対する情報というのは、実地において有用であろう。

　Rでは、ヨンキー・タプストラ検定のためのライブラリがあるので、それを利用する。マニュアルはここ 。もっとも、インストールした後は、？ JT.testで良いけど。

source("http://bioconductor.org/biocLite.R")
biocLite("SAGx") # パッケージのインストール
library(SAGx) #読み込み

　ここで、添付しているようなサンプルデータがあるとしよう。35匹のマウスについて、24種類の遺伝子の発現を調べたものである。35匹のマウスは、とある重症度の指標に基づいて、3群に分けられている（第1群: 0, 第2群: I-II, 第3群: III-IV）。第1群にはMouse 1-11、第2群には Mouse 12-22、第3群にはMouse 23-35のマウスが属する。それぞれ24種類の遺伝子について、これら3群で一定の傾向性があるかどうか検定したい。

　まずは、テキストファイル(データ間はタブ区切り)からデータの読み込み。

dat <- read.table("JT-sample.txt", header=TRUE, row.names=1, sep="¥t")

マウス3群に対して、適当にラベルを与える。ここでは1,2,3にした。

label <- ordered(c(rep(1, 11), rep(2, 11), rep(3, 13))) # それぞれ11匹、11匹、13匹

実際の検定は、JT.testでOK。 引数のalternativeについては、傾向性は、増加と減少、共に考えられるのでtwo-sidedにした。

jt <- JT.test(data = dat, class = label, labs = c("0", "I-II", "III-IV"), alternative = "two-sided")

この結果のP値を取り出すには、jt$p.valueでOK。結果をテキストファイルに落とすために、以下のように見やすいように適当にデータを作成しよう。

tmp <- cbind(as.matrix(colnames(dat), 1, 24), sprintf("%0.3g", jt$p.value), rank(jt$p.value))

最初に、元のテキストでは見出しの最初の行にあった遺伝指名(Gene 1～Gene 24)を取り出して、それを列にするための指定をしている。もっと簡単な方法がありそうだけど、知らないので。2番目に有効数字3桁ほどに丸めたP値、最後にP値のランク情報を付けておく。これをテキストに保存するには、

write.table(tmp, "JT-sample P.txt", sep = "¥t", append=F, quote=F, row.names=F) 

でOK。この結果は、以下の通り。

V1	V2	V3
Gene.1	0.00181	9
Gene.2	1.95e-05	5
Gene.3	0.477	23
Gene.4	0.0715	15
Gene.5	0.000127	7.5
Gene.6	0.0957	17
Gene.7	2.05e-07	1
Gene.8	0.134	19
Gene.9	0.296	22
Gene.10	0.00258	10
Gene.11	0.0174	13
Gene.12	0.000127	7.5
Gene.13	0.0957	16
Gene.14	0.209	20
Gene.15	0.0205	14
Gene.16	0.596	24
Gene.17	0.22	21
Gene.18	1.26e-06	2
Gene.19	2.48e-06	3
Gene.20	6.86e-06	4
Gene.21	0.00462	11
Gene.22	3.11e-05	6
Gene.23	0.00641	12
Gene.24	0.112	18

　ランキング1位の遺伝子は、Gene.7で、そのP値は、2.05e-07とかなり小さい。この遺伝子の実際のデータ分布を視覚化してみよう。傾向を見せるには、いわゆる箱髭図が良いだろう。表示させて、PNGファイルに落とす。

boxplot(dat$Gene.7 ~ label, names=c("0", "I-II", "III-IV"))
dev.print(png, file="JT-Gene7.png", width=350, height=350)

　こんな感じで、簡単に使えます。ヨンキー・タプストラ検定は、どうも有意なP値が出やすい気がする。

　

　前回は、エクセルのコピーした数値群をクリップボードから読み取って、それをいじるということをしてみたが、相手にするデータが多いときに、逐一それをやっていると日が暮れてしまう。なので、テキストファイルにしたエクセルデータ(タブ区切り)を読み取って、P値などの計算を自動化させて、結果をファイルに出力するようにしてみたい。

　添付のようなサンプルデータがあるとする。32人の患者さんについて、ある治療前後の23個の遺伝子の発現量を調べたものだ(あくまでサンプルデータ)。形式は、第1列に患者IDがあって、CR1_pre, CR2_pre, ..., CR1_post, CR2_post, ..., NR1_pre, NR2_pre, ..., NR1_post, NR2_post という具合に、最初の10行分が治療で完治した(CR)患者さん5人の治療前後のデータで、続いて、11人分の治療に不応答(NR)だった、22行の治療前後のデータが並んでいる。第1行は、調べた23個の遺伝子がGene 1, Gene2, ...と並んでいる。

　このテキストデータは、エクセルで作成してテキストエディタに貼り付けて保存しただけだ。このデータをRに読む込ませるには、以下のようにする。

data <- read.table("sample.txt", header=TRUE, row.names=1, sep="¥t")

以下で、読み込んだデータの行数と列数を表示してくれる。

dim(data) # 結果は[1] 32 23

各列のデータは、例えばdata[1:5, 10] # 10番目の遺伝子の、CR群治療前の発現データみたいに簡単にアクセスできるので、治療前後の対応のあるpaired t-testをして、ファイルに出力してみよう。

data <- read.table("sample.txt", header=T, row.names=1, sep="¥t") 
for (i in 1:23){ # 1列から23列まで同じ操作を繰り返す
 CR_PRE <- data[1:5, i]
 CR_POST <- data[6:10, i]
 NR_PRE <- data[11:21, i]
 NR_POST <- data[22:32, i]
 p1 <- t.test(CR_PRE, CR_POST, paired=T) # 対応のあるt検定
 p2 <- t.test(NR_PRE, NR_POST, paired=T)
 p_cr[i] = sprintf("%0.3g", p1$p.value) # 有効数字3桁で十分かな
 p_nr[i] = sprintf("%0.3g", p2$p.value)
}
out = cbind(colnames(data), p_cr, p_nr) # 列の遺伝子名を取り出し、第1列へ
write.table(out, "foo.txt", sep="¥t", append=F, quote=F, row.names=F) # foo.txtというテキストファイルに結果を出力

出力された結果は、以下の通り。

	p_cr	p_nr
Gene.1	0.175	0.821
Gene.2	0.168	0.0689
Gene.3	0.0254	0.664
Gene.4	0.222	0.102
Gene.5	0.482	0.66
Gene.6	0.00177	0.783
Gene.7	0.342	0.511
Gene.8	0.0933	0.604
Gene.9	0.199	0.169
Gene.10	0.24	0.767
Gene.11	0.324	0.528
Gene.12	0.31	0.5
Gene.13	0.239	0.897
Gene.14	0.721	0.0492
Gene.15	0.398	0.102
Gene.16	0.21	0.347
Gene.17	0.574	0.0541
Gene.18	0.484	0.19
Gene.19	0.118	0.38
Gene.20	0.338	0.355
Gene.21	0.681	0.991
Gene.22	0.0658	0.579
Gene.23	0.0155	0.816

ふむふむ。p_crには、CR群のP値が収納されているけど、0.05以下を抜き出したいときは、以下でOK。

p_cr[p_cr < 0.05] # 結果は、[1] "0.025" "0.002" "0.015"

ただ、これだけだと、数値のみでどの遺伝子なのかが分からない。上で定義したoutというデータフレーム？に対しても、この手の条件抽出が出来るようだ。

> out[p_cr < 0.05,] # 最後にカンマを付ける
               p_cr    p_nr   
[1,] "Gene.3"  "0.025" "0.664"
[2,] "Gene.6"  "0.002" "0.783"
[3,] "Gene.23" "0.015" "0.816"

なるほどねぇ。いやぁ、便利だ。

　さて、ここでは、対応のあるデータ群の検定ということで何も考えずにt検定を行ったが、もともとt検定はパラメトリックなデータ群に対して使われるものである。はたして、sample.txtのデータは、正規分布に従っているのだろうか？ 例えば、CR群でP値が一番小さかった3番目の遺伝子の、基本統計量(とりあえず全サンプルで)を求めてみると、

> gene3 <- data$Gene.3
> summary(gene3)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0708  0.4340  0.7050  1.6300  1.0900 14.3000

とても正規分布に従っているとは思えない。実際にヒストグラムを描いてPNG画像として出力してみよう。

hist(gene3)
dev.print(png, file="gene3.png", width=350, height=350)

これで、以下のPNG画像が作業ディレクトリに保存される。簡単ですなぁ。

　こういうデータ分布に対するt検定の結果は、いかにも心許ない。こういう場合は、対数変換をすると良い場合がある。いつでも対数変換をすれば良いという訳ではないけれど、この場合は、対数変換の適応にはそれなりの根拠があるものとしよう。Rでは底が２の対数変換はlog2()で簡単に与えられる。これで同様にヒストグラムを出力してみよう。

hist(log2(gene3)) # log-transformation (base 2)
dev.print(png, file="gene3-2.png", width=350, height=350)

　もともとのサンプルサイズが小さいので、厳密なことは言えないけど、明らかに対数変換をした方が正規分布に近付いている。なので、対数変換を施した全データについて、再度対応のあるt検定をやってみよう。

data <- read.table("sample.txt", header=T, row.names=1, sep="¥t")
data <- log2(data) # 全データを対数変換
for (i in 1:23){
 CR_PRE <- data[1:5, i]
 CR_POST <- data[6:10, i]
 NR_PRE <- data[11:21, i]
 NR_POST <- data[22:32, i]
 p1 <- t.test(CR_PRE, CR_POST, paired=T)
 p2 <- t.test(NR_PRE, NR_POST, paired=T)
 p_cr[i] = sprintf("%0.3g", p1$p.value)
 p_nr[i] = sprintf("%0.3g", p2$p.value)
}
out = cbind(colnames(data), p_cr, p_nr)
write.table(out, "foo-log2.txt", sep="¥t", append=F, quote=F, row.names=F)

　この結果は、以下の通り。

	p_cr	p_nr
Gene.1	0.269	0.273
Gene.2	0.0335	0.0874
Gene.3	0.0476	0.397
Gene.4	0.224	0.0539
Gene.5	0.462	0.995
Gene.6	0.00279	0.921
Gene.7	0.873	0.296
Gene.8	0.0271	0.79
Gene.9	0.108	0.894
Gene.10	0.0179	0.381
Gene.11	0.0296	0.982
Gene.12	0.131	0.711
Gene.13	0.253	0.643
Gene.14	0.68	0.0136
Gene.15	0.42	0.0989
Gene.16	0.166	0.964
Gene.17	0.671	0.283
Gene.18	0.44	0.186
Gene.19	0.0358	0.333
Gene.20	0.423	0.687
Gene.21	0.764	0.945
Gene.22	0.0686	0.435
Gene.23	0.000343	0.79

　さっきと結構違うねぇ。同様にP値が0.05以下を抜き出してみると、

> out[p_cr < 0.05,]
                p_cr       p_nr    
[1,] "Gene.2"  "0.0335"   "0.0874"
[2,] "Gene.3"  "0.0476"   "0.397" 
[3,] "Gene.6"  "0.00279"  "0.921" 
[4,] "Gene.8"  "0.0271"   "0.79"  
[5,] "Gene.10" "0.0179"   "0.381" 
[6,] "Gene.11" "0.0296"   "0.982" 
[7,] "Gene.19" "0.0358"   "0.333" 
[8,] "Gene.23" "0.000343" "0.79"  

ふーん。ここでは、対応あるデータということでt検定を用いたが、もともとデータ分布が、ノンパラメトリックなものであるならば、Wilcoxon signed rank testを最初から使うのも手だ。まぁ、あまり個人的にこの検定法は好きではないけれど。これは、次回のメモに取っておこう。

---

Rは、言わずと知れた、高度な統計解析が出来るフリーのソフト。Rは、まだ初心者なので、備忘録も兼ねてメモメモメモ。

	Pre	Post
NR	2595	3034
NR	1526	1576
NR	420	488
NR	2454	3054
NR	2171	2282
NR	2571	2034
NR	623	644
NR	430	433
NR	635	652
R	5382	2657
R	4503	2213
R	4312	3509
R	3392	3000

　上のようなエクセルの13行2列のテーブルがあって、2列に渡る数値の部分をコピー。これをこの形式のままRに読み込むには、

> data <- matrix(scan("clipboard"), ncol=2, byrow=T)

とする。すると、

> data
      [,1] [,2]
 [1,] 2595 3034
 [2,] 1526 1576
 [3,]  420  488
 [4,] 2454 3054
 [5,] 2171 2282
 [6,] 2571 2034
 [7,]  623  644
 [8,]  430  433
 [9,]  635  652
[10,] 5382 2657
[11,] 4503 2213
[12,] 4312 3509
[13,] 3392 3000

　さて、ここでノンレスポンス群(NR)とレスポンス群(R)とのそれぞれについて、PREとPOSTとで対応のあるt検定をしたい。NR群(最初の9行)のPREのデータは、以下のようにして取得できる。

> NR_PRE <- data[1:9, 1]

> NR_PRE
[1] 2595 1526  420 2454 2171 2571  623  430  635

POSTのデータは同様に、

> NR_POST <- data[1:9, 2]

対応のあるt検定は、

> t.test(NR_PRE, NR_POST, paired=T)

        Paired t-test

data:  NR_PRE and NR_POST
t = -0.8163, df = 8, p-value = 0.438
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -328.0993  156.5438 
sample estimates:
mean of the differences 
              -85.77778 

　レスポンス群についても同様に、

> R_PRE <- data[10:13, 1] # レスポンス群のPREを取得

> t.test(R_PRE, R_POST, paired=T)

        Paired t-test

data:  R_PRE and R_POST 
t = 2.7491, df = 3, p-value = 0.0708
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -244.7542 3349.7542 
sample estimates:
mean of the differences 
                 1552.5

　レスポンス群のサンプルサイズが小さいので、有意差は出ないがその傾向にはある。似たようなことを繰り返しやるのであれば、ファイル自体を読み込ませて、もっと包括的にP値を取得した方が早いのは間違いないけど、まぁ練習として。R初心者だしね。あるいは、以下のような関数を定義してやるのもそう労力ではない。

ptt <- function(){ # paired t-test via Excel clipboard
 data <- matrix(scan("clipboard"), ncol=2, byrow=T)
 NR_PRE <- data[1:9, 1]
 NR_POST <- data[1:9, 2]
 NR =  t.test(NR_PRE, NR_POST, paired=T)
 
 R_PRE <- data[10:13, 1]
 R_POST <- data[10:13, 2]
 R =  t.test(R_PRE, R_POST, paired=T)
 
 cat(NR$p.value, "¥t", R$p.value, "¥n") # それぞれのP値の出力
}

　この関数を定義しておけば、後は、エクセルの範囲をコピーした状態で、ptt()と打ち込めば、それぞれの群のP値が出力される。

Rは、非常に柔軟な言語で、10人いれば10通りの記述が出来る点がおいしいところ。早く使いこなしたいなぁ。

---