Bioinformaticsのお勉強

ウィーン大学で生命科学者をしている今鉄男と申します。生命の尊さ、美しさがどのようなメカニズムで生じるのかに興味があります。科学の方法論を用いて、このような問いに応えたい、私はこう思って医学生物学の基礎研究のトレーニングを受けてきました。生命を科学的手法を用いて理解を試みる上で、genomeを始めとした種々の大量データの処理が必要不可欠であることを痛感しました。また、生命科学が物理学、数学、統計学、有機化学などの種々の学問と深い関わりを持つことを実感しました。そのため、このブログは広範囲の学問領域に関しての記事を載せています。日々の学習内容を文書に書き残し、それを読み返すことによって、体系化された知識を身に付けることを目標としています。どうぞよろしくお願いします。

コンソート声明に基づいたランダム化試験のデザイン

無作為ランダム化試験は臨床医学研究の分野では信頼性の高いエビデンスを確立する上で、必須となってきています。

ランダム化比較試験の統一化しようとの試みで、CONSORT statment(Consolidated Statndards of Reporting Trials)がという宣言が出されています。

以下に、そのCONSORT statmentに基づいたランダム化試験のデザインのフローチャートを示します。

しっかりとしたクライテリアの下で品質の保証された臨床研究が行われていって欲しいものです。

一つの母集団からのサンプル抽出を味わう

$ vim command.in

dat <- rnorm(10000, mean=0, sd=1)
#母集団のヒストグラムに母集団の平均値を記述(green)
png("120730_hist_1.png")
hist(dat, breaks=100)
abline(v=mean(dat), col="green")
dev.off()

#サンプルの平均値をmagendaで記述。
#サンプルの平均値のばらつきを見極める。
png("120730_hist_2.png")
hist(dat, breaks=100)
abline(v=mean(dat), col="green")
for(i in 1:10)
{
abline(v=mean(dat[sample(1000,10)], col="magenda"))
}
dev.off()

$ fg
> source("command.in")

サンプルの平均値は母平均のそばを漂っているイメージである。
この事実から言えることは、サンプル
サンプリングした標本に対しての統計学的なデータの要約方法について記述する。

#サンプリング(1000個の母集団の要素中から10個を抽出）

> sample <- dat[sample(1000,10)]
> sample
[1] 0.3197940 1.8119807 -0.2420026 -0.2343419 -0.6559658 0.5846391
[7] -0.7373418 0.5727909 0.7129338 0.0694458

#標本数

> length(sample)

[1] 10

#平均値の計算

> mean(sample)

[1] 0.2201932

#標本標準偏差の計算
#Rにはデフォルトでは不偏分散
> variance <- function(x){var(x)*(length(x)-1)/length(x)}
> sqrt(variance(sample))
[1] 0.7191055

#標本の95%信頼区域を計算(ついでに平均値が0であるかの検定も行う)

> t.test(sample, mu=0, alternative="two.sided")

One Sample t-test

data: sample
t = 0.9186, df = 9, p-value = 0.3823
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
-0.3220500 0.7624364
sample estimates:
mean of x
0.2201932

以上より、以下の記述が可能です。

○(平均(標準偏差)[範囲])
グループA(n=20)
データ(無次元) 0.22(0.72)[-0.32 - 0.76 ]

○(平均±標準偏差)[範囲])
グループA(n=20)
データ(無次元) 0.22 ± 0.72[-0.32 - 0.76 ]

#ちなみに、95%信頼区間が0をまたいでいる[-0.32 - 0.76 ]ということは、t検定により、母平均が0とは異なるとまでは言えないということと同値である。かりに、p値が0.05を下回るようならば、95%信頼区間が正負いずれかの方向に偏ることになる。

臨床的に有意な差とは何か　〜サンプル数は多すぎても少なくてもいけない〜

２群の平均値の比較において大切なことは、「有意差が出る」ことではなくて「実質的に（ex.臨床的に、実験的に）有意な差を検出する」ことです。

例えば、A,B二つの集団の身長の平均値の差について考えます。

二群の平均値の差が0.1cm程度しかなかったとき、その差を膨大なサンプリングによって検出することに果たして意味があるのでしょうか。いやありません。

サンプル数は多すぎても問題なのです。

大切なのは、実験を行う前の段階で、どの程度の平均値の差なら実質的に有意であるかを明らかにした上で、必要最小限のサンプルを行うことです。

以下に、サンプル数とp値のシミュレーションの結果を書きます。

$ vim comman.in

result1 <- c(1:99)

for(i in 1:99)

{

#2群の平均値の差が10cmの場合

dat1 <- rnorm(mean=175.0, n=i+1, sd=5)

dat2 <- rnorm(mean=185.0, n=i+1, sd=5)

result1[i] <- t.test(dat1,dat2,alternative="two.side",mu=0, paired=FALSE,var.equal=FALSE,conf.level=0.95)$p.value

}

