• ECCシステムでは「アプリケーション」→「R.app」にインストールされている
• 統合TVによる解説動画 統計解析ソフト「R」の使い方 〜導入編〜
• 教科書的なもの R-Tips

# アレイ解析のパッケージlimmaのインストール
# 初回のみ必要で、2回目からは実行しなくてもよい
source("http://bioconductor.org/biocLite.R")
biocLite("limma")

# アレイ解析のパッケージlimmaを読み込む
library(limma)

# 数値化データが置いてある場所を指定
datapath <- "/home08/sirna/Microarray/2015-06-12"

# 数値化データのファイル名を記載したsample_list.txtを読み込む。こんなの：
#
# ===== sample_list.txt =====
# SampleNumber	FileName	Condition
# 1	US91503671_253949442637_S01_GE1_105_Dec08_1_1.txt	mock_1
# 2	US91503671_253949442637_S01_GE1_105_Dec08_1_2.txt	RNA
# 3	US91503671_253949442637_S01_GE1_105_Dec08_1_3.txt	2OMe3
# 4	US91503671_253949442637_S01_GE1_105_Dec08_1_4.txt	2OMe5
# 5	US91503671_253949442637_S01_GE1_105_Dec08_2_1.txt	mock_2
# 6	US91503671_253949442637_S01_GE1_105_Dec08_2_2.txt	2OMe7
# 7	US91503671_253949442637_S01_GE1_105_Dec08_2_3.txt	LNA3
# 8	US91503671_253949442637_S01_GE1_105_Dec08_2_4.txt	LNA7
# ===== sample_list.txt =====
#
targets <- readTargets("sample_list.txt", path=datapath)

# フィルターを定義。GeneSpringで設定した内容と同等
myFlagFun <- function(x){
	control   <- x$ControlType      == 0
	present   <- x$gIsPosAndSignif  == 1
	uniform   <- x$gIsFeatNonUnifOL == 0
	wellabove <- x$gIsWellAboveBG   == 1
	saturated <- x$gIsSaturated     == 0
	popnol    <- x$gIsFeatPopnOL    == 0
	y <- as.numeric(
		control & present & uniform & wellabove & saturated & popnol
	)
	y
}

# 生データ読み込み。GとRが必須だが1色法なのでRはダミーデータを与えておく
RG <- read.maimages(targets, path=datapath,
	columns = list(G="gProcessedSignal",R="gProcessedSignal"),
	annotation = c("Row","Col","FeatureNum","ControlType","ProbeName","SystematicName"),
	wt.fun=myFlagFun
)

# 操作しやすいようにカラム名を変更
colnames(RG$G)       <- RG$targets$Condition
colnames(RG$R)       <- RG$targets$Condition
colnames(RG$weights) <- RG$targets$Condition

# シグナル強度の分布をプロット
plotDensities(RG)

# サンプル間でquantile正規化
RG$G <- normalizeBetweenArrays(RG$G, method="quantile")
RG$R <- normalizeBetweenArrays(RG$R, method="quantile")

# シグナル強度の分布をプロット
plotDensities(RG)

# log2を取る
RG$G <- log2(RG$G)
RG$R <- log2(RG$R)

# 同一プローブが複数ある場合に平均を取る。
# RG.avg <- avereps(RG, ID=RG$genes$ProbeName)

# プロットを作成
sample1name <- "mock_1"
sample2name <- "RNA"

xa <- NULL
ya <- NULL
xp <- NULL
yp <- NULL

# フィルターを通った「使える」値
pos_P <- grep(1,RG$weight[,c(sample1name)] * RG$weight[,c(sample2name)])

# フィルターではじかれた値
pos_A <- grep(0,RG$weight[,c(sample1name)] * RG$weight[,c(sample2name)])

xa <- RG$G[pos_A,c(sample1name)]
ya <- RG$G[pos_A,c(sample2name)]

xp <- RG$G[pos_P,c(sample1name)]
yp <- RG$G[pos_P,c(sample2name)]

### XYプロットを描画する。###
plot(xa,ya,
	xlim=c(0,20),
	ylim=c(0,20),
	xlab="log2(sample1)",
	ylab="log2(sample2)",
	title(paste(sample1name, " vs ", sample2name)),
	col="lightblue",pch=20,cex=0.5
)
par(new=T)
plot(xp,yp,
	xlim=c(0,20),ylim=c(0,20),
	xlab="",ylab="",
	col="gray50",pch=20,cex=0.5
)
lines(c(-100,100),c(-100,100),col="brown")
lines(c(-100,100),c( -99,101),col="brown")
lines(c(-100,100),c(-101, 99),col="brown")
### ここまで ###

