Rでお手軽にアレイ解析

1. はじめに
2. Rの基本的な使い方
3. Affymetrix GeneChipデータの正規化
4. クラスタリング
5. 発現変動遺伝子の同定（ランキング）

1. はじめに

R-2.6.0-win32.exeをクリックして「実行する」を選択
聞かれるがままに「次へ」などを押してとにかくインストールを完了させる
インストールが無事完了したら、デスクトップに出現する「R2.6.0」アイコンをダブルクリックして起動
以下を、「R コンソール画面上」でコピー&ペーストする（どこからダウンロードするか？と聞かれるので、一番近いサイトを指定）
#------    ここから    ------
source("http://www.iu.a.u-tokyo.ac.jp/~kadota/R/R_functions.R")#おまじない
source("http://www.bioconductor.org/getBioC.R")                #おまじない
getBioC()                                                      #おまじない
install.packages("som")                                        #somパッケージのインストール
install.packages("qvalue")                                     #qvalueパッケージのインストール
install.packages("permtest")                                   #permtestパッケージのインストール
install.packages("pvclust")                                    #pvclustパッケージのインストール
install.packages("cclust")                                     #cclustパッケージのインストール
install.packages("e1071")                                      #e1071パッケージのインストール
install.packages("GeneTS")                                     #GeneTSパッケージのインストール
install.packages(c("corpcor", "longitudinal", "locfdr"))       #corpcor, longitudinal, locfdrパッケージのインストール
install.packages("st", dependencies=TRUE)                      #stパッケージのインストール
source("http://bioconductor.org/biocLite.R")                   #おまじない
biocLite("RankProd")                                           #RankProdパッケージのインストール

#------    ここまで    ------

2. Rの基本的な使い方

#------    ここから    ------
1+2           #1+2を計算
hoge <- 4     #hogeに4を代入
hoge          #hogeの中身を表示
sqrt(hoge)    #4のルートを計算
sqrt(4)       #4のルートを計算(当然同じ意味)
#の後の文章はコメントなので何を書いてもいいですよ。

#------    ここまで    ------

以下をコピー&ペーストしてください
#------    ここから    ------
x <- c(1, 3, 7, 9, 12, 30)
mean(x)     #平均値を計算
median(x)   #中央値を計算

#------    ここまで    ------

以下をコピー&ペースト
#------    ここから    ------
library(affy)
tukey.biweight(x)   #Tukey biweight値を計算

#------    ここまで    ------

3. Affymetrix GeneChipデータの正規化

1. 「ファイル」メニューから「ディレクトリの変更」を選び、「ブラウズ」ボタンを押して「GSE7197」ディレクトリを選択して（「C:\Documents and Settings\kadota\デスクトップ\GSE7197」などとなった状態で）OKボタンを押す
2. R code for DFW (PDF)のページ中にある二ページ分のcodeをコピペ（ペースト先はもちろんR console画面）。
   
3. 「リターン」キーを押す
4. 以下をコピペ

   #------    ここから    ------
   library(affy)       #パッケージの読み込み
   data = ReadAffy()   #CELファイルの読み込み
   data.out = expresso(data, bgcorrect.method="none", normalize.method="quantiles", pmcorrect.method="pmonly", summary.method="dfw")   #DFW法を実行し、結果をdata.outに保存。
   write.table(data.frame(exprs(data.out) ,check.names=FALSE), "data_dfw.txt",sep="\t",col.names=NA,quote=FALSE)                       #結果をdata_dfw.txtに保存
   
   #------    ここまで    ------
   

サンプルID    : 実際のサンプル名                      --> 変更後のサンプル名
----------------------------------------------------------------------------
GSM173080.CEL : Rice MH63 Control2                    --> Control2
GSM173086.CEL : Rice MH63 Control1                    --> Control1
GSM173089.CEL : Rice MH63 Control3                    --> Control3
GSM173091.CEL : Rice OsSRT1 RNAi transgenic line LM1  --> OsSRT1_LM1
GSM173093.CEL : Rice OsSRT1 RNAi transgenic line LM2  --> OsSRT1_LM2
GSM173094.CEL : Rice OsSRT1 RNAi transgenic line LM3  --> OsSRT1_LM3

1. 「ファイル」メニューから「ディレクトリの変更」を選び、「ブラウズ」ボタンを押して「GSE7197」ディレクトリを選択して（「C:\Documents and Settings\kadota\デスクトップ\GSE7197」などとなった状態で）OKボタンを押す
2. 以下をコピペ

   #------    ここから    ------
   library(affy)                               #パッケージの読み込み
   data = ReadAffy()                           #CELファイルの読み込み
   eset_rma <- rma(data)                       #RMA法を実行し、結果をeset_rmaに保存。
   write.exprs(eset_rma, file="data_rma.txt")  #結果をdata_rma.txtに保存
   
   #------    ここまで    ------
   