result2 <- c(1:99)

for(i in 1:99)

{

#2群の平均値の差が５cmの場合

dat1 <- rnorm(mean=175.0, n=i+1, sd=5)

dat2 <- rnorm(mean=180.0, n=i+1, sd=5)

result2[i] <- t.test(dat1,dat2,alternative="two.side",mu=0, paired=FALSE,var.equal=FALSE,conf.level=0.95)$p.value

}

result3 <- c(1:99)

for(i in 1:99)

{

#2群の平均値の差が1cmの場合

dat1 <- rnorm(mean=175.0, n=i+1, sd=5)

dat2 <- rnorm(mean=176.0, n=i+1, sd=5)

result3[i] <- t.test(dat1,dat2,alternative="two.side",mu=0, paired=FALSE,var.equal=FALSE,conf.level=0.95)$p.value

}

n <- c(2:100)

png("120731_n.png")

par(mfrow=c(1,3))

plot(n,result1)

plot(n,result2)

plot(n,result3)

dev.off()

$ R

> source("command.in")

平均値の差が小さくなると、統計的に有意であると判断されるのに多くのサンプルが必要であることがわかりました。

次に、平均値の差が1cmの時に、さらにサンプル数を多くして検証します。

対照群として、平均値の差がない群の比較を提示しています。

> cntrl-z

$ vim command.in.2

result4 <- c(1:999)

for(i in 1:999)

{

dat1 <- rnorm(mean=175.0, n=i+1, sd=5)

dat2 <- rnorm(mean=176.0, n=i+1, sd=5)

result4[i] <- t.test(dat1,dat2,alternative="two.side",mu=0, paired=FALSE,var.equal=FALSE,conf.level=0.95)$p.value

}

result5 <- c(1:99)

for(i in 1:999)

{

dat1 <- rnorm(mean=175.0, n=i+1, sd=5)

dat2 <- rnorm(mean=175.0, n=i+1, sd=5)

result5[i] <- t.test(dat1,dat2,alternative="two.side",mu=0, paired=FALSE,var.equal=FALSE,conf.level=0.95)$p.value

}

n <- c(2:1000)

png("120731_n_2.png")

par(mfrow=c(1,2))

plot(n,result4)

plot(n,result5)

dev.off()

$ fg

> source("command.in.2")

確かに、1000サンプル近く抽出すればそれなりに有意差が出てくるようですが、この作業に何の意味があるのでしょうか（いやない）。

あらかじめ有意であると見なすことのできる差を想定した上でその差を検出するのに必要最低限なサンプルを確保することが王道ではないかと思います。

やたらめったら、サンプルを増やしてp<0.05としたところで、その差が実際のところどの程度の影響を持つのか考える必要があろうかと思います。

一標本からのサンプルの抽出についての検討

一標本からのサンプルの抽出について検討しました。

今、身長の平均値175cm, 標準偏差が5cmの成人男性の集団（合計1000人)を考えます。
(対象の集団は175±5cm(n=1000)とも表現可)

Rを使って、この集団のデータを擬似的に作成します。

#Rの起動
$ R
#データの取得
> dat <- rnorm(mean=175, n=1000, std=5)

#データの要約を表示

> summary(dat)
Min. 1st Qu. Median Mean 3rd Qu. Max.
160.3 171.5 175.2 175.1 178.4 192.5

###解釈###
最大値が192.5cm, 最小値が160.3cm, 平均値が175.2cm
#このデータは偏りがないはずであるが、偏りがもしもある場合は、中央値と四分位範囲、範囲について議論するのがよいとされる。
第一四分位数が171.5cm, 第三分位数が178.4cm(四分位範囲：171.5-178.4cm)
（中央値（四分位範囲）[範囲]): 175.2(171.5 - 178.4)[160.3 - 192.5] cm

ここで、この集団から無作為に10人抽出して身長を測定することを考えます。
この作業から、母集団の平均値をどの程度推定できるか検討してみます。

#ランダムに50人のサンプル抽出

> sample <- dat[sample(1000,10)]

#標本平均を算出

> mean

[1] 174.9052

#ヒストグラムに母平均値を示す垂線を重ねて描写

png("120730_hist_1.png")

hist(dat, breaks=100)

abline(v=mean(dat), col="green")

dev.off()

#繰り返し標本抽出と標本平均の算出を行う。

> mean(dat[sample(1000,10)])

[1] 175.2295

> mean(dat[sample(1000,10)])

[1] 173.3754

> mean(dat[sample(1000,10)])

[1] 179.5982

> mean(dat[sample(1000,10)])

[1] 177.0066

> mean(dat[sample(1000,10)])

[1] 173.4215

