====== Rによるデータ解析 ====== ===== はじめに ===== * ECCシステムでは「アプリケーション」→「R.app」にインストールされている * 統合TVによる解説動画 [[http://togotv.dbcls.jp/20090618.html|統計解析ソフト「R」の使い方 〜導入編〜]] * 教科書的なもの [[http://cse.naro.affrc.go.jp/takezawa/r-tips/r.html|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") ### ここまで ###