4. クラスタリング

#------    ここから    ------
library(pvclust)                                                                       #パッケージの読み込み
data <- read.table("data_dfw2.txt", header=TRUE, row.names=1, sep="\t")                #解析したいファイルの読み込み
data.pv <- pvclust(data, method.hclust="average", method.dist="correlation", nboot=10) #クラスタリングの実行
plot(data.pv)                                                                          #樹形図（デンドログラム）の表示

#------    ここまで    ------

#------    ここから    ------
library(pvclust)                                                                      #パッケージの読み込み
data <- read.table("data_dfw2.txt", header=TRUE, row.names=1, sep="\t")               #解析したいファイルの読み込み
data.pv <- pvclust(data, method.hclust="single", method.dist="correlation", nboot=10) #クラスタリングの実行
plot(data.pv)                                                                         #樹形図（デンドログラム）の表示

#------    ここまで    ------

#------    ここから    ------
library(pvclust)                                                                        #パッケージの読み込み
data <- read.table("data_dfw2.txt", header=TRUE, row.names=1, sep="\t")                 #解析したいファイルの読み込み
data.pv <- pvclust(data, method.hclust="complete", method.dist="correlation", nboot=10) #クラスタリングの実行
plot(data.pv)                                                                           #樹形図（デンドログラム）の表示

#------    ここまで    ------

#------    ここから    ------
library(pvclust)                                                                      #パッケージの読み込み
data <- read.table("data_dfw2.txt", header=TRUE, row.names=1, sep="\t")               #解析したいファイルの読み込み
data.pv <- pvclust(data, method.hclust="average", method.dist="correlation", nboot=30)#クラスタリングの実行
plot(data.pv)                                                                         #樹形図（デンドログラム）の表示
pvrect(data.pv, alpha=0.95, pv="bp")                                                  #ブートストラップ確率が95%より高い頑健なクラスターを四角で囲む
pvpick(data.pv, alpha=0.95, pv="bp")                                                  #得られる頑健なクラスター（ブートストラップ確率 > 95%）を構成するサンプル名を表示

#------    ここまで    ------

#------    ここから    ------
library(pvclust)                                                       #パッケージの読み込み
data <- read.table("data_dfw2.txt", header=TRUE, row.names=1, sep="\t")#解析したいファイルの読み込み
data.sub <- data[1:100,]                                               #最初の100遺伝子を抽出し、data.subに格納
data.pv <- pvclust(t(data.sub), method.hclust="average", method.dist="correlation", nboot=10)#クラスタリングの実行
plot(data.pv)                                                          #樹形図（デンドログラム）の表示
pvrect(data.pv, alpha=0.95, pv="bp")                                   #ブートストラップ確率が95%より高い頑健なクラスターを四角で囲む
pvpick(data.pv, alpha=0.95, pv="bp")                                   #得られる頑健なクラスター（ブートストラップ確率 > 95%）を構成するサンプル名を表示

#------    ここまで    ------

#------    ここから    ------
data <- read.table("data_dfw2.txt", header=TRUE, row.names=1, sep="\t") #解析したいファイルの読み込み
data.sub <- data[1:100,]                                                #最初の100遺伝子を抽出し、data.subに格納
data.dist <- dist(data.sub, method = "euclidean")                       #各サンプル間の距離を計算。デフォルトはユークリッド距離
data.hclust <- hclust(data.dist, method = "average")                    #average-linkage clustering（群平均法）を実行
plot(data.hclust, main="hclust")                                        #結果を表示

#------    ここまで    ------

#------    ここから    ------
library(pvclust)                                                                      #パッケージの読み込み
data <- read.table("data_agilent.txt", header=TRUE, row.names=1, sep="\t")            #解析したいファイルの読み込み
data.pv <- pvclust(data, method.hclust="average", method.dist="correlation", nboot=10)#クラスタリングの実行
plot(data.pv)                                                                         #樹形図（デンドログラム）の表示

#------    ここまで    ------

#------    ここから    ------
data <- read.table("data_agilent.txt", header=TRUE, row.names=1, sep="\t")#解析したいファイルの読み込み
data.dist <- dist(t(data), method = "euclidean")                          #各サンプル間の距離を計算。デフォルトはユークリッド距離
data.hclust <- hclust(data.dist, method = "average")                      #average-linkage clustering（群平均法）を実行
plot(data.hclust, main="hclust")                                          #結果を表示

#------    ここまで    ------

