Embedded Files
1. はじめに
生物における情報の比較は、数値の比較になりますので、統計学が極めて重要です。
Rは大規模データ解析が可能であるのみならず、統計学の基礎的な理解をサポートしてくれるフリーのソフトウエアです。
また、下記のサイトに入るとブラウザ上でRを使用することができるようですので、ダウンロードすることなくできます（その場合、自己責任で行ってください）。
https://rextester.com/l/r_online_compiler
https://paiza.io/ja/languages/r
2. Rのダウンロード
CRANよりソフトウエアをダウンロードします。
ミラーサイトのリストからJapanにあるダウンロード先の一つを選んでください。
ダウンロード先のノートPCに合してLinux、Mac、Windowsの選択を行ってください。
baseを選択して、ダウンロードしてください。
インストールの際、使用する言語は日本語を選択してください。
起動時オプションのカスタマイズでは、いいえ(デフォルトのまま)を選択してください。
defaultのまま作業を続け、デスクトップにアイコンが現れれば、完了です。
あるいは、左下のスタートをクリックして、インストールされたリストの中にRがあれば、完了です。
3. Rの起動と入力
Rを起動させてください。
終了するときはq()としてください(本講義では作業スペースの保存は必要ありませんので、いいえを選択してください)。
Rでは大文字と小文字を区別しますので、aとAは違っていると認識します。
また、スペースは無視されますので、q()もq( )も同じです。
c(1, 2, 3)と打ち込んでください。cはconcatenate(連結)を意味しています。
文字は"a"のように""を付けます。例えば、c("A", "b", "abc")は3つの文字からなるベクターを意味しています。もし、""を外して入力した場合、定義していないオブジェクトしてエラーが表示されます。
ただ、Tはtruth(真)、Fはfalse(偽)を意味していますので、そのように認識されます。
例えば、c(1, 2, 3) > 1と入力してみてください。
FALSE TRUE TRUEと表示されます。これは3つの要素の中で1より大きいものが2番目、3番目にあることを意味しています。
あるデータに対して、ある条件を満たすものを抽出する際にはT or Fの情報が必須となるため、これらの文字はあらかじめそのように定義されています。
seq(1, 9)と打ち込んでください。1から9までの整数が表示されます。seqはsequence(連続)を意味しています。
また、defaultでは間隔が1となっていますが、間隔を0.5に変更する場合にはseq(1, 9, 0.5)とします。大きな数から小さな数の場合にはseq(9, 1)としてください。
間隔が1の場合には1:9とすることもできます。
同じパターンを繰り返す(replicate)場合にはrep()を使います。
rep(c(1, 2, 3), 3)
はc(1, 2, 3)を3回繰り返すこと、
rep(c(1, 2, 3), 1:3)
はc(1, 2, 3)の1番目を1回、2番目を2回、3番目を3回繰り返すこと、
rep(1:3, c(3, 2, 1))
は1:3の1番目を3回、2番目を2回、3番目を1回繰り返すことを意味しています。
4. データの収納と四則演算
先ほど、aやAを文字ではなく、定義する必要があることを言いましたが、aやAなどにデータを収納することができます。
例えば、a <- c(1, 2, 3)と打ち込んでください。次に、aと打ち込んでください。aがc(1, 2, 3)のベクターと定義されたことがわかります。
<-は=としても同じ意味です(すなわち、Rでは=は数式のイコールを意味していません。数式のイコールは==となります)。
また、aを異なることで定義しますと、前の定義は残らず、エラー表記もなしに、次の定義へ更新されますので、注意してください。
a+aは各要素の同じ位置における和、
a-aは各要素の同じ位置における差、
a*aは各要素の同じ位置における積、
a/aは各要素の同じ位置における商を示します。
c(a, a)はベクターaとベクターaの連結を示し、和ではないことを認識してください。
sqrt(a)は平方根(square root)、log(a)は自然対数、log2(a)は底が2の対数を意味します。
a > 2はaの要素の中で2より大きいものを判別、
a >= 2はaの中で2以上のものを判別、
a == 2は2と同じものを判別、
a != 2は2ではないものを判別することを意味しています。
aから特定条件を満たすものを抽出するときはa[a > 2]のように[ ]を使います。
a[1]-a[3]はベクターaの1番目から3番目を引くことを意味しています。
差の大きさを知りたい場合にはabs(a[1]-a[3])とすれば、絶対値(absolute value)を知ることができます。
Rにおいてはabs()のようにすでに定義されているコマンドがありますので、よく使うものを示します。
a <- c(1, 20, 300, 4000)を例にすると、
sum of elements（要素の総和）sum(a)
sample size（サンプルサイズ）length(a)
平均はsum(a)/length(a)のことですが、mean(a)
中央値はmedian(a)、要素数が偶数の場合には上位半分の底と下位半分の頂の平均を示します。
最大値はmax(a)
最小値はmin(a)
最大値と最小値の表示はrange(a)
不偏分散は平方和/自由度sum((a-mean(a))^2)/(length(a)-1)と計算できますが、var(a)でわかります。
標準偏差(standard deviation)はsqrt(var(a))と計算できますが、sd(a)でわかります。
5. 標準正規分布
Rでは関数を定義することもできます。
例えば、f(x) = x2はkansu1 <- function(x) x^2
f(x, y) = x2 + y3はkansu2 <- function(x, y) x^2+y^3
標準正規分布(standard normal distribution)はf(x) = (1 / sqrt(2π))exp(-x2/2)の式で表現できます。
これをRを用いて定義するとseikibunpu <- function(x) 1/sqrt(2*pi)*exp(-x^2/2)と表現できます。
ここではseikibunpuという関数を定義しました。seikibunpuと打ち込むと式が出てきます。
この関数がどのようなものか見るためには、seikibunpu(0)、seikibunpu(1)などの結果を見たり、plot(seikibunpu, -5, 5)として見ることができます。
もちろん、標準正規分布はRにすでに登録されていますので、dnormを使います。
先ほどのseikibunpuの形状とdnormを比較してみましょう。
a <- seq(-5, 5, 0.2)
b <- seikibunpu(a)
plot(a, b)
により、seikibunpuの様子を見ることができます。
ここにdnormの形状を重ね書きしましょう。
重ね書きをするには、par(new=T)と入力してください。
その後、plot(a, dnorm(a), type="l", col="red")とすれば、赤色の線(line)で重ね書きされます。
完全に一致していることが確認できます。
なお、plotのタイプには"b"、"c"、"h"、"o"、"s"、"S"などがあります。
色を変えて重ねて見ましょう。
par(new=T)
plot(a, dnorm(a), type="b", col="blue")
par(new=T)
plot(a, dnorm(a), type="c", col="white")
par(new=T)
plot(a, dnorm(a), type="h", col="lightblue")
par(new=T)
plot(a, dnorm(a), type="o", col="yellow")
par(new=T)
plot(a, dnorm(a), type="s", col="pink")
par(new=T)
plot(a, dnorm(a), type="S", col="orange")
Rで使用できる色の名前はcolors()で見ることができます。
6. データの分布を表す
a <- rnorm(1000)
は標準正規分布に従う乱数1000個を発生させaに収納することを意味しています。
rnorm(1000, mean=0, sd=1)におけるmean=0、sd=1がdefaultになっています。
課題１ (1) aのsum, mean, sdを出してください。 (2) plot(1:1000, a)を行って乱数の様子を見てください。次に、-1以上1以下の数値にある比率を求めるためのRにおける実行文を示し、その値を示してください。
hist(a)はaの分布をヒストグラム(histogram)で表します。
hist(a, breaks=seq(-4, 4, 0.2), main="rnorm", xlab="")はヒストグラムを-4から4の範囲で0.2間隔表記、タイトルはrnorm、X軸ラベルはなしで表すことを意味しています。
hist(a, prob=T)はY軸を頻度から確率密度への変更を示し、
density(a)は密度評価(density estimation)を意味し、
lines(density(a))は近似曲線(fitted curve)を描きます(近似曲線を描くにはY軸を確率密度にしなればなりません)。近似の様子を変えるには、
lines(density(a, bw=0.5), col="red")
lines(density(a, bw=0.1), col="blue")
などとbw(bandwidth)を変化させることにより可能です。
quantile(a)はクォンタイル、四分位を、
summary(a)はaの分布の様子の要約を示します。
quantile(a, prob=c(0.025, 0.975))はaの分布における上位2.5%および下位2.5%の値を示します。
なお、正規分布における理論値(theoretical value)はqnorm(c(0.025, 0.975))です。
qqnorm(a)はQ-Qプロットを示し、正規分布にaが従っていると直線になります。
分布に対しての最初の文字(例えば、normの最初についている文字)は、
dがdensity、pがprobability、qがquantile、rがrandomを意味しています。
boxplot(a)は箱ひげ図を示し、箱の下側が第一四分位、上側が第三四分位、中央は中央値を示します。
ひげの領域はdefaultでは箱の幅の1.5倍となっていますが、最大値や最小値がそれ以内の場合にはそこで切ります。
ひげの領域から外れた値を外れ値などと表現することがあります。
ひげの長さを変える際には、boxplot(a, range=2)などrangeを使います。
7. ブートストラップ法
リサンプリングによって母集団の分布を推定する方法の一つです。
例えば、a <- c(1, 10, 100, 1000, 10000)をオリジナルデータとします。
データ数が5であることだけがわかっているとし、5回の復元無作為抽出(random sampling with replacement)により一つの仮想データを作成します。この作業を1000回してみます。
データの平均値を知ることを目的とします。1000の仮想データをdに収納することを考えます。そこにはどのような数値が入るかわかりませんので、dには何もない状況で定義だけします。
d <- c()またはd <- numeric()とすれば、dは定義されたことになります。
for(i in 1:1000) d[i] <- mean(sample(a, 5, replace=T))により、dには1000の平均値が入ります。
hist(d)により、それらの様子を見てみましょう。
mean(d)とオリジナルデータの平均値mean(a)を比較してください。
8. 集合
a <- 1:50
b <- 10:60
union(a, b)はaとbの結合集合を意味し、
intersect(a, b)は共通集合を意味しています。
setdiff(a, b)は差集合を意味しています。
runif(50, 1, 100)は1から100の範囲で一様分布に従う乱数を50個生成しますが、それらは小数点以下を含みます。
整数部分だけにしたい場合には、trunc(runif(50, 1, 100))とします。
intersect(runif(50, 1, 100), runif(50, 1, 100))とintersect(trunc(runif(50, 1, 100)), trunc(runif(50, 1, 100)))を比べてください。
a <- trunc(runif(50, 1, 100))
runifやrnormはそれぞれの分布に従う無限の母集団から、ある個数のサンプリングを行うということです。
それを何度繰り返しても同じ要素からなるもにになりません。
ですので、
intersect(runif(50), rnif(50))
と
d1 <- runif(50); d2 <- runif(50); intersect(d1, d2) #異なることを実行させるときには、改行していますが、セミコロン（;）を使うこともできます。また、シャープ（#）を入れると、それ以降はRでは実行文以外と処理されますので、この部分をコピペしてもエラーは生じません。
上記２つのintersect( )は違っています。
前者は実行するごとにメンバーが違っていますが、後者はd1, d2で固定されたメンバーとなっているからです。
d1やd2は一様分布に従う無限の母集団からサンプリングされた世界で１つの標本ということになります。
前者を何度も繰り返して様子を見ることは意味がありますが、後者は繰り返しても無意味です。
さて、上記のaの場合に、
a[a>20]は20より大きな要素の抽出、
a[a>20 & a<=30]は20より大きく30以下の要素の抽出、
a[a<=20 | a>30]は20以下あるいは30より大きな要素の抽出を意味します。
a[1] <- 30はaの1番目を30に換えること、
rev(a)は反転(reverse)を、order(a)はaの中で小さい値から何番目に存在しているかを示します。defaultはdecreasing=Fです。
sort(a)は要素の値が小さいものから並び換えますので、a[order(a)]に同じです。
課題２ (1) a <- runif(1000) ; b <- runif(1000) とした場合、両者は一様分布に従う乱数が１０００個入っていますが、その値は違っているはずです。そこで、a - b を行うと、正の数と負の数はほぼ同じになり、完全に一致している場合はほぼゼロであるので、それぞれが５００に近いと考えられます。そのことを示してください。実行文とその結果を示してください。 (2) 次にaに対して1000回の復元無作為抽出を行いAに収納し、bに対して1000回の復元無作為抽出を行いBに収納し、(1)と同様にしてA - Bを行って、正と負の数を出してください。実行文とその結果を示してください。
9. 行列とデータフレーム
matrix(c(1, 2, 3, 4, 5, 6), 2, 3)は要素1, 2, 3, 4, 5, 6を2行3列に行列させます。
matrix(c(1, 2, 3, 4, 5, 6), 2)でも同じ結果です。
matrix(c(1, 2, 3, 4, 5, 6), nrow=2)という意味で、matrix(c(1, 2, 3, 4, 5, 6), ncol=3)でも同じです。
matrix(c(1, 2, 3, 4, 5, 6), 2, byrow=T)は要素を並べる際に行を埋めて次の行へということを示します。
a <- c(1, 2, 3, 4, 5)
b <- c(2, 4, 6, 8, 10)
ab <- data.frame(a, b)
はベクトルaとbをデータフレーム化、aとbは列となります。
ab[ ,1]はデータフレームabの1列目を抽出、
ab[ ,2]は2列目を抽出、
ab[1, ]は1行目を抽出、
ab[3, ]は3行目を抽出することを意味しています。
a <- matrix(1:6, 2, byrow=T)
rownames(a) <- c("A", "B")は1行目をA、2行目をBと名付けます。
colnames(a) <- letters[1:3]は1列目、2列目、3列目をa, b, cと名付けます。
t(a)は転置関数(transposition function)を意味します。
cbind(A=1:3, B=4:6)は列で結合、
rbind(A=1:3, B=4:6)は行で結合します。
a <- matrix(1:10, ncol=2)
a[ ,1]は1列目を抽出、
a[ ,1]>3は1列目の要素で3より大きいもの、
a[a[ ,1]>3, ]は1列目の要素で3より大きいものを持つ行を抽出、
a[a[ ,1]>3, ][ ,2]は1列目の要素で3より大きいものを持つ行を抽出し、その2列目の要素を抽出、a[a[,1]>3, 2]でも同じです。
apply(a, 1, mean)は各行の平均、apply(a, 2, median)は各列の中央値、apply(a, c(1,2), sqrt)はすべての要素の平方根を示し、sqrt(a)と同じです。
行列やデータフレームの大きさは
dim(a)で行数、列数を知ることができます。
行数だけを知りたい場合には
nrow(a)
列数だけを知りたい場合には
ncol(a)とします。
例えば、
d1 <- runif(100)
d2 <- runif(100)
D <- data.frame(d1, d2)
において、同じ行における1列目が2列目よりも小さい値である行の数を求めるには、
nrow(D[D[,1] < D[,2], ])
整数部分だけを示すことと行列を使って次のようなことができます。
例えば、白球、赤球、青球がそれぞれ100個あり、それらを無作為に1列に並べるとします。
その際に白球と白球の間に赤および青球が存在する場合の数をRで計算する場合、
a <- sample(c(rep(0,100), rep(1,100), rep(0.001,100)), 300)
A <- data.frame(1:300, a)
b <- A[A[,2]==0, 1]
d <- c()
for(i in 1:99) d[i] <- sum(A[b[i]:b[i+1], 2])
length(d[d>1 & d>trunc(d)])
また、乱数を発生させた際には、小数点以下が入りますので、それを整数として扱う場合があります。
例としてショットガンシークエンスについて述べます。
ショットガンシークエンスとは、DNAシークエンスにおいてランダムに塩基配列を決定して、それらをつなぎ合わせる方法です。
大量並列型DNAシークエンスはこの方法です。
例えば、4M塩基対からなるバクテリアの塩基配列を決めるとします。
一回のシークエンスは200塩基とすると、過不足なく完全に効率的にシークエンスしたとすると20000回で終了しますが、実際には重複する領域がありますので、20000回シークエンスした際の断片の末端の位置をシミュレーションすると、
a <- trunc(runif(20000, 1, 4*10^6))
これらの位置間の距離がすべて200塩基以下の場合に全長がシークエンスできたことになります。
b <- sort(a)
d <- c()
for(i in 1:19999) d[i] <- b[i+1] - b[i]
max(d)
length(d[d > 200])によってギャップの数がわかりますので、それにプラス1した数がコンティグ数ということになります。
当然、20000回のシークエンスで200塩基を下回ることはありません。
では何回シークエンスすれば200塩基以下になるか試してみてください。
課題３ 上記のシークエンスにおいて２０万回行った場合にはコンティグ数がいくつになるか答えてください。その際、Rでの実行文も示してください。実行文を示さない場合には採点対象外となります。
DNAシークエンスということが出てきましたので、ここでさらにDNAの塩基組成（GC含量）の計算について述べます。
DNAは二重らせん構造をとっていますが、二つの鎖はアデニン（A）とチミン（T）およびグアニン（G）とシトシン（C）がそれぞれ水素結合してニ重鎖となっています。
よって、AとTの数は等しく、GとCの数も等しくなっています。
このことはワトソンとクリックによるDNAのニ重らせん構造の発表より以前にシャルガフ（Erwin Chargaff）によって示されています。
ゲノムDNAを比較する際（特にバクテリアのゲノムDNAの場合）、全塩基に占めるGとCの割合が使われることがあり、GC含量と呼ばれています。
もちろん、GC含量とAT含量の和は1（100%）となっています。
それでは実際にゲノムのGC含量を調べてみましょう。
NCBIのトップページにアクセスしてください。
All Databasesに対してMixia osmundaeと打ち込み、Searchをクリックしてください。
そうすると、これまでに異なる２つのグループがMixia osmundaeのゲノム決定を行ったことがわかります。
ここで使用するのはMosmundae 2.0というものです。
このサイトにたどり着いた場合には、WGS projectのBABT02をクリックしてください。
BABT02000001-BABT02000283をクリックすると下のような画面になります。
ここより、283のコンティグ塩基配列にアクセスすることができます。
例えば、BABT02000007は286369塩基、128のタンパク質コード領域を持っていることがわかります。
Viewのカラムに表記されているGenBankをクリックしてください。
どのようなタンパク質がどの領域にコードされているかがわかります。
次に塩基配列情報をダウンロードします。
右上にあるsendをクリックし、FASTAを選択します。
上記のサイトにたどり着きましたか？ほかのグループもゲノムを決定したので、データベースには２つあり、混乱したひとがいるかもしれません。
このサイトは
https://www.ncbi.nlm.nih.gov/nuccore/BABT02000007.1
です。たどり着けなかったひとは、ここにアクセスしてください。
このファイルをダウンロードすると、sequence.fastaというファイルができます。
このファイルをメモ帳などで開き、FATSA形式の塩基配列であることを確認してください。
ここではGC含量のデータ解析をしますので、FASTAファイルにおける塩基配列以外の情報はすべて削除してください。
次に、AとTに対して0 [数字の0とスペース]に置換し、CとGに対して1 [数字の1とスペース]に置換してください。
ファイルが以下のようになれば、それをgc.txtとしてフォルダgenomeに保存してください。
次に遺伝子の位置情報を得るためにBABT02000007のサイトに戻ってください。
先ほどFASTAを選択したところにおいて、Feature Tableを選択してください。
ファイルをダウンロードすると、sequenceというファイルができます。
このファイルをExcelで開いてください。
各列においてフィルター設定し、列Cにおいて「gene」と表記されている行だけを表記してください。
そこに表記されているものがこのコンティグに載っている遺伝子の位置情報です。
2列からなる位置情報をコピーしてメモ帳にペーストし、gene.txtとしてフォルダgenomeに保存してください。
Rにおいてgc.txtおよびgene.txtを読み込みます。
ファイルの読み込みのためには、Rを実行しているフォルダ内にファイルがあることが必要です。
あるいは、ファイルの存在している場所にディレクトリを変更する必要があります。
たとえば、デスクトップに置いている場合には、Rの左端にあるファイルをクリックして、ディレクトリの変更を選択し、デスクトップに変更してください。
gc <- scan("gc.txt")
gene <- read.table("gene.txt")
Webサイト上で行っているひとは、上記のことができません。
工夫すれば、geneは読み込めますが、gcは容量オーバーでできません。
各遺伝子の長さのデータを取ります。
size <- c( )
for(i in 1:nrow(gene)) size[i] <- length(gc[gene[i,1]:gene[i,2]])
次に各遺伝子におけるGとCの数のデータを取ります。
GC <- c( )
for(i in 1:nrow(gene)) GC[i] <- sum(gc[gene[i,1]:gene[i,2]])
よって、各遺伝子のGC含量の分布は
hist(GC/size)
boxplot(GC/size)
この分布が正規分布に近いかどうかを知るためには
qqnorm(GC/size)
平均値は
mean(GC/size)
中央値は
median(GC/size)
課題４ BABT02000007と同様にしてBABT02000004のmean(GC/size)を計算してください。 なお、BABT02000004は下記からアクセスしてください。
https://www.ncbi.nlm.nih.gov/nuccore/BABT02000004.1
10. 頑強性
実験データにノイズが入ることは一般的です。よって、それらをいかに取り除く、あるいは影響を最小限にとどめるための工夫を行わなければなりません。
例えば、同じ実験を何回か行った際にその代表値(representative value)として相応しいものは何か考えてみましょう。
最大値や最小値を代表値とすることは抵抗を感じる場合が多いと思いますが、平均値と中央値ではどうでしょうか？
例えば、10回の繰り返し実験と行ったとします。その実験値がa <- runif(10)(0と1の間の一様分布に従う10つの乱数)としますと、その平均値はmean(a)であり、中央値はmedian(a)となります。
仮に、10回の中で10というノイズが一度観測された場合には、b <- c(runif(9), 10)となり、mean(b)、median(b)となります。
中央値の方が平均値よりもノイズによる影響が小さいことがわかります。この場合、中央値の方が平均値よりも頑強性(robustness)があるといいます。
実験データ間における違い(実際には実験データを与える母集団間における違い)があるかどうかを調べる際（検定）、t検定などのパラメトリックな方法とウィルコクソン順位和検定(マン・ホイットニーU検定)などのノンパラメトリックな方法があります。Rでは様々な検定を行うコマンドがありますので、この講義でも説明いたします。
実験データに基づき、コントロールとトリートメントが有意に(significantly)に違うか、有意差はないかが結果として記載される必要がありますが、その場合には検定を行う必要があります。
すなわち、生物情報学の最も基礎にある領域は統計学(statistics)です。
11. t検定とウィルコクソン順位和検定
ここに2のサンプルがあったとします。これらのサンプル間における検定を行うとします。
検定には帰無仮説(null hypothesis)と対立仮説(alternative hypothesis)があり、帰無仮説はサンプルの母集団に差がない(同じ母集団)、対立仮説は差があるということになります。
私たちに見えているものは、サンプル（標本）であり、母集団ではありませんが、検定するのは、母集団が同じかどうかという視点で行うということです、
t検定ではt、ウィルコクソン順位和検定(マン・ホイットニーU検定)ではUという検定統計量（下記の例において確認してください）を算出します。
その算出した検定統計量が起こり得るすべての場合における検定統計量の分布においてどの位置にあるかを知ります。
その位置が起こり得る確率が低い方から積算してどれほどの確率で生じるかを計算し、その総和(p値)が0.05よりも小さい場合には帰無仮説を棄却し、0.05以上の場合には棄却しません。
Wilcoxonは1945年にウィルコクソン順位和検定を発表しました(Biometrics Bulletin 1, 80-83)。
MannとWhitneyは1947年にU検定を発表しました(Annals of Mathematical Statistics 18, 50-60)。
のちにこれら2つの検定は同じことがわかりました。その方法を示します。
例えば、サンプル甲のデータが1, 2, 4, 7, 9であり、サンプル乙のデータが3, 5, 6, 8であるとき、これらの母集団に有意な差があるかどうかを検定します。
帰無仮説は「サンプル甲および乙の母集団に差がない」であり、対立仮説は「サンプル甲および乙の母集団に差がある」です。
サンプル甲と乙の間で各々の数の大小を比べます（例えば、サンプル甲の1はサンプル乙の4つの数のすべてより小です）。
いま、5×4=20の比較を行い、サンプル甲が大の場合の数は8、小の場合の数は12となります。
そこで、検定統計量Uを「大、小の場合の数の少ない方とする」と定義します（サンプル間において偏りが大きいほど、Uは0に近くなり、偏りが少ないほど10に近づくことがイメージできます）。
この場合、U=8となります。
さて、帰無仮説に基づくすべてのパターンを考えますと、1から9の数字を5個取り出す（あるいは4個取り出す）場合の数、すなわち、126通りあります。
帰無仮説のもと、これら126通りは等確率で生じます。これらをUの値で整理すると、
この表において、U=0とはどちらかのサンプルがすべて大あるいは小となる場合となり、母集団が同じの場合では、最も起こり難い場合になります。
そこで、U=8が生じる起こりやすさを考えますと、U=0からU=8までの確率を足し合わせることにし、これがp値となります。
よって、p値=0.016+0.016+0.032+0.048+0.079+0.095+0.127+0.142+0.175=0.73
有意水準（危険率）が0.05の場合、U=0またはU=1の場合に、棄却域に入り、帰無仮説を棄却することとなります。
U=8の場合は、p > 0.05ですので、帰無仮説は棄却されません。
よって、サンプル甲と乙の母集団に差がないことを棄却できないということになります（差があると言ってはいけないということ）。
課題５ グループ１のサンプル数が１０、順位和が１１３であり、グループ２のサンプル数が１５、順位和が２１２の場合、検定統計量Uを求めてください。 なお、U1 = R1 - (N1^2 + N1)/2, U2 = R2 - (N2^2 + N2)/2 であり、U1, U2の数の少ない方をUとする。RおよびNは順位和およびサンプル数を示す。
次にt検定について述べます。
[0, 1]間の一様分布を母集団として、10のサンプルサイズの2つの標本間におけるt検定を行います。
t検定は1908年、Gosset(発表はペンネームStudentとして)によりBiometrika (vol. 6, pp. 1-25)に発表されました。
Rでは、2つの標本の分散が等しい場合と等しくない場合によって、異なりますので、注意してください。
defaultの設定では、2つの分散が等しくないとして、Welch 2標本検定が行われます。
a <- runif(10); b <- runif(10) #runif(10)は一様分布に従う無限の母集団から10個をサンプリングすることです。よって、aとbにおいて一致している要素は一つもありません。
t.test(a, b)と打ち込んで確認してください。
分散が等しいかどうかについては、F検定を行います。Rでは、var.test(a, b)としてください。
F値はvar(a)/var(b)のことであり、自由度(9, 9)のF分布（分布の様子は２つの自由度によって決まります、ここではサンプル数が10と10なので、自由度はそれぞれ9となります）において、この値の位置を知ることができればよいわけです。
自由度(9, 9)のF分布の様子を見てみましょう。
x <- seq(0, 5, 0.01)
plot(x, df(x, 9, 9))
var(a)/var(b)のときもvar(b)/var(a)のときもp値は同じになります（もちろん、自由度が違っている場合には、一方が(n, m)の分布であれが、他方は(m, n)の分布においてです）。
var.test(a, b)とvar.test(b, a)で確認してください。
有意水準(危険率、significance level)を5%とした場合、F検定において、p値が0.05よりも小さい場合には帰無仮説（2つの標本の母集団に違いはない）を棄却します。すなわち、分散は等しくないとして、t検定を行います。
F検定の結果、p値が0.05より大きい場合には帰無仮説を棄却できませんので、等分散として、t検定を行います。すなわち、t.test(a, b, var.equal=T)とします。
等分散ではない場合と結果を比べてください。
乱数を使っていますので、それぞれのコンピュータ上で異なる結果となっているかと思いますが、95%の確率でp値が0.05より大きくなり、帰無仮説(2つの標本の母集団に違いはない)を棄却できない結果のはずです。
有意水準5%とは、同じ母集団であっても5%の確率で違うと判断してしまうことを意味していますので当然ですね。
そのことを確認するために、同じt検定を1000回行い、p値が0.05より小さい場合の数を数えてみます。
d <- c()
for(i in 1:1000) d[i] <- t.test(runif(10), runif(10))$p.value #この実行によって、無限の母集団からの10個のサンプリングを２回して、それらを比べるということを1000回します。どの組みあわせにおいても、一つの要素も一致しているものはありません（正確に言えば、天文学的な確率の低さです）。しかし、母集団は同じであるというところがミソです。もう一つ、$を使うと、結果から特定のものだけを抽出できます。ここでは、p値だけを抽出しています。
length(d[d<0.05]) #p値が0.05より小さい場合の数を数えます
同様なチェックをウィルコクソン順位和検定（マン・ホイットニーU検定）でもしてみましょう。
ノンパラメトリックな方法ですので、等分散等の確認をする必要はありません。
wilcox.test(runif(10), runif(10))とすれば結構です。この場合においてもp値の確認をしてみましょう。
d <- c()
for(i in 1:1000) d[i] <- wilcox.test(runif(10), runif(10))$p.value
length(d[d<0.05])
それでは異なる母集団からの標本の場合はどうでしょうか？
runif(10)とrnorm(10)を検定してみましょう。
dT <- c()
for(i in 1:1000) dT[i] <- t.test(runif(10), rnorm(10))$p.value
length(dT[dT<0.05])
dW <- c()
for(i in 1:1000) dW[i] <- wilcox.test(runif(10), rnorm(10))$p.value
length(dW[dW<0.05])
棄却されなければならないはずですが、結構棄却されていないことがわかります。
課題６ 上記で実行したことは、dTがrunif(10)とrnorm(10)におけるt検定を１０００回行った際の５％棄却域に入った数を示し、dWがウィルコクソン順位和検定での数です。すなわち、サンプルサイズは１０ということです。そこで、サンプルサイズを１００にした場合に同様に１０００回の検定を行った際、その際の５％棄却域に入った数を産出し、考察してください。その際の実行文とその結果を示し、その上で考察してください。
有意水準(危険率)が真の仮説を棄却する確率を示し、その値を小さくすればするほど、偽の仮説が棄却されない場合が増します。検出力(statistical power)とは偽の仮説が棄却される確率を言います。
t検定のようにパラメトリックな検定においては、他の条件が定まれば、あと一つの条件を求めることができます。
例えば、検出力0.8、有意水準0.05、標準偏差2の分布において0.5の差を検定するために必要な標本の大きさは
power.t.test(power=0.8, sig.level=0.05, sd=2, delta=0.5)により253が求められ、
標本の大きさ 200、有意水準0.05、標準偏差2の分布において0.5の差を検定するときの検出力は
power.t.test(n=200, sig.level=0.05, sd=2, delta=0.5)により0.703が求められます。
12. 累積相対度分布による比較
よく目にする分布は大半が相対度分布です。
例えば、釣鐘型になる正規分布などを思い描いてください。
この相対度をすべて累積していくと最終的には100%となりますので、累積相対度分布は0%から単調に増加して100%になります。
経験累積度数関数(empirical cumulative distribution function)を使ってa <- runif(50)の累積度数分布は
plot(ecdf(a))
b <- rnorm(50, 0.5)の累積度数分布は
plot(ecdf(b))
で見ることができます。
そこで異なる分布を累積相対度分布で比較することをコルモゴロフとスミルノフが発表しました(コルモゴロフ・スミルノフ検定)。
例えば、aとbの累積相対度分布で比較するには、
ks.test(a, b)
この結果と、t.test(a, b)やwilcox.test(a, b)の結果と比べて見てください。
13. 分散分析
標本が３つ以上ある場合にはどうすればよいでしょうか？ 3標本以上が同じ母集団を持つかどうかを検定するにはどのようにすればよいでしょうか？
その場合には最初に分散分析(analysis of variance、ANOVA)を行います。
例えば、グループA、B、Cの要素がc(2, 4, 9)、c(1, 3, 5, 7)、c(2, 3, 4)の場合、それぞれのグループにおける平均値は5, 4, 3であり、全体では4です。
総平方和は
(2-4)2+(4-4)2+(9-4)2+(1-4)2+(3-4)2+(5-4)2+(7-4)2+(2-4)2+(3-4)2+(4-4)2=54
グループ間平方和は
(5-4)2+(5-4)2+(5-4)2+(4-4)2+(4-4)2+(4-4)2+(4-4)2+(3-4)2+(3-4)2+(3-4)2=6
グループ内平方和は
(2-5)2+(4-5)2+(9-5)2+(1-4)2+(3-4)2+(5-4)2+(7-4)2+(2-3)2+(3-3)2+(4-3)2=48
となります。すなわち、総平方和＝グループ間平方和＋グループ内平方和となっています。
この場合の分散分析表を示します。
すなわち、自由度(2, 7)のF分布におけるF値0.44の位置を知ればよいということになります。
d1 <- c(2, 4, 9)
d2 <- c(1, 3, 5, 7)
d3 <- c(2, 3, 4)
data <- c(d1, d2, d3)
group <- factor(rep(1:3, c(3, 4, 3)))
分散分析表はanova(lm(data~group))により見ることができます。
Prのところの数値がP値を示しています。
t検定のときと同様に、デフォルトでは等分散のサンプルとして扱われていないため、等分散として扱ってよいかどうかを検定します。
bartlett.test(data~group)を行います。
等分散として扱える場合にはoneway.test(data~group, var.equal=T)、扱えない場合にはoneway.test(data~group)
データの順位に基づく検定として、クラスカル・ウォリス(Kruskal-Wallis)検定があり、Rではkruskal.test(data~group)となります。
Kruskal WH, Wallis WA (1952) Use of ranks in one-criterion variance analysis. Journal of the American Statistical Association 47, 583-621.
どの群間に差があるかを調べる場合
a <- c(1, 2, 3, 4, 5, 6, 7, 8, 10)
b <- factor(c(1, 1, 1, 2, 2, 2, 3, 3, 3))
bartlett.test(a~b)
oneway.test(a~b, var.equal=T)
pairwise.t.test(a, b)
pairwise.t.test(a, b, pool.sd=F)は共通の標準偏差を用いない場合
pairwise.t.test(a, b, pool.sd=F, p.adj="bonf")Holm補正(デフォルト)からBonferroni補正へ
t.test(c(1, 2, 3), c(7, 8, 10))$p.value
t.test(c(4, 5, 6), c(7, 8, 10))$p.value
t.test(c(1, 2, 3), c(4, 5, 6))$p.value
との比較からHolmおよびBonferroni補正を確認できます。
s1 <- runif(10); s2 <- runif(20); s3 <- rnorm(10)の3標本の分散分析を行いましょう。
全体の平均値はmean(c(s1, s2, s3))ですので、グループ間平方和は10*(mean(s1)-mean(c(s1, s2, s3)))^2+20*(mean(s2)-mean(c(s1, s2, s3)))^2+10*(mean(s3)-mean(c(s1, s2, s3)))^2であり、グループ内平方和はsum((s1-mean(s1))^2)+sum((s2-mean(s2))^2)+sum((s3-mean(s3))^2)です。
自由度はグループ間では3-1=2、グループ内（誤差）では10-1+20-1+10-1=37ですので、平均平方はグループ間が(10*(mean(s1)-mean(c(s1, s2, s3)))^2+20*(mean(s2)-mean(c(s1, s2, s3)))^2+10*(mean(s3)-mean(c(s1, s2, s3)))^2)/2、グループ内が(sum((s1-mean(s1))^2)+sum((s2-mean(s2))^2)+sum((s3-mean(s3))^2))/37となります。
よって、F値は((10*(mean(s1)-mean(c(s1, s2, s3)))^2+20*(mean(s2)-mean(c(s1, s2, s3)))^2+10*(mean(s3)-mean(c(s1, s2, s3)))^2)/2)/((sum((s1-mean(s1))^2)+sum((s2-mean(s2))^2)+sum((s3-mean(s3))^2))/37)となります。
Rにおいて、一元配置分散分析を行うには、data <- c(s1, s2, s3); group <- factor(rep(1:3, c(10, 20, 10)))として、anova(lm(data~group))で分散分析表の表示できますので、先ほどのF値と一致しているかを確認してください。
oneway.test(data~group)としても結果を見ることができます。この場合、F値が異なっていることを確認してください。これはt.test()のときと同様で、3標本の分散が等しくないことがdefaultとなっているためです。
oneway.test(data~group, var.equal=T)として比較してください。
二元配置分散分析ではデータに2つのファクターがかかわります。
例えば、染色体における遺伝子の数が染色体の種類とワトソン鎖・クリック鎖のDNA構造上の場所どちらに関連しているかどうかを調べるような場合です。酵母Saccharomyces cerevisiaeの16本の染色体のワトソン鎖とクリック鎖にはそれぞれ下記の数の遺伝子がコードされています。
watson <- c(51, 187, 68, 371, 137, 64, 275, 153, 94, 183, 159, 244, 240, 206, 278, 244)
click <- c(46, 219, 90, 377, 136, 62, 250, 128, 111, 173, 151, 268, 221, 187, 258, 216)
この場合、染色体の種類とDNAの構造の2つのファクターがありますので、
data <- c(watson, click)
wc <- factor(c(rep(1, 16), rep(2, 16))) #ワトソン鎖、クリック鎖の場合分け
chr <- factor(rep(1:16, 2)) #染色体ごとの場合分け
anova(lm(data~wc+chr))
となり、染色体の種類において偏りがあることがわかります（Prの列に書いている数値がp値を示します）。
課題７ 上記では染色体毎における違いに有意差が認められました。しかし、それは染色体毎に長さが違っているからではないかと考えられるので、各染色体における遺伝子の密度を計算し、それに対して同様の分散分析をしてください。なお、各染色体のサイズ（Kbp）は
size <- c(230, 813, 317, 1532, 577, 270, 1091, 563, 440, 746, 667, 1078, 924, 784, 1091, 948)
です。上記の遺伝子数をどのようにすれば遺伝子密度になるか考えて、実行文とその結果を示してください。
14. 対応のあるデータ解析
次に対応のあるデータの場合について考えてみます。
遺伝子1から遺伝子20までの遺伝子発現レベルがコントロール条件でa <- runif(20)、トリートメント条件でb <- runif(20)であった場合、発現レベルに違いがあるかどうかを検定するにはt.test(a, b, var.equal=T)やwilcox.test(a, b)でよいでしょうか？
この場合、遺伝子1のデータ値としてa[1]とb[1]、遺伝子2のデータ値としてa[2]とb[2]・・・とそれぞれ対応しています。
よって、t.test(a, b, paired=T)、wilcox.test(a, b, paired=T)とする必要があります。
独立したデータと対応のあるデータでは、t検定の際の自由度が半分になっていることを確認し、p値の結果が違うことも確認してください。
どうして違うか考察してください。
これらの結果はt.test(a-b)およびwilcox.test(a-b)に等しいことがわかります。
もし、aおよびbがランダム（無関係）であれば、どちらが大きいかは半々ですので、a-bが0に対して大きくなる場合と小さくなる場合は半々に近いはずです。
しかし、トリートメントにより、遺伝子発現レベルが上がったり下がったりしている場合にはa-bは正または負に偏ります。
上記の検定はそのようなことを検定しているといえます。
15. 回帰と相関
ある処理を行うことで遺伝子発現レベルが上がったとしても、一般的にはそのレベルには大小があります。
遺伝子発現レベルに差があった遺伝子がどれほど同じようにに変化したどうかを知るにはどのようにすれば良いでしょうか？
例えば、ここに異なる20の遺伝子における遺伝子発現を比較する実験を行ったとし、処理前の遺伝子発現レベルをd1、処理後をd2とします。
d1 <- runif(20)およびd2 <- runif(20, 10, 20)の場合、すべての遺伝子発現レベルは上昇しましたが、同じ程度ではありません。
D1 <- sort(d1)およびD2 <- sort(d2)の場合、すべての遺伝子発現レベルは上昇し、それらは同じ程度です。
ここで注意してほしいことは、d1, d2, D1, D2はそれぞれサンプルとしての例です。
sort( )というステップを踏むというために書いているのではありません。
それぞれ違ったサンプルとして作成しただけです。
最初に、t.testおよびwilcox.testではd1とd2およびD1とD2の間で同じような結果となることを確認してください。
次に、plot(d1, d2)およびplot(D1, D2)を実行して、散布図を見てください。
今、コントロールの遺伝子発現レベル(x)が与えられたときにトリートメント後の遺伝子発現レベル(y)を予測できるかを考えてみます。
yとxの関係が直線上に乗ると仮定した場合、y = a + bxと書くことができます。
しかし、実際にはこの直線より外れたところに実測値が存在していますので、この実測値とこの直線から予測された数値との差が存在しています。
この差を最小にする、すなわち、残差平方和(sum of residual squared)= ∑(y[i] – a – bx[i])2を最小にするa、bを求め、このa、bをy切片、傾きとする直線を回帰直線(linear regression line)と呼びます。
参考までに
∑(y[i] - a - bx[i])2
= ∑y[i]2 + ∑a2 + ∑b2x[i]2 - 2∑ay[i] + 2∑abx[i] - 2∑bx[i]y[i]
= ∑y[i]2 + na2 + b2∑x[i]2 - 2a∑y[i] + 2ab∑x[i] - 2b∑x[i]y[i]
aで偏微分して
2na - 2∑y[i] + 2b∑x[i] = 0
よって
a = (∑y[i] - b∑x[i])/n
同様にbで偏微分して
2b∑x[i]2 + 2a∑x[i] - 2∑x[i]y[i] = 0
よって
b∑x[i]2 + (∑x[i]∑y[i] - b∑x[i]∑x[i])/n – ∑x[i]y[i] = 0
b(∑x[i]2 – ∑x[i]∑x[i]/n) + ∑x[i]∑y[i]/n - ∑x[i]y[i] = 0
以上より
b = (∑x[i]y[i] – ∑x[i]y[i]/n)/(∑x[i]2 – ∑x[i]∑x[i]/n)
a = ∑y[i]/n – (Sxy/Sxx)∑x[i]/n
Sxy = ∑x[i]y[i] – ∑x[i]y[i]/n
Sxx = ∑x[i]2 – ∑x[i]∑x[i]/n
遺伝子発現が上昇した先の例の場合、lm(d2~d1)およびlm(D2~D1)により、y切片および傾きを得ることができます。
また、plot(d1, d2)後にabline(lm(d2~d1))とすることにより、回帰直線を散布図上に書き込むことができます。
d1とd2における回帰直線と、plot(D1, D2)後のabline(lm(D2~D1))における回帰直線を比較すると、モデル化は後者の方が優れていることがわかるはずです。
summary(lm(d2~d1))およびsummary(lm(D2~D1))におけるy切片および傾きの標準誤差やp値を比較して確認してください。
次に、xとyの相関を見ましょう。Rではcor()によりピアソンの積率相関係数（r = ∑((x[i]-xの平均)(y[i]-yの平均))/sqrt(∑(x[i]-xの平均)2・∑(y[i]-yの平均)2)）を求めることができます。
cor(d1, d2)とcor(D1, D2)を比べてください。
なお、cor(d1, d2)*sd(d2)/sd(d1)およびcor(D1, D2)*sd(D2)/sd(D1)がそれぞれの回帰直線の傾きに等しいことを確認してください。
スピアマンの順位相関係数（ピアソンの積率相関係数のx[i], y[i]がそれぞれにおける順位）の場合は、cor(d1, d2, method="s")です。
ケンドールの順位相関係数はxy上の点1と点2における(x1 - x2)×(y1 - y2)の値が正の場合に1、負の場合に-1として、n(n-1)/2のすべての2点間の関係を調べ、総和をとります。
この総和をn(n-1)/2で割ると、そのときの値は-1以上1以下となり、正の相関がある場合には1、負の相関がある場合には-1に近づきます。
Rでは、cor(d1, d2, method="k")です。
相関係数に関する検定
帰無仮説: 相関なし(相関係数が0)
ピアソンの積率相関係数の有意性検定
cor.test(x, y, method="pearson")
method="p"でも可、default
検定統計量(statistic): r×sqrt(n-2)/sqrt(1-r2)
自由度(degree of freedom)n-2のt分布
スピアマンの順位相関係数の有意性検定
cor.test(x, y, method="spearman")
method="s"でも可
ケンドールの順位相関係数の有意性検定
cor.test(x, y, method="kendall")
method="k"でも可
課題８ xy平面上に４つの点があるとします。それらの相関係数がゼロとなるような（その際、回帰直線の傾きがゼロになっています）４つの点の具体的な例を示してください。
回帰直線は主成分分析(principal component analysis)における第一成分とは異なります。主成分分析では各点と直線の距離を考えていますので、当然、yの値で考えている回帰直線とは異なります。
参考までに主成分分析を上記の2つ場合で行った結果は、
summary(prcomp(data.frame(d1, d2)))
summary(prcomp(data.frame(D1, D2)))
により見ることができます。
標準偏差、寄与率、累積寄与率が表示されます。
プロットする場合には、
biplot(prcomp(data.frame(d1, d2)))
biplot(prcomp(data.frame(D1, D2)))
としてください。
16. クラスタリング
次に系統樹の作成をRをつかってやってみましょう。距離行列に基づく階層的クラスタリングを行ってみます。
d1 <- runif(20)
d2 <- runif(20)
を使って説明します。
散布図とクラスタリング結果を見比べるため、
par(mfcol=c(1,2))として、
plot(d1, d2, type="n")
text(d1, d2) をしてください。
次に、各点から各点への距離の算出は
dist(data.frame(d1, d2))
で行います。
この距離行列に基づき、階層的クラスタリング(hierarchical clustering)を行います。
plot(hclust(dist(data.frame(d1, d2))))
と打ち込んでください。
階層的クラスタリングの具体的な例として、平均連結法について説明します。
下記のような距離行列があったとします。
A B C
B 3
C 2 1
D 4 4 7
距離が最も短いBとCが最初にクラスターを形成します。これらは、共通祖先より0.5の距離をもって位置しているとします。
そうすると次の距離行列は次のようになります。
A BC
BC 2.5
D 4 5.5
BCとの距離が平均距離として出しています。この中では、AとBCが最も短いので、クラスターを形成します。それぞれの地点までの距離は1.25です。
そうすると次の距離行列は次のようになります。
ABC
D 5
最後にABCとDが共通祖先より2.5の距離で連結して、クラスタリングが完成します。
このような方法は、類似した塩基配列やアミノ酸配列を比較して、系統樹を書く際に使用されます。
距離行列法といわれる方法ですが、興味のあるひとは、近隣結合法などを調べてみてください。
課題９ 階層的クラスタリングを行った際、複数の系統樹が作成される場合があります。それはどのような時か、平均連結法を例にして答えてください。
この階層的クラスタリングに対して、非階層的クラスタリングの例として、K-means(K平均)法を用いてクラスタリングを行ってみましょう。
この方法ではクラスター数を予め指定します。例えば、5としましょう。すると、kmeans(as.matrix(data.frame(d1, d2)), 5)とすれば結構です。
散布図において分かり易く見てみましょう。
plot(data.frame(d1, d2), col=kmeans(as.matrix(data.frame(d1, d2)), 5)$cluster)
とすれば、5つのグループが色分けされます。
この操作を何回か繰り返すと、異なるクラスタリングを行うことが確認できます。
「ゲノム比較解析」の「４．ゲノムシグネチャ」および「５．遺伝子アノテーションに基づく系統解析」へ進んでください。
17. 分割表を用いた検定
次に分割表を用いた検定について述べます。
例）DNA塩基配列比較のところで述べましたが、アミノ酸をコードしている領域においてはアミノ酸を変えない同義的サイトとアミノ酸を変える非同義的サイトがあります。
実際に生じた同義的塩基置換と非同義的塩基置換の頻度は同じかどうかを検定する場合、分割表を使います。
同義的サイト数がS、非同義的サイト数がNとし、実際に塩基置換が生じた数が同義的サイトにおいてs、非同義的サイトでnであったとしますと、次の分割表となります。
例えば、S=193、N=710、s=99、n=387の場合、
フィッシャー正確確率検定
帰無仮説: 同義的・非同義的塩基置換の頻度は同じ
対立仮説: 同義的・非同義的塩基置換の頻度は同じではない
帰無仮説のもと、この分割表が得られる確率を求めます。
合計: 903を193と710に分ける場合の数で903!/193!/710!
塩基置換: 486を99と387に分ける場合の数で 486!/99!/387!
非塩基置換: 417を94と323に分ける場合の数で 417!/94!/323!
よって(486!/99!/387!)×(417!/94!/323!)/(903!/193!/710!)
さらに、これより極端なものを全て計算し、それらを足し合わせます。
fisher.test(matrix(c(99, 94, 387, 323), 2, byrow=T))
カイ2乗検定
(要素 - 期待値)2/期待値 = χ2
帰無仮説に基づくと、塩基置換サイトにおいては
同義的サイト数 486×193/903、非同義的サイト数 486×710/903
非塩基置換サイトにおいては
同義的サイト数 417×193/903、非同義的サイト数 417×710/903
が期待値です。
これらの期待値に基づき、それぞれのχ2を求めることができ、それらの和0.6298を得ることができます。
この場合は自由度1((2-1)*(2-1))のχ2分布を適用します。
chisq.test(matrix(c(99, 94, 387, 323), 2, byrow=T))
chisq.test(matrix(c(99, 94, 387, 323), 2, byrow=T), correct=F)
連続性の補正(修正)なし、上記のχ2の値はこの場合です。
なお、2×2の分割表の場合にはmatrix(c(99, 94, 387, 323), 2)としても横に同義的サイト数と非同義的サイト数、縦に塩基置換数と非塩基置換数を書いたことになりますので、検定では同じ結果になります。
fisher.test(matrix(c(99, 94, 387, 323),2))
chisq.test(matrix(c(99, 94, 387, 323),2))
1000行2列からなるデータにおいて、1列と2列の値が有意に違っている行を抜き出すことをやってみましょう。
a <- trunc(runif(1000, 1, 100))
b <- trunc(runif(1000, 1, 1000))
ab <- data.frame(a, b)
d <- c()
for(i in 1:(nrow(ab)-1)) d[i] <- fisher.test(matrix(c(ab[i,1], ab[i,2], sum(ab[c(1:(i-1), (i+1):nrow(ab)), 1]), sum(ab[c(1:(i-1), (i+1):nrow(ab)), 2])), 2))$p.value
d1000 <- fisher.test(matrix(c(ab[1000,1], ab[1000,2], sum(ab[1:999, 1]), sum(ab[1:999, 2])), 2))$p.value
abd <- data.frame(a, b, c(d, d1000))
abd[abd[,3] < 0.05, ]
上記は場合分けして面倒ですね。下記でもできます。
for(i in 1:nrow(ab)) d [i] <- fisher.test(matrix(c(ab[i,1], ab[i,2], sum(ab[-i,1]), sum(ab[-i,2])),2))$p.value
abd <- data.frame(a, b, d)
abd[abd[,3] < 0.05, ]
課題１２ 上記はfisher.test()の場合でしたが、chisq.test()の場合はどのようになるか行ってください。その際の実行文と有意差があった行数を示してください（抜き出すのではなく、行数を答えてください）。
「ゲノム比較解析」の「２．DNAの塩基組成（GC含量）」に進んでください。
Page updated
Google Sites
Report abuse