#ヒストグラムの上に標本平均の垂線を重ねて描写(計5回のサンプリング)

png("120730_hist_2.png")

hist(dat, breaks=100)

abline(v=mean(dat), col="green")

abline(v=mean(dat[sample(1000,10)], col="magenda"))

dev.off()

sampleオブジェクト(対象の母集団から50人の標本抽出を行ったもの)の平均値が、175と比べて統計的に有意に異なるかを検定する。(一標本t検定)

> t.test(sample, mu=175, alternative="two.sided")

One Sample t-test

data: sample

t = -0.1224, df = 49, p-value = 0.9031

alternative hypothesis: true mean is not equal to 175

95 percent confidence interval:

173.3477 176.4626

sample estimates:

mean of x

174.9052

この結果、p=0.9031<0.05であり、統計的に有意に母集団の平均が175cmと異なるとは言えない。(母平均が175cmであるとは言えないが。)

そして、母平均は173.3477cmから176.4626cmの間に95%の確率で存在する。

(95%信頼区域： 173.3477 - 176.4626 cm)

Cent OSサーバーにRをインストールする

Cent OSはyumでRをインストールができません。

そのため、サードパーティーを利用してRをインストールしてあげる必要があります。

epelというサードパーティーリポジトリにRが収録されているそうです。

今日はそのための方法を書きます。

#OSのバージョンを調べる

$ cat /etc/redhat-release
CentOS release 5.7 (Final)

#管理者権限で、epelリポジトリを登録する

$ sudo rpm -Uvh http://dl.fedoraproject.org/pub/epel/5/i386/epel-release-5-4.noarch.rpm

#epelレポジトリの追加を確認する

$ sudo yum repolist

repo id repo name status

CRAN http://cran.md.tsukuba.ac.jp/ 5

base CentOS-5 - Base 3,591

epel Extra Packages for Enterprise Linux 5 - x86_64 7,107

extras CentOS-5 - Extras 273

updates CentOS-5 - Updates 770

#あとはyumを使ってインストールするだけ

$ sudo yum install R

#最後はepelレポジを無効化しておきます。
#サードパーティーは必要に応じてつか

$ sudo vim /etc/yum.repos.d/epel.repo

1 [epel]

2 name=Extra Packages for Enterprise Linux 5 - $basearch
3 #baseurl=http://download.fedoraproject.org/pub/epel/5/$basearch
4 mirrorlist=http://mirrors.fedoraproject.org/mirrorlist?repo=epel-5&arch=$basearch
5 failovermethod=priority
6 enabled=1　　　　########ここを0に変更#########
7 gpgcheck=1
8 gpgkey=file:///etc/pki/rpm-gpg/RPM-GPG-KEY-EPEL

#以下のオプションをつければ、epelを無効化している状態でも必要に応じて、使用することができる。

–enablerepo=epel

#epelレポジトリの無効化を確認する

$ sudo yum repolist

repo id repo name status

CRAN http://cran.md.tsukuba.ac.jp/ 5

base CentOS-5 - Base 3,591

extras CentOS-5 - Extras 273

updates CentOS-5 - Updates 770

#Rの動作をチェックします。

$ R --version

R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows"

#つい最近のバージョンが入っていることがわかりました。

#起動します。

$ R

R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows"

ISBN 3-900051-07-0

Platform: x86_64-redhat-linux-gnu (64-bit)

Rは、自由なソフトウェアであり、「完全に無保証」です。

一定の条件に従えば、自由にこれを再配布することができます。

配布条件の詳細に関しては、'license()'あるいは'licence()'と入力してください。

Rは多くの貢献者による共同プロジェクトです。

詳しくは'contributors()'と入力してください。

また、RやRのパッケージを出版物で引用する際の形式については

'citation()'と入力してください。

'demo()'と入力すればデモをみることができます。

'help()'とすればオンラインヘルプが出ます。

'help.start()'でHTMLブラウザによるヘルプがみられます。

'q()'と入力すればRを終了します。

#うまくいきました。

登録: 投稿 (Atom)

Bioinformaticsのお勉強

コンソート声明に基づいたランダム化試験のデザイン

一つの母集団からのサンプル抽出を味わう

臨床的に有意な差とは何か　〜サンプル数は多すぎても少なくてもいけない〜

一標本からのサンプルの抽出についての検討

Cent OSサーバーにRをインストールする

自己紹介

過去のブログ♪♪♪

コンソート声明に基づいたランダム化試験のデザイン

一つの母集団からのサンプル抽出を味わう

臨床的に有意な差とは何か 〜サンプル数は多すぎても少なくてもいけない〜

一標本からのサンプルの抽出についての検討

Cent OSサーバーにRをインストールする

臨床的に有意な差とは何か　〜サンプル数は多すぎても少なくてもいけない〜