Rでお手軽にアレイ解析
- はじめに
- Rの基本的な使い方
- Affymetrix GeneChipデータの正規化
- クラスタリング
- 発現変動遺伝子の同定（ランキング）
1. はじめに
統計解析フリーソフトRの最新版（2007/10/16現在、R-2.6.0 for Windows）を東京大学のFTPサイトからダウンロードしたい場合：
- 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の基本的な使い方
Step 1：
Rを起動して以下をコピー&ペーストしてみてください。自分のマイクロアレイデータ解析を実際にRを用いて行う上で、このページの情報をもとに自分がどのような作業をすれば目的を達成できるかの参考になるはずです。
#------ ここから ------
1+2 #1+2を計算
hoge <- 4 #hogeに4を代入
hoge #hogeの中身を表示
sqrt(hoge) #4のルートを計算
sqrt(4) #4のルートを計算(当然同じ意味)
#の後の文章はコメントなので何を書いてもいいですよ。
#------ ここまで ------
Step 2：
遺伝子x のプローブレベルのデータがx = (1, 3, 7, 9, 12, 30)のベクトルで与えられるとき、どうやってその算術平均値や中央値を得るのかやってみます。
以下をコピー&ペーストしてください
#------ ここから ------
x <- c(1, 3, 7, 9, 12, 30)
mean(x) #平均値を計算
median(x) #中央値を計算
#------ ここまで ------
次に、中央値よりもより頑健なTukey biweight値を計算するためにはどうすればいいのでしょうか？
以下をコピー&ペースト
#------ ここから ------
library(affy)
tukey.biweight(x) #Tukey biweight値を計算
#------ ここまで ------
Step 3：
実際に自分のマイクロアレイデータを用いて解析をしようと思い始めたとき、数行を一度にコピー&ペーストすると、様々なパラメータを変更することができずにつまづくことがあるかもしれません。
そのような場合、例えば「ワードパッド」や「メモ帳」で予め作成しておいた一連のコマンド群からなるrun.txtを用意しておき、
その中身をコピー&ペーストすることでrun.txtの中身を実行します。
参考URL1：R Tips(竹澤様)
3. Affymetrix GeneChipデータの正規化
AffymetrixのGeneChipテクノロジーは
Perfect Match/Mismatch (PM/MM)戦略をとっており、データ解析時に必要な「遺伝子発現行列」を得るための前処理（NormalizationやSummarization）が必要です。
ここでは、「Rice MH63 Control」3サンプルと「Rice OsSRT1 RNAi transgenic line LM」3サンプルの実験
（Huang et al., Plant Physiol., 2007；
生データ自体はここからも取得可能）を行って得られた生データ（.CELファイル）を入力として、
DFWアルゴリズム
（とRMAアルゴリズム）を適用して出力データである
「（51,279遺伝子×計6サンプルからなる）遺伝子発現行列」を得るための手順を説明します。
ここでは、生データ（.CEL）ファイルを含むディレクトリ（ディレクトリ名：GSE7197）がデスクトップにあると仮定しています。
DFWの場合（こちらがお勧め）：
- 「ファイル」メニューから「ディレクトリの変更」を選び、「ブラウズ」ボタンを押して「GSE7197」ディレクトリを選択して（「C:\Documents and Settings\kadota\デスクトップ\GSE7197」などとなった状態で）OKボタンを押す
- R code for DFW (PDF)のページ中にある二ページ分のcodeをコピペ（ペースト先はもちろんR console画面）。
- 「リターン」キーを押す
- 以下をコピペ
#------ ここから ------
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に保存
#------ ここまで ------
（数分以上かかりますが）これでデスクトップ上の「GSE7197」ディレクトリ中に、
出力ファイル（ファイル名：data_dfw.txt）ができています。
Excelで開いて確認してみてください。このときサンプルIDと実際のサンプル名の対応関係は以下のようになっていますので、サンプル名が分かりやすいように
data_dfw2.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
RMAの場合：
- 「ファイル」メニューから「ディレクトリの変更」を選び、「ブラウズ」ボタンを押して「GSE7197」ディレクトリを選択して（「C:\Documents and Settings\kadota\デスクトップ\GSE7197」などとなった状態で）OKボタンを押す
- 以下をコピペ
#------ ここから ------
library(affy) #パッケージの読み込み
data = ReadAffy() #CELファイルの読み込み
eset_rma <- rma(data) #RMA法を実行し、結果をeset_rmaに保存。
write.exprs(eset_rma, file="data_rma.txt") #結果をdata_rma.txtに保存
#------ ここまで ------
（数分以上かかりますが）これでデスクトップ上の「GSE7197」ディレクトリ中に、
出力ファイル（ファイル名：data_rma.txt）ができています。
Excelで開いて確認してみてください。
4. クラスタリング
4-1. サンプルのクラスタリング
data_dfw2.txtのサンプルクラスタリング（どのサンプルとどのサンプルが似ているか？）を行います。
a) 類似度：1-相関係数（method.dist="correlation"）、方法：群平均法（method.hclust="average"）の場合
以下をコピペ
#------ ここから ------
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) #樹形図（デンドログラム）の表示
#------ ここまで ------
b) 類似度：1-相関係数（method.dist="correlation"）、方法：最短距離法（method.hclust="single"）の場合
以下をコピペ
#------ ここから ------
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) #樹形図（デンドログラム）の表示
#------ ここまで ------
c) 類似度：1-相関係数（method.dist="correlation"）、方法：最長距離法（method.hclust="complete"）の場合
以下をコピペ
#------ ここから ------
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) #樹形図（デンドログラム）の表示
#------ ここまで ------
図は「ファイル」メニューの「別名で保存」を選ぶことで任意の形式で保存することができます。
nbootで指定している数字（10）はブートストラップのリサンプリング回数（デフォルトは1000）です。
実際の解析ではここの数字を（マシンパワーにも依存しますがせめて100くらいは...）増やして計算を行い、クラスターのブートストラップ確率を得ることで、頑健なクラスターを知ることができます。
a)のやり方でリサンプリング回数を増やして（nboot=30）再計算を行い、頑健なクラスターを調べてみましょう。
以下をコピペ
#------ ここから ------
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%）を構成するサンプル名を表示
#------ ここまで ------
4-2. 遺伝子のクラスタリング
（私見ですが）Rで遺伝子のクラスタリングはやってもわけがわかりません。また、たった「6サンプルのクラスタリング」であれだけの時間がかかることを考えると「数万遺伝子のクラスタリング」が多大な時間を消費してしまうことが容易に想像つくと思います。
よって、ブートストラップ確率などを得るのは非現実的です。
ここでは51,279遺伝子×6サンプルからなるdata_dfw2.txtデータの最初の100遺伝子を抽出して、遺伝子のクラスタリングを行うやり方を（一応）示します。
a) 類似度：1-相関係数（method.dist="correlation"）、方法：群平均法（method.hclust="average"）の場合（pvclust関数）
以下をコピペ
#------ ここから ------
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%）を構成するサンプル名を表示
#------ ここまで ------
b) 類似度：ユークリッド距離（method="euclidean"）、方法：群平均法（method="average"）の場合（hclust関数）
以下をコピペ
#------ ここから ------
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") #結果を表示
#------ ここまで ------
Tips：dist関数で、他にどんな類似度を指定できるかを知りたいときは、"help(dist)"と打てば分かります。
4-3. Agilentデータのクラスタリング
Agilentデータのクラスタリングも、Affymetrixのときとなんら変わりません。例えば43734遺伝子×12サンプルのlog ratioデータ（ファイル名：data_agilent.txt）がデスクトップ上の「GSE7197」ディレクトリ上にあるとします。
a) pvclust関数を用いたサンプルのクラスタリングの場合
以下をコピペ
#------ ここから ------
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) #樹形図（デンドログラム）の表示
#------ ここまで ------
b) hclust関数を用いたサンプルのクラスタリングの場合
以下をコピペ
#------ ここから ------
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") #結果を表示
#------ ここまで ------
c) hclust関数を用いた（最初の90）遺伝子のクラスタリングの場合
以下をコピペ
#------ ここから ------
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") #結果を表示
#------ ここまで ------
d) 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に基づく方法）
a) Affymetrix GeneChipデータの場合（51,279遺伝子×6サンプルからなるdata_dfw2.txtデータ）
目的：最初の3列とその後の3列のサンプル間で発現が異なる遺伝子をランキングしたい
以下をコピペ
#------ ここから ------
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列の値（統計量の低い方の値）の順位
b) Agilentデータの場合（43734遺伝子×12サンプルのlog ratioデータ; data_agilent.txt）
目的：最初の6列とその後の6列のサンプル間で発現が異なる遺伝子をランキングしたい
(log(A1/ref), log(A2/ref), ..., log(A6/ref), log(B1/ref), log(B2/ref), ..., log(B6/ref)のようなサンプルの順番だとする)
以下をコピペ
#------ ここから ------
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系の方法で解析
a) Affymetrix GeneChipデータの場合（51,279遺伝子×6サンプルからなるdata_dfw2.txtデータ）
目的：最初の3列とその後の3列のサンプル間で発現が異なる遺伝子をshrinkage t-statistic、a moderated t-statistic、SAMでランキングしたい
以下をコピペ
#------ ここから ------
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の方法で得られた「統計量」と「その順位」を表しています。
b) Agilentデータの場合（43734遺伝子×12サンプルのlog ratioデータ; data_agilent.txt）
目的：最初の6列とその後の6列のサンプル間で発現が異なる遺伝子をshrinkage t-statistic、a moderated t-statistic、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の方法で得られた「統計量」と「その順位」を表しています。
参考URL2：Rでマイクロアレイデータ解析
Last updated on 2007.10.19