#------    ここから    ------
data <- read.table("data_agilent.txt", header=TRUE, row.names=1, sep="\t")#解析したいファイルの読み込み
data.sub <- data[1:90,]                                                   #最初の90遺伝子を抽出し、data.subに格納
data.dist <- dist(data.sub, method = "euclidean")                         #各サンプル間の距離を計算。デフォルトはユークリッド距離
data.hclust <- hclust(data.dist, method = "average")                      #average-linkage clustering（群平均法）を実行
plot(data.hclust, main="hclust")                                          #結果を表示

#------    ここまで    ------

#------    ここから    ------
data <- read.table("data_agilent.txt", header=TRUE, row.names=1, sep="\t")#解析したいファイルの読み込み
data.dist <- dist(data, method = "euclidean")                             #各サンプル間の距離を計算。デフォルトはユークリッド距離
data.hclust <- hclust(data.dist, method = "average")                      #average-linkage clustering（群平均法）を実行
plot(data.hclust, main="hclust")                                          #結果を表示

#------    ここまで    ------

5. 発現変動遺伝子の同定（ランキング）

5-1. Rank productで解析（Fold Changeに基づく方法）

#------    ここから    ------
library(RankProd)                                                                                 #パッケージの読み込み
data <- read.table("data_dfw2.txt", header=TRUE, row.names=1, sep="\t")                           #data_dfw2.txtファイルの読み込み
data.cl <- c(rep(0, 3), rep(1, 3))                                                                #最初の3列を"0(class1に相当)"、後の3列を1(class2に相当)"としたベクトルdata.clを作成
RP.out <- RP(data, data.cl, num.perm = 10, logged = TRUE, na.rm = FALSE, plot = FALSE, rand = 123)#Rank Product関数（RP）の実行（実際の解析ではnum.perm = 10の数値を100とかに設定すべし）
min_RPs <- apply(RP.out$RPs, 1, min)                                                              #統計量の小さいほうを抽出
tmp <- cbind(rownames(data), data, RP.out$pfp, RP.out$RPs, RP.out$RPrank, min_RPs, rank(min_RPs, ties.method = "min"))#出力用にオリジナル遺伝子発現行列（data）の右側に「FDR」,「統計量」,「統計量の順位」,「統計量の低い値」, 「統計量の低い値の順位」の情報を付加したtmpを用意。
write.table(tmp, "data_dfw2_rankprod.txt", sep = "\t", append=F, quote=F, row.names=F)            #結果をファイル（ファイル名：data_dfw2_rankprod.txt）に保存。

#------    ここまで    ------

Rank productを実行して得られたdata_dfw2_rankprod.txtをExcelで開いてみましょう。
H列：J列(or L列)で昇順にソートし、上位x位までを発現変動遺伝子とみなしたときのFalse Discovery Rate(FDR)
I列：K列(or M列)で昇順にソートし、上位x位までを発現変動遺伝子とみなしたときのFalse Discovery Rate(FDR)
J列：class2(OsSRT1_LM)の発現レベルが高い順にみたときの統計量
K列：class1(Control)の発現レベルが高い順にみたときの統計量
L列：J列の統計量の順位
M列：K列の統計量の順位

例えば「J列で昇順にソートし、上位264位までを発現変動遺伝子とみなしたときのFDRは1%(0.01)未満である」という解釈になります。
また、N, O列は後で行うt検定系の解析結果との比較を容易にするために用意したものです。
N列：J列とK列の二つの統計量の最小値（低い方の値）
O列：N列の値（統計量の低い方の値）の順位

#------    ここから    ------
library(RankProd)                                                                                 #パッケージの読み込み
data <- read.table("data_agilent.txt", header=TRUE, row.names=1, sep="\t")                        #data_agilent.txtファイルの読み込み
data.cl <- c(rep(0, 6), rep(1, 6))                                                                #最初の6列を"0(class1に相当)"、後の6列を1(class2に相当)"としたベクトルdata.clを作成
RP.out <- RP(data, data.cl, num.perm = 10, logged = TRUE, na.rm = FALSE, plot = FALSE, rand = 123)#Rank Product関数（RP）の実行（実際の解析ではnum.perm = 10の数値を100とかに設定すべし）
min_RPs <- apply(RP.out$RPs, 1, min)                                                              #統計量の小さいほうを抽出
tmp <- cbind(rownames(data), data, RP.out$pfp, RP.out$RPs, RP.out$RPrank, min_RPs, rank(min_RPs, ties.method = "min"))#出力用にオリジナル遺伝子発現行列（data）の右側に「FDR」,「統計量」,「統計量の順位」,「統計量の低い値」, 「統計量の低い値の順位」の情報を付加したtmpを用意。
write.table(tmp, "data_agilent_rankprod.txt", sep = "\t", append=F, quote=F, row.names=F)         #結果をファイル（ファイル名：data_agilent_rankprod.txt）に保存。

#------    ここまで    ------