### MAプロットを描画する。###
plot((xa+ya)/2,ya-xa,
	xlim=c(0,20),ylim=c(-7,7),
	xlab="A: (log2(sample1) + log2(sample2))/2",
	ylab="M: log2(sample2) - log2(sample1)",
	title(paste(sample1name, " vs ", sample2name)),
	col="lightblue",pch=20,cex=0.5
)
par(new=T)
plot((xp+yp)/2,yp-xp,
	xlim=c(0,20),ylim=c(-7,7),
	xlab="",ylab="",
	col="gray50",pch=20,cex=0.5
)
lines(c(-100,100),c( 0, 0),col="brown")
lines(c(-100,100),c(-1,-1),col="brown")
lines(c(-100,100),c( 1, 1),col="brown")
### ここまで ###

# 6つのサンプルともフィルターを通過したもの（フラグの積が0でない）を得る
pos_detected <- grep(1,
	RG$weight[,c("mock_1")] *
	RG$weight[,c("RNA"   )] *
	RG$weight[,c("2OMe3" )] *
	RG$weight[,c("2OMe5" )] *
	RG$weight[,c("mock_2")] *
	RG$weight[,c("2OMe7" )] *
	RG$weight[,c("LNA3"  )] *
	RG$weight[,c("LNA7"  )]
)

# サンプルの指定
x_name <- "mock_1"
y_name <- "RNA"

# 遺伝子リストを読み込む
genelist <- read.delim("/home08/sirna/Microarray/2015-06-12/genelist/guide_seed_genelist.txt")

# 両方のサンプルでフィルターを通過したもののシグナル値を得る。
x <- RG$G[pos_detected,c(x_name)]
y <- RG$G[pos_detected,c(y_name)]

# 参考：mock_1, mock_2の平均をとる場合は、上記のxのかわりに
# x1 <- RG$G[pos_detected,c("mock_1")]
# x2 <- RG$G[pos_detected,c("mock_2")]
# x  <- (x1 + x2)/2
# のようにできる。

# 列の名前をつける
colnames(genelist) <- c("RefSeqID","NumHits")

# RefSeqIDは因子としてデータフレームに格納されているので、文字列に変換する。
genelist$RefSeqID <- as.vector(genelist$RefSeqID)

# アレイデータのSystematicName（RefSeqIDが入っている）と遺伝子リストのRefSeqIDを紐付ける。
arraydata <- cbind(RG$gene[pos_detected,],RG$G[pos_detected,])
merged    <- merge(arraydata, genelist, by.x="SystematicName", by.y="RefSeqID", all.x=T, sort=F)

# mock_1（1班）と RNA（2班）の値で、seedをもつものを取り出す。
x_seedmatch <- subset(merged, NumHits > 0)[,c(x_name)]
y_seedmatch <- subset(merged, NumHits > 0)[,c(y_name)]

# 参考：mock_1, mock_2の平均をとる場合は、上記のx_seedmatchのかわりに
# x1_seedmatch <- subset(merged, NumHits > 0)[,c("mock_1")]
# x2_seedmatch <- subset(merged, NumHits > 0)[,c("mock_2")]
# x_seedmatch  <- (x1_seedmatch + x2_seedmatch)/2
# のようにできる。

### MAプロットを描画する。###
plot((x+y)/2,y-x,
	xlim=c(0,20),ylim=c(-4,4),
	xlab="A: (log2(sample1) + log2(sample2))/2",
	ylab="M: log2(sample2) - log2(sample1)",
	title(paste(sample1name, " vs ", sample2name)),
	col="lightblue",pch=20,cex=0.5
)
par(new=T)
plot((x_seedmatch+y_seedmatch)/2,y_seedmatch-x_seedmatch,
	xlim=c(0,20),ylim=c(-4,4),
	xlab="",ylab="",
	col="blue",pch=20,cex=0.5
)
lines(c(-100,100),c( 0, 0),col="brown")
lines(c(-100,100),c(-1,-1),col="brown")
lines(c(-100,100),c( 1, 1),col="brown")
### ここまで ###

### 累積度数曲線を描画する ###
plot(
	ecdf( y - x ),
	verticals=TRUE,
	do.p=FALSE,
	xlim=c(-1,1), ylim=c(0,1),
	col.hor="black", col.vert="black",
	main=paste(x_name, " vs ", y_name), xlab="", ylab=""
)
par(new=T)  # 重ねあわせる。リフレッシュしない。
plot(
	ecdf( y_seedmatch - x_seedmatch ),
	verticals=TRUE,
	do.p=FALSE,
	xlim=c(-1,1), ylim=c(0,1),
	col.hor="red", col.vert="red",
	main="", xlab="", ylab="",
	axes=F
)
lines(c(0,0),c(-1,2),col="brown")
lines(c(-10,10),c(0.1,0.1),col="lightblue")
lines(c(-10,10),c(0.2,0.2),col="lightblue")
lines(c(-10,10),c(0.3,0.3),col="lightblue")
lines(c(-10,10),c(0.4,0.4),col="lightblue")
lines(c(-10,10),c(0.5,0.5),col="lightblue")
lines(c(-10,10),c(0.6,0.6),col="lightblue")
lines(c(-10,10),c(0.7,0.7),col="lightblue")
lines(c(-10,10),c(0.8,0.8),col="lightblue")
lines(c(-10,10),c(0.9,0.9),col="lightblue")
### ここまで ###