script01.nb.html | |
File Size: | 838 kb |
File Type: | html |
まずはホームディレクトリを設定します。
home <- [配布したフォルダをコピーしたパス]
setwd(home)
library(data.table)
Warning message:
In scan(file = file, what = what, sep = sep, quote = quote, dec = dec, :
EOF within quoted string
t検定について説明します。以下のように、10個の遺伝子に対して二つの条件下(例えば、Tumor / Normal)で5サンプルずつの遺伝子発現データ考えます。目標は、条件間で発現量の異なる遺伝子群を同定することです。
set.seed(100)
n <- 10
p <- 5
cond1 <- matrix(rnbinom(n*p,2,mu = 300),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(1,p))))
cond2 <- matrix(rnbinom(n*p,2,mu = 100),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(2,p))))
data_example <- cbind(cond1,cond2)
Condition1 | Condition1 | Condition1 | Condition1 | Condition1 | Condition2 | Condition2 | Condition2 | Condition2 | Condition2 | |
---|---|---|---|---|---|---|---|---|---|---|
Gene1 | 123 | 91 | 178 | 86 | 327 | 92 | 86 | 75 | 46 | 329 |
Gene2 | 419 | 250 | 648 | 210 | 264 | 210 | 93 | 43 | 89 | 62 |
Gene3 | 277 | 119 | 397 | 148 | 538 | 87 | 58 | 87 | 46 | 33 |
Gene4 | 306 | 198 | 530 | 139 | 428 | 163 | 246 | 54 | 96 | 135 |
Gene5 | 243 | 211 | 172 | 152 | 47 | 102 | 3 | 233 | 58 | 49 |
Gene6 | 205 | 299 | 620 | 202 | 206 | 127 | 64 | 128 | 52 | 284 |
Gene7 | 125 | 56 | 361 | 549 | 80 | 345 | 62 | 189 | 39 | 45 |
Gene8 | 513 | 14 | 217 | 451 | 179 | 89 | 37 | 294 | 125 | 169 |
Gene9 | 71 | 525 | 36 | 276 | 83 | 24 | 141 | 70 | 29 | 166 |
Gene10 | 288 | 151 | 1005 | 108 | 275 | 37 | 96 | 228 | 7 | 195 |
二つの条件間の違いを見る方法はいくつか考えられます。一つは、log fold change (LFC)を計算することです。遺伝子\(i\)に対するLFCとは
\[\mbox{logFC}_i=\mbox{log2}(M_{i,1} / M_{i,2})\]
(\(M_{i,j}\)は遺伝子\(i\)条件\(j\)の発現量の平均です。
M1 <- apply(cond1,1,function(x)mean(x))
M2 <- apply(cond2,1,function(x)mean(x))
(logFC <- log2((M1 + 1) / (M2 + 1)))
Gene1 Gene2 Gene3 Gene4 Gene5 Gene6 Gene7 Gene8 Gene9 Gene10
0.3557164 1.8390281 2.2314946 1.2001075 0.8831863 1.2195792 0.7797122 0.9395588 1.1951303 1.6894567
LFCはどの程度、違うかを表す指標ですが、どの閾値以上であれば差がある
かの明確な基準はありません。また遺伝子ごとにデータのばらつき方は異なり、得られたLFCが本当に偶然かどうかも見定める必要がありますがそれはできません。そのために、統計的仮説検定という枠組みを用いて、ばらつきを確率的に考えて、本当にその差があるのかを評価します。t-検定は、二群間の平均の差を検定する代表的な検定手法です。
早速、t-検定を実行してみましょう。遺伝子1に対してt-検定を実行してみます。
x <- cond1[1,]
y <- cond2[1,]
(result_t <- t.test(x,y,var.equal = TRUE))
Two Sample t-test
data: x and y
t = 0.51975, df = 8, p-value = 0.6173
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-121.6617 192.4617
sample estimates:
mean of x mean of y
161.0 125.6
p値はいくらかというと、
result_t$p.value
[1] 0.6173102
となり、有意水準5%で考えると有意ではありません。Rでt-検定を行うのは非常に簡単です。ただし、その仕組みを知っていないと、あまり意味も分からずにp値を振り回したり、振り回されてしまいます。さきほどの数行で何をやったのかの意味を理解してみようと思います。
統計的仮説検定では、まず帰無仮説と対立仮説を立てます。帰無仮説とは最終的に否定
したい事柄を、対立仮説にはそれと対比をなす事柄を考えます。今回の場合だと、条件1と条件2で、発現量が異なる
遺伝子群を同定したいので、仮説を以下のようにします。 \[\mbox{帰無仮説}:M_{i,1}=M_{i,2}\] \[\mbox{対立仮説}:M_{i,1}\neq M_{i,2}\]
次に、帰無仮説が本当ではないことを確かめるために、帰無仮説を定量化した検定統計量
を求めます。
t-統計量
\[\mbox{t-統計量} = \frac{\mbox{群間の平均の違い}}{\mbox{各群のばらつき}}\]
先ほどのt統計量は、
result_t$statistic
t
0.5197481
です。
このt-統計量は単なる数値ではありません。自分の確率的なばらつきかたを表す確率分布
を背負っています。特に帰無仮説の確率分布なので、帰無分布
と呼びます。仮説検定のポイントはこの帰無分布を正しく理解することです。t統計量の従う分布はt分布という名前がついており、その形はサンプル数のみによって決まります。Rでは簡単にt分布からデータを生成させることができるので、実際に見てみましょう。
data_t <- rt(10000,df = n - 2)
hist(data_t,freq = F,xlab="t",ylab="Probability")
上図はt分布からt統計量を10000個、仮想的に生成したときの頻度分布です。全ての頻度の和を取ると、1になるように正規化しています。これが先ほど実行したt検定の帰無分布です。帰無仮説は二群間の平均の違い=0、すなわち帰無仮説に対応するt統計量が0なので、理論上はゼロになるはずですが、分布上は、確率的なばらつきを持っています。別の言い方をすると、もしも、10000回(実際にはそんなことはできませんが)、繰り返しt統計量を計算した時のばらつき方を表しています。
繰り返しになりますが、この分布は帰無仮説に対応したt統計量のばらつきを表しています。では、対立仮説に対応するt統計量の値(t=-2.026503)はどのあたりかというと、
hist(data_t,freq = F,xlab="t",ylab="Probability")
abline(v=result_t$statistic,col="red")
確かに帰無仮説に対応する0ではないですが、誤差の範囲内と言えそうです。これが確率的に
誤差の範囲かどうかを判定するには、有意水準
という閾値を設けます。具体的に有意水準5%といったときには、
図中青色の線が有意水準に対応する閾値を表しており、その両側に対応する確率の和が有意水準に対応する5%となっています。有意水準5%とは、偶然0ではない(帰無仮説が正しい)とする閾値以上となる確率が高々5%であると言っています。仮に間違って帰無仮説が間違っているといったとしても、100回中5回程度だということです。裏を返せば、もしもその閾値を越えていれば、それは帰無仮説の下では、滅多に起こりえないことが起こったということになります。したがって、5%という誤差はあるけれども、今回の帰無仮説は成り立たない
すなわち有意ではないと考えるわけです。
今回は、5%の有意水準を越えていないので、帰無仮説は5%の有意水準で棄却
できません。
ではp値とは何かというと、図中の赤色の線の右側の確率の和です(もし負であれば左側)。先ほど、p値は0.613でしたが、右側だけの有意水準に対応する領域は0.025(2.5%)なので、この値を用いても先ほどと同様に赤色の線が青色の線を越えていないということができます。つまり、p値は有意水準に対応する閾値を越えているかを判断するための値です。
全ての遺伝子に対して、検定してみましょう。
pvals <- NULL
for(i in 1:n){
x <- cond1[i,]
y <- cond2[i,]
result_t <- t.test(x,y,var.equal = TRUE)
pvals[i] <- result_t$p.value
}
pvals
[1] 0.61731016 0.01664043 0.01819819 0.05069772 0.17883478 0.08861926 0.40741858 0.22987814 0.27743056 0.17355748
次に有意水準5%の説明が本当かどうかを確かめてみましょう。確かめ方は簡単です。同じ確率分布から二群の遺伝子発現量を生成させて、5%有意水準でのt検定を10000回繰り返したときに、間違って帰無仮説を棄却
する回数がおよそ500回であれば5%の意味は正しいと言えます。
set.seed(100)
n <- 10000
p <- 5
cond1_many <- matrix(rnbinom(n*p,2,mu = 100),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(1,p))))
cond2_many <- matrix(rnbinom(n*p,2,mu = 100),ncol=p,dimnames=list(paste0("Gene",1:n),paste0("Condition",rep(2,p))))
pvals_many <- NULL
for(i in 1:n){
x <- cond1_many[i,]
y <- cond2_many[i,]
result_t <- t.test(x,y,alternative = "two.sided",var.equal = TRUE)
pvals_many[i] <- result_t$p.value
}
sum(pvals_many <= 0.05) / n
[1] 0.044
確かに、5%くらいの確率で帰無仮説を誤って棄却
しているようです。
今回は、両側
検定と呼ばれる方法で検定を行いました。それがalternative="two.sided
の部分です。両側の意味は、対立仮説が\(M_{i,1} < M_{i,2}\)と\(M_{i,1}> M_{i,2}\)(つまり\(M_{i,1}\neq M_{i,2}\))の両方を含めているということです。片側検定の場合は、どちらかを検定することになります。つまり、
\[\mbox{帰無仮説}:M_{i,1}>M_{i,2}\] \[\mbox{対立仮説}:M_{i,1}\neq M_{i,2}\]
を検定します。このときの注意は、有意水準の考え方です。両側検定で有意水準5%といった場合には、二つのしきい値は、外側の和が5%になるようなしきい値となりますが、片側検定ですと、片方のしきい値しか用いなので、片方のしきい値の外側の和が5%となるようなしきい値となります。
alternative="greater"
とすれば、対立仮説は\(M_{i,1} > M_{i,2}\)となり、alternative="less"
とすれば対立仮説は\(M_{i,1} < M_{i,2}\)となります。t.test
関数の基本設定ではalternative="two.sided"
となっています。
x <- cond1[1,]
y <- cond2[1,]
(result_t <- t.test(x,y,greater="greater",var.equal = TRUE))
Two Sample t-test
data: x and y
t = 0.51975, df = 8, p-value = 0.6173
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-121.6617 192.4617
sample estimates:
mean of x mean of y
161.0 125.6
var.equal
は2群のばらつきを同じと仮定するか否かを設定します。今回は同じと仮定したのでvar.equal=TRUE
としています。もしvar.equal=FALSE
とすると、ウェルチの検定という方法になります。状況に応じて使い分けてください。
ところで、有意水準5%で検定を100回繰り返すと、有意水準の定義から、確率的には5回、間違って帰無仮説を棄却することになります。遺伝子数は2万近くありますから、遺伝子ごとに検定を行うと20000回の検定を行うことになり、理論上は$20000\times 0.05=1000$
回は間違って、発現の差異がないのに、あると判断してしまうことになります。つまり、擬陽性の遺伝子を1000個程度含んでしまうことになります。統計的仮説検定を繰り返し行うことを多重検定
と言い、擬陽性をいかにして減らすかという問題は統計学では多重検定補正
といって、重要な問題です。擬陽性率をコントロールして擬陽性を減らすために、いくつかの方法が提案されています。研究論文で広く用いられている方法の一つはBenjamini and HochbergによるFalse discovery rate (FDR)です。FDRとは間違って棄却した帰無仮説のうち、本当は帰無仮説が正しい割合です。FDR=5%といった場合には、本当は発現差異がないのに、あると判断してしまった遺伝子1000個のうち、\(1000\times 0.05=50\)個の遺伝子は本当は発現差異がないという仮定で、p値を補正します。詳しい原理はここでは省略します。FDRによる多重検定補正を行うには、p.adjust
関数を用います。まず、10個の遺伝子に対するp値が以下のように与えられたとします。
pvals
[1] 0.61731016 0.01664043 0.01819819 0.05069772 0.17883478 0.08861926 0.40741858 0.22987814 0.27743056 0.17355748
5%有意な遺伝子は
pvals<=0.05
[1] FALSE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
です。
これに対して、FDRによる多重検定補正を行うには以下のように実行します。
(qvals <- p.adjust(pvals,method="BH"))
[1] 0.61731016 0.09099094 0.09099094 0.16899241 0.29805796 0.22154816 0.45268731 0.32839735 0.34678820 0.29805796
FDR5%の結果は、以下のようになります。
qvals <= 0.05
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
p値による結果では有意だったものが、FDRによる多重検定補正で有意ではなくなりました。多重検定補正により有意となる基準が厳しくなるので、実際の解析では擬陽性の遺伝子の数を減らすことができます。