Rank productを実行して得られたdata_agilent_rankprod.txtをExcelで開いてみましょう。
N列：P列(or R列)で昇順にソートし、上位x位までを発現変動遺伝子とみなしたときのFalse Discovery Rate(FDR)
O列：Q列(or S列)で昇順にソートし、上位x位までを発現変動遺伝子とみなしたときのFalse Discovery Rate(FDR)
P列：class2(B)の発現レベルが高い順にみたときの統計量
Q列：class1(A)の発現レベルが高い順にみたときの統計量
R列：P列の統計量の順位
S列：Q列の統計量の順位

T, U列は、後で行うt検定系の解析結果との比較を容易にするために用意したものです。
T列：P列とQ列の二つの統計量の最小値（低い方の値）
U列：T列の値（統計量の低い方の値）の順位

5-2. t-statistic系の方法で解析

#------    ここから    ------
library(st)                                  #shrinkage t statisticを計算するためのパッケージの読み込み
library(samr)                                #SAMを計算するためのパッケージの読み込み
library(limma)                               #a moderated t statistic (empirical Bayes)を計算するためのパッケージの読み込み
data <- read.table("data_dfw2.txt", header=TRUE, row.names=1, sep="\t")             #data_dfw2.txtファイルの読み込み
data.cl <- c(rep(1, 3), rep(2, 3))           #最初の3列を"1"、後の3列を"2"としたベクトルdata.clを作成
score_st <- shrinkt.stat(t(data), data.cl)   #Shrinkage t statisticを計算し結果をscore_stに格納
score_ebayes <- modt.stat(t(data), data.cl)  #a moderated t statistic (empirical Bayes)を計算し結果をscore_ebayesに格納
score_sam <- sam.stat(t(data), data.cl)      #SAM's t-statisticを計算し結果をscore_samに格納
tmp <- NULL                                  #おまじない
tmp <- cbind(rownames(data), data, as.vector(score_st), rank(-abs(score_st)), as.vector(score_ebayes), rank(-abs(score_ebayes)), as.vector(score_sam), rank(-abs(score_sam)))   #「それぞれの計算方法で得られた統計量」と「その順位」をShrinkage, moderated, SAMの順に右のカラムに結合した結果をtmpに格納。
write.table(tmp, "data_dfw2_tstat.txt", sep = "\t", append=F, quote=F, row.names=F) #tmpの中身をdata_dfw2_tstat.txtというファイル名で保存。

#------    ここまで    ------


以上を実行して得られたdata_dfw2_tstat.txtファイルをExcelで開いてみましょう。
H, I列がそれぞれshrinkage t-statisticの方法で得られた「統計量」と「その順位」を、
J, K列がそれぞれa moderated t-statisticの方法で得られた「統計量」と「その順位」を、
L, M列がそれぞれSAMの方法で得られた「統計量」と「その順位」を表しています。

#------    ここから    ------
library(st)                                  #shrinkage t statisticを計算するためのパッケージの読み込み
library(samr)                                #SAMを計算するためのパッケージの読み込み
library(limma)                               #a moderated t statistic (empirical Bayes)を計算するためのパッケージの読み込み
data <- read.table("data_agilent.txt", header=TRUE, row.names=1, sep="\t")             #data_agilent.txtファイルの読み込み
data.cl <- c(rep(1, 6), rep(2, 6))           #最初の6列を"1"、後の6列を"2"としたベクトルdata.clを作成
score_st <- shrinkt.stat(t(data), data.cl)   #Shrinkage t statisticを計算し結果をscore_stに格納
score_ebayes <- modt.stat(t(data), data.cl)  #a moderated t statistic (empirical Bayes)を計算し結果をscore_ebayesに格納
score_sam <- sam.stat(t(data), data.cl)      #SAM's t-statisticを計算し結果をscore_samに格納
tmp <- NULL                                  #おまじない
tmp <- cbind(rownames(data), data, as.vector(score_st), rank(-abs(score_st)), as.vector(score_ebayes), rank(-abs(score_ebayes)), as.vector(score_sam), rank(-abs(score_sam)))   #「それぞれの計算方法で得られた統計量」と「その順位」をShrinkage, moderated, SAMの順に右のカラムに結合した結果をtmpに格納。
write.table(tmp, "data_agilent_tstat.txt", sep = "\t", append=F, quote=F, row.names=F) #tmpの中身をdata_agilent_tstat.txtというファイル名で保存。

#------    ここまで    ------


以上を実行して得られたdata_agilent_tstat.txtファイルをExcelで開いてみましょう。
N, O列がそれぞれshrinkage t-statisticの方法で得られた「統計量」と「その順位」を、
P, Q列がそれぞれa moderated t-statisticの方法で得られた「統計量」と「その順位」を、
R, S列がそれぞれSAMの方法で得られた「統計量」と「その順位」を表しています。