(Rで)マイクロアレイデータ解析(last modified 2023/05/24, since 2005)

in_f <- "sample14.txt"                 #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- 0.2                           #CVの閾値を指定

#必要なパッケージをロード
library(genefilter)                    #パッケージの読み込み

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
obj <- genefilter(data, cv(param, Inf))#条件を満たすかどうかを判定した結果をobjに格納
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
dim(data)                              #確認してるだけです

#ファイルに保存
tmp <- cbind(rownames(data), data)     #保存したい情報をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample14.txt"                 #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- 0.2                           #CVの閾値を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
hoge <- apply(data,1,sd)/apply(data,1,mean)#CVを計算した結果をhogeに格納
obj <- (hoge >= param)                 #条件を満たすかどうかを判定した結果をobjに格納
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
dim(data)                              #確認してるだけです

#ファイルに保存
tmp <- cbind(rownames(data), data)     #保存したい情報をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample14.txt"                 #入力ファイル名を指定してin_fに格納
out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param <- 0.1                           #フィルタリングしたい下位x%を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
hoge <- apply(data,1,sd)/apply(data,1,mean)#CVを計算した結果をhogeに格納
hoge2 <- quantile(hoge, probs=param)   #paramで指定した下位の位置情報をhogeに格納
obj <- (hoge >= hoge2)                 #条件を満たすかどうかを判定した結果をobjに格納
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
summary(hoge)                          #確認してるだけです
dim(data)                              #確認してるだけです

#ファイルに保存
tmp <- cbind(rownames(data), data)     #保存したい情報をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample14.txt"                 #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- 0.2                           #閾値を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
hoge <- apply(data, 1, var)            #分散を計算した結果をhogeに格納
obj <- (hoge >= param)                 #条件を満たすかどうかを判定した結果をobjに格納
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
summary(hoge)                          #確認してるだけです
dim(data)                              #確認してるだけです

#ファイルに保存
tmp <- cbind(rownames(data), data)     #保存したい情報をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample14.txt"                 #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- 0.1                           #フィルタリングしたい下位x%を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
hoge <- apply(data, 1, var)            #分散を計算した結果をhogeに格納
hoge2 <- quantile(hoge, probs=param)   #paramで指定した下位の位置情報をhogeに格納
obj <- (hoge >= hoge2)                 #条件を満たすかどうかを判定した結果をobjに格納
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
summary(hoge)                          #確認してるだけです
dim(data)                              #確認してるだけです

#ファイルに保存
tmp <- cbind(rownames(data), data)     #保存したい情報をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2_BAT.txt"           #入力ファイル名を指定してin_fに格納
out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param <- 0.2                           #フィルタリングしたい下位x%を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
hoge <- apply(data, 1, var)            #分散を計算した結果をhogeに格納
hoge2 <- quantile(hoge, probs=param)   #paramで指定した下位の位置情報をhogeに格納
obj <- (hoge >= hoge2)                 #条件を満たすかどうかを判定した結果をobjに格納
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
summary(hoge)                          #確認してるだけです
dim(data)                              #確認してるだけです

#ファイルに保存
tmp <- cbind(rownames(data), data)     #保存したい情報をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

param1 <- "rat2302.db"                 #パッケージ名を指定

#必要なパッケージのインストール
source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
biocLite(param1, suppressUpdates=TRUE) #指定したパッケージのインストール
	

param1 <- "rat2302.db"                 #パッケージ名を指定

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#本番(利用可能なアノテーション情報を表示)
keytypes(eval(parse(text=param1)))     #利用可能なアノテーション情報を表示
	

in_f1 <- "data_rma_2.txt"              #入力ファイル名を指定してin_f1に格納(発現データ)
in_f2 <- "hoge3_GPL1355.txt"           #入力ファイル名を指定してin_f2に格納(Gene symbolとIDの対応表のデータ)
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定

#前処理(IDとGene symbolとの対応関係を含む情報を入手)
sym <- read.table(in_f2, header=TRUE, row.names=1, sep="\t", quote="")#in_f2で指定したファイルの読み込み
IDs <- as.vector(sym[,1])              #Gene symbol情報をベクトルに変換し、IDsに格納
names(IDs) <- rownames(sym)            #IDsを行名で対応づけられるようにしている
uniqID <- unique(IDs)                  #non-redundant ID情報を抽出し、uniqIDに格納
uniqID <- uniqID[uniqID != ""]         #uniqIDの中から指定したIDがないものを除く
uniqID <- uniqID[!is.na(uniqID)]       #uniqIDの中から指定したIDが"NA"のものを除く
uniqID <- uniqID[!is.nan(uniqID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

#本番
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")#in_f1で指定したファイルの読み込み
hoge <- t(apply(as.matrix(uniqID), 1, function(i, d = data, s = IDs, p = param) {#uniqIDを一つずつ処理する
                apply(d[which(s == i), ], 2, p, na.rm = TRUE)#dataの中から現在のgene symbol(i)と同じprobesを全て抽出し、その平均値(mean)を返す
        }, data, IDs, param))          #apply関数でdata,IDs,paramを使えるように代入
rownames(hoge) <- uniqID               #non-redundant IDをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #指定したIDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f1 <- "data_rma_2.txt"              #入力ファイル名を指定してin_f1に格納(発現データ)
in_f2 <- "hoge3_GPL1355.txt"           #入力ファイル名を指定してin_f2に格納(Gene symbolとIDの対応表のデータ)
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定

#前処理(IDとGene symbolとの対応関係を含む情報を入手)
sym <- read.table(in_f2, header=TRUE, row.names=1, sep="\t", quote="")#in_f2で指定したファイルの読み込み
IDs <- as.vector(sym[,1])              #Gene symbol情報をベクトルに変換し、IDsに格納
names(IDs) <- rownames(sym)            #IDsを行名で対応づけられるようにしている
uniqID <- unique(IDs)                  #non-redundant ID情報を抽出し、uniqIDに格納
uniqID <- uniqID[uniqID != ""]         #uniqIDの中から指定したIDがないものを除く
uniqID <- uniqID[!is.na(uniqID)]       #uniqIDの中から指定したIDが"NA"のものを除く
uniqID <- uniqID[!is.nan(uniqID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

#本番
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")#in_f1で指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(uniqID)){            #non-redundant ID数分だけループを回す
  obj <- is.element(IDs, uniqID[i])    #条件を満たすかどうかを判定した結果をobjに格納(IDsベクトル内でi番目のuniqIDと一致する位置をTRUE)
  hoge <- rbind(hoge, apply(data[names(IDs[obj]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
  #hoge <- rbind(hoge, apply(data[which(IDs == uniqID[i]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
}                                      #non-redundant ID数分だけループを回す
rownames(hoge) <- uniqID               #non-redundant IDをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #指定したIDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f1 <- "sample13_7vs7.txt"           #入力ファイル名を指定してin_f1に格納(発現データ)
in_f2 <- "hoge6_GPL887.txt"            #入力ファイル名を指定してin_f2に格納(Gene symbolとIDの対応表のデータ)
out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定

#前処理(IDとGene symbolとの対応関係を含む情報を入手)
sym <- read.table(in_f2, header=TRUE, row.names=1, sep="\t", quote="")#in_f2で指定したファイルの読み込み
IDs <- as.vector(sym[,1])              #Gene symbol情報をベクトルに変換し、IDsに格納
names(IDs) <- rownames(sym)            #IDsを行名で対応づけられるようにしている
uniqID <- unique(IDs)                  #non-redundant ID情報を抽出し、uniqIDに格納
uniqID <- uniqID[uniqID != ""]         #uniqIDの中から指定したIDがないものを除く
uniqID <- uniqID[!is.na(uniqID)]       #uniqIDの中から指定したIDが"NA"のものを除く
uniqID <- uniqID[!is.nan(uniqID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

#本番
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")#in_f1で指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(uniqID)){            #non-redundant ID数分だけループを回す
  obj <- is.element(IDs, uniqID[i])    #条件を満たすかどうかを判定した結果をobjに格納(IDsベクトル内でi番目のuniqIDと一致する位置をTRUE)
  hoge <- rbind(hoge, apply(data[names(IDs[obj]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
  #hoge <- rbind(hoge, apply(data[which(IDs == uniqID[i]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
}                                      #non-redundant ID数分だけループを回す
rownames(hoge) <- uniqID               #non-redundant IDをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #指定したIDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2.txt"               #入力ファイル名を指定してin_fに格納
out_f <- "hoge4.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定
param1 <- "rat2302.db"                 #パッケージ名を指定
param2 <- "SYMBOL"                     #欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

#必要なパッケージのインストール(必要に応じて#を外して実行すべし)
#source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
#biocLite(param1, suppressUpdates=TRUE)#指定したパッケージのインストール

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#probe IDと指定したIDとの対応関係を含む情報を入手
hoge <- sub(".db", param2, param1)     #用いたいオブジェクト名情報を取得(hgu133aSYMBOLに相当)
hoge <- eval(parse(text=hoge))         #hogeを文字列としてではなくオブジェクトとして取り扱う
IDs <- unlist(as.list(hoge))           #全probesに対応する指定したID情報を抽出し、IDsに格納
uniqID <- unique(IDs)                  #non-redundant ID情報を抽出し、uniqIDに格納
uniqID <- uniqID[uniqID != ""]         #uniqIDの中から指定したIDがないものを除く
uniqID <- uniqID[!is.na(uniqID)]       #uniqIDの中から指定したIDが"NA"のものを除く
uniqID <- uniqID[!is.nan(uniqID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

#本番
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(uniqID)){            #non-redundant ID数分だけループを回す
  obj <- is.element(IDs, uniqID[i])    #条件を満たすかどうかを判定した結果をobjに格納(IDsベクトル内でi番目のuniqIDと一致する位置をTRUE)
  hoge <- rbind(hoge, apply(data[names(IDs[obj]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
  #hoge <- rbind(hoge, apply(data[which(IDs == uniqID[i]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
}                                      #non-redundant ID数分だけループを回す
rownames(hoge) <- uniqID               #non-redundant IDをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #指定したIDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample5.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge5.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定
param1 <- "hgu133a.db"                 #パッケージ名を指定
param2 <- "SYMBOL"                     #欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

#必要なパッケージのインストール(必要に応じて#を外して実行すべし)
#source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
#biocLite(param1, suppressUpdates=TRUE)#指定したパッケージのインストール

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#probe IDと指定したIDとの対応関係を含む情報を入手
hoge <- sub(".db", param2, param1)     #用いたいオブジェクト名情報を取得(hgu133aSYMBOLに相当)
hoge <- eval(parse(text=hoge))         #hogeを文字列としてではなくオブジェクトとして取り扱う
IDs <- unlist(as.list(hoge))           #全probesに対応する指定したID情報を抽出し、IDsに格納
uniqID <- unique(IDs)                  #non-redundant ID情報を抽出し、uniqIDに格納
uniqID <- uniqID[uniqID != ""]         #uniqIDの中から指定したIDがないものを除く
uniqID <- uniqID[!is.na(uniqID)]       #uniqIDの中から指定したIDが"NA"のものを除く
uniqID <- uniqID[!is.nan(uniqID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

#本番
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(uniqID)){            #non-redundant ID数分だけループを回す
  obj <- is.element(IDs, uniqID[i])    #条件を満たすかどうかを判定した結果をobjに格納(IDsベクトル内でi番目のuniqIDと一致する位置をTRUE)
  hoge <- rbind(hoge, apply(data[names(IDs[obj]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
  #hoge <- rbind(hoge, apply(data[which(IDs == uniqID[i]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
}                                      #non-redundant ID数分だけループを回す
rownames(hoge) <- uniqID               #non-redundant IDをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #指定したIDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2.txt"               #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定
param1 <- "rat2302.db"                 #パッケージ名を指定
param2 <- "ENTREZID"                   #欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

#必要なパッケージのインストール(必要に応じて#を外して実行すべし)
#source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
#biocLite(param1, suppressUpdates=TRUE)#指定したパッケージのインストール

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#probe IDと指定したIDとの対応関係を含む情報を入手
hoge <- sub(".db", param2, param1)     #用いたいオブジェクト名情報を取得(rat2302ENTREZIDに相当)
hoge <- eval(parse(text=hoge))         #hogeを文字列としてではなくオブジェクトとして取り扱う
IDs <- unlist(as.list(hoge))           #全probesに対応する指定したID情報を抽出し、IDsに格納
uniqID <- unique(IDs)                  #non-redundant ID情報を抽出し、uniqIDに格納
uniqID <- uniqID[uniqID != ""]         #uniqIDの中から指定したIDがないものを除く
uniqID <- uniqID[!is.na(uniqID)]       #uniqIDの中から指定したIDが"NA"のものを除く
uniqID <- uniqID[!is.nan(uniqID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

#本番
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(uniqID)){            #non-redundant ID数分だけループを回す
    hoge <- rbind(hoge, apply(data[which(IDs == uniqID[i]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
}                                      #non-redundant ID数分だけループを回す
rownames(hoge) <- uniqID               #non-redundant IDをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #指定したIDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample5.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定
param1 <- "hgu133a.db"                 #パッケージ名を指定
param2 <- "ENTREZID"                   #欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

#必要なパッケージのインストール(必要に応じて#を外して実行すべし)
#source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
#biocLite(param1, suppressUpdates=TRUE)#指定したパッケージのインストール

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#probe IDと指定したIDとの対応関係を含む情報を入手
hoge <- sub(".db", param2, param1)     #用いたいオブジェクト名情報を取得hgu133aENTREZIDに相当)
hoge <- eval(parse(text=hoge))         #hogeを文字列としてではなくオブジェクトとして取り扱う
IDs <- unlist(as.list(hoge))           #全probesに対応する指定したID情報を抽出し、IDsに格納
uniqID <- unique(IDs)                  #non-redundant ID情報を抽出し、uniqIDに格納
uniqID <- uniqID[uniqID != ""]         #uniqIDの中から指定したIDがないものを除く
uniqID <- uniqID[!is.na(uniqID)]       #uniqIDの中から指定したIDが"NA"のものを除く
uniqID <- uniqID[!is.nan(uniqID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

#本番
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(uniqID)){            #non-redundant ID数分だけループを回す
    hoge <- rbind(hoge, apply(data[which(IDs == uniqID[i]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
}                                      #non-redundant ID数分だけループを回す
rownames(hoge) <- uniqID               #non-redundant IDをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #指定したIDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2.txt"               #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定
param1 <- "rat2302.db"                 #パッケージ名を指定
param2 <- "ENSEMBL"                    #欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

#必要なパッケージのインストール(必要に応じて#を外して実行すべし)
#source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
#biocLite(param1, suppressUpdates=TRUE)#指定したパッケージのインストール

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#probe IDと指定したIDとの対応関係を含む情報を入手
hoge <- sub(".db", param2, param1)     #用いたいオブジェクト名情報を取得(rat2302ENSEMBLに相当)
hoge <- eval(parse(text=hoge))         #hogeを文字列としてではなくオブジェクトとして取り扱う
IDs <- unlist(as.list(hoge))           #全probesに対応する指定したID情報を抽出し、IDsに格納
uniqID <- unique(IDs)                  #non-redundant ID情報を抽出し、uniqIDに格納
uniqID <- uniqID[uniqID != ""]         #uniqIDの中から指定したIDがないものを除く
uniqID <- uniqID[!is.na(uniqID)]       #uniqIDの中から指定したIDが"NA"のものを除く
uniqID <- uniqID[!is.nan(uniqID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

#本番
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(uniqID)){            #non-redundant ID数分だけループを回す
    hoge <- rbind(hoge, apply(data[which(IDs == uniqID[i]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
}                                      #non-redundant ID数分だけループを回す
rownames(hoge) <- uniqID               #non-redundant IDをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #指定したIDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample5.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定
param1 <- "hgu133a.db"                 #パッケージ名を指定
param2 <- "UNIGENE"                    #欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

#必要なパッケージのインストール(必要に応じて#を外して実行すべし)
#source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
#biocLite(param1, suppressUpdates=TRUE)#指定したパッケージのインストール

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#probe IDと指定したIDとの対応関係を含む情報を入手
hoge <- sub(".db", param2, param1)     #用いたいオブジェクト名情報を取得(hgu133aUNIGENEに相当)
hoge <- eval(parse(text=hoge))         #hogeを文字列としてではなくオブジェクトとして取り扱う
IDs <- unlist(as.list(hoge))           #全probesに対応する指定したID情報を抽出し、IDsに格納
uniqID <- unique(IDs)                  #non-redundant ID情報を抽出し、uniqIDに格納
uniqID <- uniqID[uniqID != ""]         #uniqIDの中から指定したIDがないものを除く
uniqID <- uniqID[!is.na(uniqID)]       #uniqIDの中から指定したIDが"NA"のものを除く
uniqID <- uniqID[!is.nan(uniqID)]      #uniqIDの中から指定したIDが"NaN"のものを除く

#本番
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(uniqID)){            #non-redundant ID数分だけループを回す
    hoge <- rbind(hoge, apply(data[which(IDs == uniqID[i]),], 2, param, na.rm=TRUE))#dataの中からi番目の指定したIDと同じprobesを全て抽出し、そのparamで指定した代表値をhogeの一番下の行に追加
}                                      #non-redundant ID数分だけループを回す
rownames(hoge) <- uniqID               #non-redundant IDをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #指定したIDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f1 <- "data_rma_2.txt"              #入力ファイル名を指定してin_f1に格納(発現データ)
in_f2 <- "GPL1355-14795_symbol.txt"    #入力ファイル名を指定してin_f2に格納(Gene symbolとIDの対応表のデータ)
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定

#前処理(IDとGene symbolとの対応関係を含む情報を入手)
sym <- read.table(in_f2, header=TRUE, row.names=1, sep="\t", quote="")#in_f2で指定したファイルの読み込み
symbols <- as.vector(sym[,1])          #Gene symbol情報をベクトルに変換し、symbolsに格納
names(symbols) <- rownames(sym)        #symbolsをIDで対応づけられるようにしている
unique_sym <- unique(symbols)          #symbolsの中からnon-redundantな情報のみを抽出し、unique_symに格納
unique_sym <- unique_sym[unique_sym != ""]#unique_symの中から、Gene symbolがないものを除く
unique_sym <- unique_sym[!is.na(unique_sym)]#unique_symの中から、Gene symbolが"NA"のものを除く
unique_sym <- unique_sym[!is.nan(unique_sym)]#unique_symの中から、Gene symbolが"NaN"のものを除く

#本番
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")#in_f1で指定したファイルの読み込み
hoge <- t(apply(as.matrix(unique_sym), 1, function(i, d = data, s = symbols, p = param) {#unique_symを一つずつ処理する
                apply(d[which(s == i), ], 2, p, na.rm = TRUE)#dataの中から現在のgene symbol(i)と同じprobesを全て抽出し、その平均値(mean)を返す
        }, data, symbols, param))      #apply関数でdata,symbols,paramを使えるように代入
rownames(hoge) <- unique_sym           #non-redundant gene symbolsをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #遺伝子IDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f1 <- "data_rma_2.txt"              #入力ファイル名を指定してin_f1に格納(発現データ)
in_f2 <- "GPL1355-14795_symbol.txt"    #入力ファイル名を指定してin_f2に格納(Gene symbolとIDの対応表のデータ)
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定

#前処理(IDとGene symbolとの対応関係を含む情報を入手)
sym <- read.table(in_f2, header=TRUE, row.names=1, sep="\t", quote="")#in_f2で指定したファイルの読み込み
symbols <- as.vector(sym[,1])          #Gene symbol情報をベクトルに変換し、symbolsに格納
names(symbols) <- rownames(sym)        #symbolsをIDで対応づけられるようにしている
unique_sym <- unique(symbols)          #symbolsの中からnon-redundantな情報のみを抽出し、unique_symに格納
unique_sym <- unique_sym[unique_sym != ""]#unique_symの中から、Gene symbolがないものを除く
unique_sym <- unique_sym[!is.na(unique_sym)]#unique_symの中から、Gene symbolが"NA"のものを除く
unique_sym <- unique_sym[!is.nan(unique_sym)]#unique_symの中から、Gene symbolが"NaN"のものを除く

#本番
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")#in_f1で指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(unique_sym)){        #non-redundant gene symbol数分だけループを回す
    hoge <- rbind(hoge, apply(data[which(symbols == unique_sym[i]),], 2, param, na.rm=TRUE))#dataの中からi番目のgene symbolと同じprobesを全て抽出し、その平均値(mean)をhogeの一番下の行に追加
}                                      #non-redundant gene symbol数分だけループを回す
rownames(hoge) <- unique_sym           #non-redundant gene symbolsをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #遺伝子IDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f1 <- "sample13_7vs7.txt"           #入力ファイル名を指定してin_f1に格納(発現データ)
in_f2 <- "GPL887-5640_symbol.txt"      #入力ファイル名を指定してin_f2に格納(Gene symbolとIDの対応表のデータ)
out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定

#前処理(IDとGene symbolとの対応関係を含む情報を入手)
sym <- read.table(in_f2, header=TRUE, row.names=1, sep="\t", quote="")#in_f2で指定したファイルの読み込み
symbols <- as.vector(sym[,1])          #Gene symbol情報をベクトルに変換し、symbolsに格納
names(symbols) <- rownames(sym)        #symbolsをIDで対応づけられるようにしている
unique_sym <- unique(symbols)          #symbolsの中からnon-redundantな情報のみを抽出し、unique_symに格納
unique_sym <- unique_sym[unique_sym != ""]#unique_symの中から、Gene symbolがないものを除く
unique_sym <- unique_sym[!is.na(unique_sym)]#unique_symの中から、Gene symbolが"NA"のものを除く
unique_sym <- unique_sym[!is.nan(unique_sym)]#unique_symの中から、Gene symbolが"NaN"のものを除く

#本番
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")#in_f1で指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(unique_sym)){        #non-redundant gene symbol数分だけループを回す
    hoge <- rbind(hoge, apply(data[which(symbols == unique_sym[i]),], 2, param, na.rm=TRUE))#dataの中からi番目のgene symbolと同じprobesを全て抽出し、その平均値(mean)をhogeの一番下の行に追加
}                                      #non-redundant gene symbol数分だけループを回す
rownames(hoge) <- unique_sym           #non-redundant gene symbolsをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #遺伝子IDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2.txt"               #入力ファイル名を指定してin_fに格納
out_f <- "hoge4.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定
param1 <- "rat2302.db"                 #パッケージ名を指定
param2 <- "SYMBOL"                     #欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

#必要なパッケージのインストール(利用したいパッケージが既に存在していれば2回目以降は必要なし)
source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
biocLite(param1, suppressUpdates=TRUE) #指定したパッケージのインストール

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#probe IDとGene symbolとの対応関係を含む情報を入手
hoge <- sub(".db", param2, param1)     #用いたいオブジェクト名情報を取得(rat2302SYMBOLに相当)
hoge <- eval(parse(text=hoge))         #hogeを文字列としてではなくオブジェクトとして取り扱う
symbols <- unlist(as.list(hoge))       #全31099 probesに対応するGeneSymbol(SYMBOL)情報を抽出し、symbolsに格納
unique_sym <- unique(symbols)          #non-redundantなGeneSymbol(SYMBOL)情報を抽出し、unique_symに格納
unique_sym <- unique_sym[unique_sym != ""]#unique_symの中から、Gene symbolがないものを除く
unique_sym <- unique_sym[!is.na(unique_sym)]#unique_symの中から、Gene symbolが"NA"のものを除く
unique_sym <- unique_sym[!is.nan(unique_sym)]#unique_symの中から、Gene symbolが"NaN"のものを除く

#本番
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(unique_sym)){        #non-redundant gene symbol数分だけループを回す
    hoge <- rbind(hoge, apply(data[which(symbols == unique_sym[i]),], 2, param, na.rm=TRUE))#dataの中からi番目のgene symbolと同じprobesを全て抽出し、その平均値(mean)をhogeの一番下の行に追加
}                                      #non-redundant gene symbol数分だけループを回す
rownames(hoge) <- unique_sym           #non-redundant gene symbolsをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #遺伝子IDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample18_5vs5.txt"            #入力ファイルを指定してin_fに格納
out_f <- "hoge5.txt"                   #出力ファイル名を指定してout_fに格納

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, sep="\t", quote="")#in_fで指定したファイルの読み込み

#genenameに相当する4列目の情報を抽出して加工("-ex"よりも左側の文字列のみ抽出)
hoge <- strsplit(as.character(data[,4]), "-ex", fixed=TRUE)#data[,param]を文字列に変換し、"-ex"で区切った結果をリスト形式でhogeに格納
genename <- unlist(lapply(hoge, "[[", 1))#hogeのリスト中の1番目の要素("-ex"で区切った左側部分に相当)のみ抽出してgenenameに格納
unique_genename <- unique(genename)    #non-redundantなgenename情報を抽出し、unique_genenameに格納

#1,6列目の情報（chrとstrand）はそのまま、5, 7-16列の情報（lengthとカウントデータ）のみsumしたいので、それぞれをサブセットに分ける
sub1 <- data[,c(1,6)]                  #1,6列目の情報のみ抽出しsub1に格納
sub2 <- data[,c(5,7:16)]               #5,7-16列目の情報のみ抽出しsub2に格納

out <- NULL                            #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(unique_genename)){   #unique_genenameの要素数分だけループを回す
    out_sub1 <- apply(sub1[which(genename == unique_genename[i]),], 2, unique, na.rm=TRUE)#sub1のところは、複数エクソンの場合は同じ情報がエクソン数分だけあることになるので、unique関数を実行した結果をout_sub1に格納
    out_sub2 <- apply(sub2[which(genename == unique_genename[i]),], 2, sum, na.rm=TRUE)#sub2のところでは、複数エクソンの場合にsum関数を実行した結果（和をとった結果）をout_sub2に格納
    out <- rbind(out, c(out_sub1, unique_genename[i], out_sub2))#「out_sub1, unique_genename[i], out_sub2」の順番で行列outの下に結果をどんどん追加
}
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample18_5vs5.txt"            #入力ファイルを指定してin_fに格納
out_f <- "hoge6.txt"                   #出力ファイル名を指定してout_fに格納

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, sep="\t", quote="")#in_fで指定したファイルの読み込み

#genenameに相当する4列目の情報を抽出して加工("-ex"よりも左側の文字列のみ抽出)
hoge <- strsplit(as.character(data[,4]), "-ex", fixed=TRUE)#data[,param]を文字列に変換し、"-ex"で区切った結果をリスト形式でhogeに格納
genename <- unlist(lapply(hoge, "[[", 1))#hogeのリスト中の1番目の要素("-ex"で区切った左側部分に相当)のみ抽出してgenenameに格納
unique_genename <- unique(genename)    #non-redundantなgenename情報を抽出し、unique_genenameに格納

#1,6列目の情報（chrとstrand）はそのまま、5, 7-16列の情報（lengthとカウントデータ）のみsumしたいので、それぞれをサブセットに分ける
sub1 <- data[,c(1,6)]                  #1,6列目の情報のみ抽出しsub1に格納
sub2 <- data[,c(5,7:16)]               #5,7-16列目の情報のみ抽出しsub2に格納

out <- NULL                            #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(unique_genename)){   #unique_genenameの要素数分だけループを回す
    tmp <- unlist(apply(sub1[which(genename == unique_genename[i]),], 2, unique, na.rm=TRUE))#sub1のところは、複数エクソンの場合は同じ情報がエクソン数分だけあることになるはずであるが、そうでない可能性があるときに見つけられるようにしている
    out_flag <- length(tmp)            #ベクトルtmpの要素数をout_flagに格納（通常はtmpの要素数が2だが、3以上のものを検出するのが目的）
    out_sub1 <- tmp[1:2]               #どんな状況になっていようと、とにかくベクトルtmpの最初の二つの要素を出力すべくout_sub1に格納
    out_sub2 <- apply(sub2[which(genename == unique_genename[i]),], 2, sum, na.rm=TRUE)#sub2のところでは、複数エクソンの場合にsum関数を実行した結果（和をとった結果）をout_sub2に格納
    out <- rbind(out, c(out_flag, out_sub1, unique_genename[i], out_sub2))#「out_sub1, unique_genename[i], out_sub2」の順番で行列outの下に結果をどんどん追加
}
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2.txt"               #入力ファイル名を指定してin_fに格納
out_f <- "hoge7.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定
param1 <- "rat2302.db"                 #パッケージ名を指定
param2 <- "ENTREZID"                   #欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

#必要なパッケージのインストール(利用したいパッケージが既に存在していれば2回目以降は必要なし)
source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
biocLite(param1, suppressUpdates=TRUE) #指定したパッケージのインストール

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#probe IDとGene symbolとの対応関係を含む情報を入手
hoge <- sub(".db", param2, param1)     #用いたいオブジェクト名情報を取得(rat2302SYMBOLに相当)
hoge <- eval(parse(text=hoge))         #hogeを文字列としてではなくオブジェクトとして取り扱う
symbols <- unlist(as.list(hoge))       #全31099 probesに対応するGeneSymbol(SYMBOL)情報を抽出し、symbolsに格納
unique_sym <- unique(symbols)          #non-redundantなGeneSymbol(SYMBOL)情報を抽出し、unique_symに格納
unique_sym <- unique_sym[unique_sym != ""]#unique_symの中から、Gene symbolがないものを除く
unique_sym <- unique_sym[!is.na(unique_sym)]#unique_symの中から、Gene symbolが"NA"のものを除く
unique_sym <- unique_sym[!is.nan(unique_sym)]#unique_symの中から、Gene symbolが"NaN"のものを除く

#本番
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(unique_sym)){        #non-redundant gene symbol数分だけループを回す
    hoge <- rbind(hoge, apply(data[which(symbols == unique_sym[i]),], 2, param, na.rm=TRUE))#dataの中からi番目のgene symbolと同じprobesを全て抽出し、その平均値(mean)をhogeの一番下の行に追加
}                                      #non-redundant gene symbol数分だけループを回す
rownames(hoge) <- unique_sym           #non-redundant gene symbolsをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #遺伝子IDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample5.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge8.txt"                   #出力ファイル名を指定してout_fに格納
param <- mean                          #代表値を指定
param1 <- "hgu133a.db"                 #パッケージ名を指定
param2 <- "ENTREZID"                   #欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

#必要なパッケージのインストール(利用したいパッケージが既に存在していれば2回目以降は必要なし)
source("http://bioconductor.org/biocLite.R")#指定したパッケージのインストール
biocLite(param1, suppressUpdates=TRUE) #指定したパッケージのインストール

#必要なパッケージをロード
library(param1, character.only=T)      #param1で指定したパッケージの読み込み

#probe IDとGene symbolとの対応関係を含む情報を入手
hoge <- sub(".db", param2, param1)     #用いたいオブジェクト名情報を取得(rat2302SYMBOLに相当)
hoge <- eval(parse(text=hoge))         #hogeを文字列としてではなくオブジェクトとして取り扱う
symbols <- unlist(as.list(hoge))       #全31099 probesに対応するGeneSymbol(SYMBOL)情報を抽出し、symbolsに格納
unique_sym <- unique(symbols)          #non-redundantなGeneSymbol(SYMBOL)情報を抽出し、unique_symに格納
unique_sym <- unique_sym[unique_sym != ""]#unique_symの中から、Gene symbolがないものを除く
unique_sym <- unique_sym[!is.na(unique_sym)]#unique_symの中から、Gene symbolが"NA"のものを除く
unique_sym <- unique_sym[!is.nan(unique_sym)]#unique_symの中から、Gene symbolが"NaN"のものを除く

#本番
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
hoge <- NULL                           #最終的に欲しい情報を格納するためのプレースホルダ
for(i in 1:length(unique_sym)){        #non-redundant gene symbol数分だけループを回す
    hoge <- rbind(hoge, apply(data[which(symbols == unique_sym[i]),], 2, param, na.rm=TRUE))#dataの中からi番目のgene symbolと同じprobesを全て抽出し、その平均値(mean)をhogeの一番下の行に追加
}                                      #non-redundant gene symbol数分だけループを回す
rownames(hoge) <- unique_sym           #non-redundant gene symbolsをhogeの行の名前として利用

#ファイルに保存
tmp <- cbind(rownames(hoge), hoge)     #遺伝子IDの列を行列hogeの左端に挿入し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f1 <- "genelist_A.txt"              #入力ファイル名を指定してin_f1に格納
in_f2 <- "genelist_B.txt"              #入力ファイル名を指定してin_f2に格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納

#データファイルの読み込み
data1 <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")#in_f1で指定したファイルの読み込み
data2 <- read.table(in_f2, header=TRUE, row.names=1, sep="\t", quote="")#in_f2で指定したファイルの読み込み

#本番
common <- intersect(rownames(data1), rownames(data2))#二つの遺伝子名のベクトル同士の積集合(intersection)をcommonに格納
obj <- is.element(rownames(data1), common)#rownames(data1)で表されるベクトル中の各要素がベクトルcommon中に含まれるか否かの情報をobjに格納（data1ではなくdata2のほうでもよい）
out <- data1[obj,]                     #objがTRUEとなる行のみ抽出した結果をoutに格納
names(out) <- rownames(data1)[obj]     #行列outに行名を付加している

#ファイルに保存
tmp <- cbind(names(out), out)          #行の名前、outを列ベクトル単位で結合し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample19.txt"                 #入力ファイル名を指定してin_fに格納

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
dist(t(data), method="euclidean")      #ユークリッド(Euclidean)距離
dist(t(data), method="manhattan")      #マンハッタン(Manhattan)距離
dist(t(data), method="maximum")        #チェビシェフ(Chebyshev)距離
dist(t(data), method="canberra")       #キャンベラ(Canberra)距離
1 - cor(data, method="pearson")        #1 - Pearson相関係数

dist(t(data), method="binary")         #ハミング(Hamming)距離
dist(t(data), method="minkowski")      #ミンコフスキー(Minkowski)距離
1 - cor(data, method="spearman")       #1 - Spearman相関係数
	

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

out <- apply(data, 1, mean)            #各行の平均値（mean）をoutに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, out)#行の名前、data、outを列ベクトル単位で結合し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

out <- apply(data, 1, median)          #各行の中央値（median）をoutに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, out)#行の名前、data、outを列ベクトル単位で結合し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納

#必要なパッケージをロード
library(affy)                          #tukey.biweight関数が含まれているaffyパッケージの読み込み

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

out <- apply(data, 1, tukey.biweight)  #各行のtukey.biweight値をoutに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, out)#行の名前、data、outを列ベクトル単位で結合し、結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
out <- colnames(data)[max.col(data)]   #最大発現量を示す組織名outに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, out)#入力データの右側にoutの情報を結合した結果をtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
out <- colnames(data)[max.col(data)]   #最大発現量を示す組織名outに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, out)#入力データの右側にoutの情報を結合した結果をtmpに格納。
tmp2 <- tmp[order(max.col(data)),]     #最大発現量を示す組織のシリアル番号順にソートした結果をtmp2に格納
write.table(tmp2, out_f, sep="\t", append=F, quote=F, row.names=F)#tmp2の中身を指定したファイル名で保存
    

in_f1 <- "sample5.txt"                 #入力ファイル名(発現データ)を指定してin_f1に格納
in_f2 <- "GDS1096_cl_heart.txt"        #入力ファイル名(テンプレート情報)を指定してin_f2に格納
out_f1 <- "data_torranded.txt"         #出力ファイル名(発現データ含むほう)を指定
out_f2 <- "data_topranked_ID.txt"      #出力ファイル名(遺伝子IDのみのほう)を指定
param1 <- 10                           #上位X個のXを指定
param2 <- "none"                       #類似度計算前の発現データのスケーリング法を指定
param3 <- "correlation"                #距離を定義する方法を指定

#必要なパッケージをロード
library(genefilter)                    #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")#in_f1で指定したファイルの読み込み
data <- as.matrix(data)                #as.matrixの意味は、「データの型を"行列として(as matrix)"dataに格納せよ」です。（read.tableで読み込んで得られたdataの型は"データフレーム"のため）
data_cl <- read.table(in_f2, sep="\t", quote="")#in_f2で指定したファイルの読み込み
template <- data_cl[,2]                #バイナリ（0 or 1）情報（2列目）のみ抽出し、templateに格納
template                               #バイナリ（0 or 1）情報の確認

#本番
tmp <- rbind(data, template)           #templateというテンプレートパターンを行列dataの最後の行に追加
template_posi <- which(rownames(tmp) == "template")#行のラベル情報が"template"に相当する行番号をtemplate_posiに格納
closeg <- genefinder(tmp, template_posi, param1, scale=param2, method=param3)#上位"param1"個の情報をclosegに格納
closeg[[1]]$indices                    #上位"param1"個の行番号を表示
closeg[[1]]$dists                      #上位"param1"個の類似度を表示
topranked <- tmp[closeg[[1]]$indices,] #上位"param1"個の遺伝子発現データを抽出し、toprankedに格納

#ファイルに保存
tmp2 <- cbind(rownames(topranked), topranked)#遺伝子IDの列を行列toprankedの左端に挿入し、結果をtmp2に格納
write.table(tmp2, out_f1, sep="\t", append=F, quote=F, row.names=F)#tmp2の中身を指定したファイル名で保存
write.table(rownames(topranked), out_f2, sep="\t", append=F, quote=F, row.names=F, col.names=F)#遺伝子IDに関する情報のみ、指定したファイル名で保存
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f1 <- "sample5.txt"                 #入力ファイル名(発現データ)を指定してin_f1に格納
in_f2 <- "GDS1096_cl_heart.txt"        #入力ファイル名(テンプレート情報)を指定してin_f2に格納
out_f1 <- "data_torranded.txt"         #出力ファイル名(発現データ含むほう)を指定
out_f2 <- "data_topranked_ID.txt"      #出力ファイル名(遺伝子IDのみのほう)を指定
param2 <- "zscore"                     #類似度計算前の発現データのスケーリング法を指定
param3 <- "correlation"                #距離を定義する方法を指定

#必要なパッケージをロード
library(genefilter)                    #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")#入力ファイル1を読み込んでdataに格納
data <- as.matrix(data)                #as.matrixの意味は、「データの型を"行列として(as matrix)"dataに格納せよ」です。（read.tableで読み込んで得られたdataの型は"データフレーム"のため）
data_cl <- read.table(in_f2, sep="\t", quote="")#入力ファイル2を読み込んでdata_clに格納
template <- data_cl[,2]                #バイナリ（0 or 1）情報（2列目）のみ抽出し、templateに格納
template                               #バイナリ（0 or 1）情報の確認

#本番
tmp <- rbind(data, template)           #templateというテンプレートパターンを行列dataの最後の行に追加
template_posi <- which(rownames(tmp) == "template")#行のラベル情報が"template"に相当する行番号をtemplate_posiに格納
param1 <- nrow(data)                   #遺伝子数をparam1に格納
closeg <- genefinder(tmp, template_posi, param1, scale=param2, method=param3)#特異的発現の度合いでランキングされた結果をclosegに格納
topranked <- tmp[closeg[[1]]$indices,] #特異的発現の度合いでランキングされた遺伝子発現データをtoprankedに格納

#ファイルに保存
tmp2 <- cbind(rownames(topranked), topranked)#遺伝子IDの列を行列toprankedの左端に挿入し、結果をtmp2に格納
write.table(tmp2, out_f1, sep="\t", append=F, quote=F, row.names=F)#tmp2の中身を指定したファイル名で保存
write.table(rownames(topranked), out_f2, sep="\t", append=F, quote=F, row.names=F, col.names=F)#遺伝子IDに関する情報のみ、指定したファイル名で保存
	

in_f <- "sample5.txt"                  #入力ファイル名を指定してin_fに格納
out_f1 <- "data_torranded.txt"         #出力ファイル名(発現データ含むほう)を指定
out_f2 <- "data_topranked_ID.txt"      #出力ファイル名(遺伝子IDのみのほう)を指定
param1 <- 5                            #上位X個のXを指定
param2 <- "range"                      #類似度計算前の発現データのスケーリング法を指定
param3 <- "manhattan"                  #距離を定義する方法を指定
param4 <- "207003_at"                  #遺伝子IDを指定

#必要なパッケージをロード
library(genefilter)                    #パッケージの読み込み

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data <- as.matrix(data)                #as.matrixの意味は、「データの型を"行列として(as matrix)"dataに格納せよ」です。（read.tableで読み込んで得られたdataの型は"データフレーム"のため）

#本番
template_posi <- which(rownames(data) == param4)#param4で指定した遺伝子IDに相当する行番号をtemplate_posiに格納
closeg <- genefinder(data, template_posi, param1, scale=param2, method=param3)#上位"param1"個の情報をclosegに格納
topranked <- data[closeg[[1]]$indices,]#上位"param1"個の遺伝子発現データを抽出し、toprankedに格納

#ファイルに保存
tmp2 <- cbind(rownames(topranked), topranked)#遺伝子IDの列を行列toprankedの左端に挿入し、結果をtmp2に格納
write.table(tmp2, out_f1, sep="\t", append=F, quote=F, row.names=F)#tmp2の中身を指定したファイル名で保存
write.table(rownames(topranked), out_f2, sep="\t", append=F, quote=F, row.names=F, col.names=F)#遺伝子IDに関する情報のみ、指定したファイル名で保存
	

in_f <- "sample5.txt"                  #入力ファイル名を指定してin_fに格納
out_f1 <- "data_torranded.txt"         #出力ファイル名(発現データ含むほう)を指定してout_f1に格納
out_f2 <- "data_topranked_ID.txt"      #出力ファイル名(遺伝子IDのみのほう)を指定してout_f2に格納
param1 <- 10                           #上位X個のXを指定
param2 <- "zscore"                     #類似度計算前の発現データのスケーリング法を指定
param3 <- "correlation"                #距離を定義する方法を指定
param4 <- 15987                        #目的遺伝子の行番号を指定

#必要なパッケージをロード
library(genefilter)                    #パッケージの読み込み

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data <- as.matrix(data)                #as.matrixの意味は、「データの型を"行列として(as matrix)"dataに格納せよ」です。（read.tableで読み込んで得られたdataの型は"データフレーム"のため）

#本番
closeg <- genefinder(data, param4, param1, scale=param2, method=param3)#上位"param1"個の情報をclosegに格納
topranked <- data[closeg[[1]]$indices,]#上位"param1"個の遺伝子発現データを抽出し、toprankedに格納

#ファイルに保存
tmp2 <- cbind(rownames(topranked), topranked)#遺伝子IDの列を行列toprankedの左端に挿入し、結果をtmp2に格納
write.table(tmp2, out_f1, sep="\t", append=F, quote=F, row.names=F)#tmp2の中身を指定したファイル名で保存
write.table(rownames(topranked), out_f2, sep="\t", append=F, quote=F, row.names=F, col.names=F)#遺伝子IDに関する情報のみ、指定したファイル名で保存
	

in_f <- "data_GSE30533_mas.txt"        #入力ファイル名を指定してin_fに格納
out_f1 <- "hoge1_G1.txt"               #出力ファイル名を指定してout_f1に格納
out_f2 <- "hoge1_G1.png"               #出力ファイル名を指定してout_f2に格納
out_f3 <- "hoge1_G2.txt"               #出力ファイル名を指定してout_f3に格納
out_f4 <- "hoge1_G2.png"               #出力ファイル名を指定してout_f4に格納
out_f5 <- "hoge1_all.png"              #出力ファイル名を指定してout_f5に格納
param_G1 <- 5                          #G1群のサンプル数を指定
param_G2 <- 5                          #G2群のサンプル数を指定
param_fig <- c(380, 420)               #ファイル出力時の横幅と縦幅を指定(単位はピクセル)
param_range <- c(1e-02, 1e+08)         #M-A plot上での軸の数値範囲を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#前処理(対数変換後のデータなので、オリジナルスケールに変換したのち、群ごとに分割)
data <- 2^data                         #データ変換
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成
G1 <- data[,data.cl==1]                #G1群のデータのみ抽出している
G2 <- data[,data.cl==2]                #G2群のデータのみ抽出している

#本番(Mean-Variance plotなど; G1群)
hoge <- G1                             #G1オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
tmp <- cbind(rownames(data), data[,data.cl==1], hoge, MEAN, VARIANCE)#保存したい情報をtmpに格納
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,#プロット
     xlim=param_range, ylim=param_range, col="blue")#プロット
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
abline(a=0, b=1, col="gray")           #y=xの直線を指定した色で追加（y=a+bxのa=0, b=1）
legend("topright", "G1", col="blue", pch=20)#凡例を作成している

hoge <- hoge[apply(hoge, 1, var) > 0,] #回帰分析のときにエラーが出ないように分散>0のもののみ抽出しているだけ
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
hoge <- as.data.frame(cbind(MEAN, VARIANCE))#回帰分析（regression analysis）を行うためののおまじない（1列目がMEAN, 2列目がVARIANCEならなる行列を作成したあとデータフレーム形式にした結果をhogeに格納）
out <- lm(VARIANCE~MEAN, data=log10(hoge))#独立変数（説明変数）をMEAN, 従属変数（目的変数）をVARIANCEとしてlog10変換したデータの線形回帰を行った結果をoutに格納
abline(out, col="black")               #回帰直線を追加
out                                    #outの簡単な中身を表示（切片(Intercept)が-0.001173, 傾き(MEAN)が1.349519であることがわかる。つまり、y=a+bxのaが切片、bが傾きに相当する）
summary(out)                           #回帰分析結果outのもう少し詳細な結果を表示している（Multiple R-squared（決定係数）の値が0.8149と1に相当近い値が得られていることから線形回帰で十分よいfittingが得られていると判断できる。また、一番下のp-valueが限りなく0に近いことは、MEANという従属変数が不要であるという帰無仮説を棄却するに値する、つまり従属変数が独立変数によって説明可能であることを意味する）
dev.off()                              #おまじない

#本番(Mean-Variance plotなど; G2群)
hoge <- G2                             #G2オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
tmp <- cbind(rownames(data), data[,data.cl==1], hoge, MEAN, VARIANCE)#保存したい情報をtmpに格納
write.table(tmp, out_f3, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

png(out_f4, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,#プロット
     xlim=param_range, ylim=param_range, col="red")#プロット
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
abline(a=0, b=1, col="gray")           #y=xの直線を指定した色で追加（y=a+bxのa=0, b=1）
legend("topright", "G2", col="red", pch=20)#凡例を作成している

hoge <- hoge[apply(hoge, 1, var) > 0,] #回帰分析のときにエラーが出ないように分散>0のもののみ抽出しているだけ
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
hoge <- as.data.frame(cbind(MEAN, VARIANCE))#回帰分析（regression analysis）を行うためののおまじない（1列目がMEAN, 2列目がVARIANCEならなる行列を作成したあとデータフレーム形式にした結果をhogeに格納）
out <- lm(VARIANCE~MEAN, data=log10(hoge))#独立変数（説明変数）をMEAN, 従属変数（目的変数）をVARIANCEとしてlog10変換したデータの線形回帰を行った結果をoutに格納
abline(out, col="black")               #回帰直線を追加
out                                    #outの簡単な中身を表示（切片(Intercept)が-0.001173, 傾き(MEAN)が1.349519であることがわかる。つまり、y=a+bxのaが切片、bが傾きに相当する）
summary(out)                           #回帰分析結果outのもう少し詳細な結果を表示している（Multiple R-squared（決定係数）の値が0.8149と1に相当近い値が得られていることから線形回帰で十分よいfittingが得られていると判断できる。また、一番下のp-valueが限りなく0に近いことは、MEANという従属変数が不要であるという帰無仮説を棄却するに値する、つまり従属変数が独立変数によって説明可能であることを意味する）
dev.off()                              #おまじない

#本番(Mean-Variance plot; G1 and G2両方)
hoge <- G1                             #G1オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
png(out_f5, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1, ann=F,#プロット
     xlim=param_range, ylim=param_range, col="blue")#プロット

par(new=T)                             #図の重ね合わせをするという宣言
hoge <- G2                             #G2オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,#プロット
     xlim=param_range, ylim=param_range, col="red")#プロット
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
abline(a=0, b=1, col="gray")           #y=xの直線を指定した色で追加（y=a+bxのa=0, b=1）
legend("topright", c("G1", "G2"), col=c("blue", "red"), pch=20)#凡例を作成している
dev.off()                              #おまじない
	

in_f <- "data_rma_2.txt"               #入力ファイル名を指定してin_fに格納
out_f1 <- "hoge2_G1.txt"               #出力ファイル名を指定してout_f1に格納
out_f2 <- "hoge2_G1.png"               #出力ファイル名を指定してout_f2に格納
out_f3 <- "hoge2_G2.txt"               #出力ファイル名を指定してout_f3に格納
out_f4 <- "hoge2_G2.png"               #出力ファイル名を指定してout_f4に格納
out_f5 <- "hoge2_all.png"              #出力ファイル名を指定してout_f5に格納
param_G1 <- 4                          #G1群(fed)のサンプル数を指定
param_G2 <- 4                          #G2群(fasted)のサンプル数を指定
param_fig <- c(380, 420)               #ファイル出力時の横幅と縦幅を指定(単位はピクセル)
param_range <- c(1e-02, 1e+08)         #M-A plot上での軸の数値範囲を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
colSums(data)                          #総リード数を表示

#前処理(対数変換後のデータなので、オリジナルスケールに変換している)
data <- 2^data                         #データ変換

#前処理(最初の8列分のみ抽出したのち、群ごとのサブセットに分けている)
data <- data[,1:8]                     #最初の8列分のデータのみ抽出している
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成
G1 <- data[,data.cl==1]                #G1群のデータのみ抽出している
colSums(G1)                            #総リード数を表示
G2 <- data[,data.cl==2]                #G2群のデータのみ抽出している
colSums(G2)                            #総リード数を表示

#本番(Mean-Variance plotなど; G1群)
hoge <- G1                             #G1オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
tmp <- cbind(rownames(data), data[,data.cl==1], hoge, MEAN, VARIANCE)#保存したい情報をtmpに格納
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,#プロット
     xlim=param_range, ylim=param_range, col="blue")#プロット
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
abline(a=0, b=1, col="gray")           #y=xの直線を指定した色で追加（y=a+bxのa=0, b=1）
legend("topright", "G1", col="blue", pch=20)#凡例を作成している

hoge <- hoge[apply(hoge, 1, var) > 0,] #回帰分析のときにエラーが出ないように分散>0のもののみ抽出しているだけ
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
hoge <- as.data.frame(cbind(MEAN, VARIANCE))#回帰分析（regression analysis）を行うためののおまじない（1列目がMEAN, 2列目がVARIANCEならなる行列を作成したあとデータフレーム形式にした結果をhogeに格納）
out <- lm(VARIANCE~MEAN, data=log10(hoge))#独立変数（説明変数）をMEAN, 従属変数（目的変数）をVARIANCEとしてlog10変換したデータの線形回帰を行った結果をoutに格納
abline(out, col="black")               #回帰直線を追加
out                                    #outの簡単な中身を表示（切片(Intercept)が-0.001173, 傾き(MEAN)が1.349519であることがわかる。つまり、y=a+bxのaが切片、bが傾きに相当する）
summary(out)                           #回帰分析結果outのもう少し詳細な結果を表示している（Multiple R-squared（決定係数）の値が0.8149と1に相当近い値が得られていることから線形回帰で十分よいfittingが得られていると判断できる。また、一番下のp-valueが限りなく0に近いことは、MEANという従属変数が不要であるという帰無仮説を棄却するに値する、つまり従属変数が独立変数によって説明可能であることを意味する）
dev.off()                              #おまじない

#本番(Mean-Variance plotなど; G2群)
hoge <- G2                             #G2オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
tmp <- cbind(rownames(data), data[,data.cl==1], hoge, MEAN, VARIANCE)#保存したい情報をtmpに格納
write.table(tmp, out_f3, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

png(out_f4, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,#プロット
     xlim=param_range, ylim=param_range, col="red")#プロット
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
abline(a=0, b=1, col="gray")           #y=xの直線を指定した色で追加（y=a+bxのa=0, b=1）
legend("topright", "G2", col="red", pch=20)#凡例を作成している

hoge <- hoge[apply(hoge, 1, var) > 0,] #回帰分析のときにエラーが出ないように分散>0のもののみ抽出しているだけ
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
hoge <- as.data.frame(cbind(MEAN, VARIANCE))#回帰分析（regression analysis）を行うためののおまじない（1列目がMEAN, 2列目がVARIANCEならなる行列を作成したあとデータフレーム形式にした結果をhogeに格納）
out <- lm(VARIANCE~MEAN, data=log10(hoge))#独立変数（説明変数）をMEAN, 従属変数（目的変数）をVARIANCEとしてlog10変換したデータの線形回帰を行った結果をoutに格納
abline(out, col="black")               #回帰直線を追加
out                                    #outの簡単な中身を表示（切片(Intercept)が-0.001173, 傾き(MEAN)が1.349519であることがわかる。つまり、y=a+bxのaが切片、bが傾きに相当する）
summary(out)                           #回帰分析結果outのもう少し詳細な結果を表示している（Multiple R-squared（決定係数）の値が0.8149と1に相当近い値が得られていることから線形回帰で十分よいfittingが得られていると判断できる。また、一番下のp-valueが限りなく0に近いことは、MEANという従属変数が不要であるという帰無仮説を棄却するに値する、つまり従属変数が独立変数によって説明可能であることを意味する）
dev.off()                              #おまじない

#本番(Mean-Variance plot; G1 and G2両方)
hoge <- G1                             #G1オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
png(out_f5, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1, ann=F,#プロット
     xlim=param_range, ylim=param_range, col="blue")#プロット

par(new=T)                             #図の重ね合わせをするという宣言
hoge <- G2                             #G2オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,#プロット
     xlim=param_range, ylim=param_range, col="red")#プロット
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
abline(a=0, b=1, col="gray")           #y=xの直線を指定した色で追加（y=a+bxのa=0, b=1）
legend("topright", c("G1", "G2"), col=c("blue", "red"), pch=20)#凡例を作成している
dev.off()                              #おまじない
	

in_f <- "data_mas.txt"                 #入力ファイル名を指定してin_fに格納
out_f1 <- "hoge3_G1.txt"               #出力ファイル名を指定してout_f1に格納
out_f2 <- "hoge3_G1.png"               #出力ファイル名を指定してout_f2に格納
out_f3 <- "hoge3_G2.txt"               #出力ファイル名を指定してout_f3に格納
out_f4 <- "hoge3_G2.png"               #出力ファイル名を指定してout_f4に格納
out_f5 <- "hoge3_all.png"              #出力ファイル名を指定してout_f5に格納
param_G1 <- 4                          #G1群(fed)のサンプル数を指定
param_G2 <- 4                          #G2群(fasted)のサンプル数を指定
param_fig <- c(380, 420)               #ファイル出力時の横幅と縦幅を指定(単位はピクセル)
param_range <- c(1e-02, 1e+08)         #M-A plot上での軸の数値範囲を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
colSums(data)                          #総リード数を表示

#前処理(対数変換後のデータなので、オリジナルスケールに変換している)
data <- 2^data                         #データ変換

#前処理(最初の8列分のみ抽出したのち、群ごとのサブセットに分けている)
data <- data[,1:8]                     #最初の8列分のデータのみ抽出している
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成
G1 <- data[,data.cl==1]                #G1群のデータのみ抽出している
colSums(G1)                            #総リード数を表示
G2 <- data[,data.cl==2]                #G2群のデータのみ抽出している
colSums(G2)                            #総リード数を表示

#本番(Mean-Variance plotなど; G1群)
hoge <- G1                             #G1オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
tmp <- cbind(rownames(data), data[,data.cl==1], hoge, MEAN, VARIANCE)#保存したい情報をtmpに格納
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,#プロット
     xlim=param_range, ylim=param_range, col="blue")#プロット
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
abline(a=0, b=1, col="gray")           #y=xの直線を指定した色で追加（y=a+bxのa=0, b=1）
legend("topright", "G1", col="blue", pch=20)#凡例を作成している

hoge <- hoge[apply(hoge, 1, var) > 0,] #回帰分析のときにエラーが出ないように分散>0のもののみ抽出しているだけ
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
hoge <- as.data.frame(cbind(MEAN, VARIANCE))#回帰分析（regression analysis）を行うためののおまじない（1列目がMEAN, 2列目がVARIANCEならなる行列を作成したあとデータフレーム形式にした結果をhogeに格納）
out <- lm(VARIANCE~MEAN, data=log10(hoge))#独立変数（説明変数）をMEAN, 従属変数（目的変数）をVARIANCEとしてlog10変換したデータの線形回帰を行った結果をoutに格納
abline(out, col="black")               #回帰直線を追加
out                                    #outの簡単な中身を表示（切片(Intercept)が-0.001173, 傾き(MEAN)が1.349519であることがわかる。つまり、y=a+bxのaが切片、bが傾きに相当する）
summary(out)                           #回帰分析結果outのもう少し詳細な結果を表示している（Multiple R-squared（決定係数）の値が0.8149と1に相当近い値が得られていることから線形回帰で十分よいfittingが得られていると判断できる。また、一番下のp-valueが限りなく0に近いことは、MEANという従属変数が不要であるという帰無仮説を棄却するに値する、つまり従属変数が独立変数によって説明可能であることを意味する）
dev.off()                              #おまじない

#本番(Mean-Variance plotなど; G2群)
hoge <- G2                             #G2オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
tmp <- cbind(rownames(data), data[,data.cl==1], hoge, MEAN, VARIANCE)#保存したい情報をtmpに格納
write.table(tmp, out_f3, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

png(out_f4, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,#プロット
     xlim=param_range, ylim=param_range, col="red")#プロット
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
abline(a=0, b=1, col="gray")           #y=xの直線を指定した色で追加（y=a+bxのa=0, b=1）
legend("topright", "G2", col="red", pch=20)#凡例を作成している

hoge <- hoge[apply(hoge, 1, var) > 0,] #回帰分析のときにエラーが出ないように分散>0のもののみ抽出しているだけ
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
hoge <- as.data.frame(cbind(MEAN, VARIANCE))#回帰分析（regression analysis）を行うためののおまじない（1列目がMEAN, 2列目がVARIANCEならなる行列を作成したあとデータフレーム形式にした結果をhogeに格納）
out <- lm(VARIANCE~MEAN, data=log10(hoge))#独立変数（説明変数）をMEAN, 従属変数（目的変数）をVARIANCEとしてlog10変換したデータの線形回帰を行った結果をoutに格納
abline(out, col="black")               #回帰直線を追加
out                                    #outの簡単な中身を表示（切片(Intercept)が-0.001173, 傾き(MEAN)が1.349519であることがわかる。つまり、y=a+bxのaが切片、bが傾きに相当する）
summary(out)                           #回帰分析結果outのもう少し詳細な結果を表示している（Multiple R-squared（決定係数）の値が0.8149と1に相当近い値が得られていることから線形回帰で十分よいfittingが得られていると判断できる。また、一番下のp-valueが限りなく0に近いことは、MEANという従属変数が不要であるという帰無仮説を棄却するに値する、つまり従属変数が独立変数によって説明可能であることを意味する）
dev.off()                              #おまじない

#本番(Mean-Variance plot; G1 and G2両方)
hoge <- G1                             #G1オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
png(out_f5, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1, ann=F,#プロット
     xlim=param_range, ylim=param_range, col="blue")#プロット

par(new=T)                             #図の重ね合わせをするという宣言
hoge <- G2                             #G2オブジェクトをhogeに格納
MEAN <- apply(hoge, 1, mean)           #各行の平均を計算した結果をMEANに格納
VARIANCE <- apply(hoge, 1, var)        #各行の分散を計算した結果をVARIANCEに格納
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,#プロット
     xlim=param_range, ylim=param_range, col="red")#プロット
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
abline(a=0, b=1, col="gray")           #y=xの直線を指定した色で追加（y=a+bxのa=0, b=1）
legend("topright", c("G1", "G2"), col=c("blue", "red"), pch=20)#凡例を作成している
dev.off()                              #おまじない
	

in_f <- "sample5.txt"                  #入力ファイル名(発現データファイル)を指定してin_fに格納
param1 <- "correlation"                #類似度(method.dist)を指定
param2 <- "average"                    #方法(method.hclust)を指定
param3 <- 20                           #resampling回数を指定（この値が大きいほどより正確にブートストラップ確率を求めることができる。実際には100とか500とか...）

#入力ファイルの読み込みとラベル情報の作成
library(pvclust)                       #パッケージの読み込み

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
out <- pvclust(data, method.hclust=param2, method.dist=param1, nboot=param3)#クラスタリングの実行
plot(out)                              #樹形図（デンドログラム）の表示

#以下は（こんなこともできますという）おまけ
#Approximately Unbiased (au) probability > 0.95を満たす頑健なクラスターを四角で囲み、その条件を満たすクラスターのメンバーをリストアップしたい場合：
param4 <- "au"                         #ブートストラップ確率計算手段を指定
param5 <- 0.95                         #ブートストラップ確率の閾値を指定
pvrect(out, alpha=param5, pv=param4)   #条件を満たすクラスターを四角で囲む
pvpick(out, alpha=param5, pv=param4)   #条件を満たすクラスターのメンバーをリストアップする
	

in_f <- "sample3.txt"                  #入力ファイル名(発現データファイル)を指定してin_fに格納

#必要なパッケージをロード
library(pvclust)                       #パッケージの読み込み

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
out <- pvclust(data)                   #クラスタリングの実行plot(out)#樹形図（デンドログラム）の表示
	

in_f <- "sample3.txt"                  #入力ファイル名(発現データファイル)を指定してin_fに格納
param1 <- "correlation"                #類似度(method.dist)を指定
param2 <- "average"                    #方法(method.hclust)を指定
param3 <- 30                           #resampling回数を指定（この値が大きいほどより正確にブートストラップ確率を求めることができる。実際には100とか500とか...）

#必要なパッケージをロード
library(pvclust)                       #パッケージの読み込み

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
out <- pvclust(t(data), method.hclust=param2, method.dist=param1, nboot=param3)#クラスタリングの実行
plot(out)                              #樹形図（デンドログラム）の表示
	

in_f <- "sample3.txt"                  #入力ファイル名を指定してin_fに格納
param <- "average"                     #方法(method)を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
data.dist <- as.dist(1 - cor(data, method="pearson"))#サンプル間の距離を計算した結果をdata.distに格納
out <- hclust(data.dist, method=param) #階層的クラスタリングを実行した結果をoutに格納
plot(out)                              #樹形図（デンドログラム）の表示
	

in_f <- "sample3.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.png"                   #出力ファイル名を指定してout_fに格納
param <- "average"                     #方法(method)を指定
param_fig <- c(500, 400)               #ファイル出力時の横幅と縦幅を指定(単位はピクセル)

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
data.dist <- as.dist(1 - cor(data, method="pearson"))#サンプル間の距離を計算した結果をdata.distに格納
out <- hclust(data.dist, method=param) #階層的クラスタリングを実行した結果をoutに格納

#ファイルに保存
png(out_f, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(out)                              #樹形図（デンドログラム）の表示
dev.off()                              #おまじない
	

in_f <- "sample3.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge3.png"                   #出力ファイル名を指定してout_fに格納
param <- "average"                     #方法(method)を指定
param_fig <- c(500, 400)               #ファイル出力時の横幅と縦幅を指定(単位はピクセル)

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
data.dist <- as.dist(1 - cor(data, method="spearman"))#サンプル間の距離を計算した結果をdata.distに格納
out <- hclust(data.dist, method=param) #階層的クラスタリングを実行した結果をoutに格納

#ファイルに保存
png(out_f, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(out)                              #樹形図（デンドログラム）の表示
dev.off()                              #おまじない
	

in_f <- "sample3.txt"                  #入力ファイル名を指定してin_fに格納
param1 <- "euclidean"                  #距離(dist)を指定
param2 <- "average"                    #方法(method)を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
data.dist <- dist(t(data), method=param1)#サンプル間の距離を計算した結果をdata.distに格納
out <- hclust(data.dist, method=param2)#階層的クラスタリングを実行した結果をoutに格納
plot(out)                              #樹形図（デンドログラム）の表示
	

in_f <- "sample3.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge5.png"                   #出力ファイル名を指定してout_fに格納
param1 <- "euclidean"                  #距離(dist)を指定
param2 <- "average"                    #方法(method)を指定
param_fig <- c(500, 400)               #ファイル出力時の横幅と縦幅を指定(単位はピクセル)

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
data.dist <- dist(t(data), method=param1)#サンプル間の距離を計算した結果をdata.distに格納
out <- hclust(data.dist, method=param2)#階層的クラスタリングを実行した結果をoutに格納

#ファイルに保存
png(out_f, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(out)                              #樹形図（デンドログラム）の表示
dev.off()                              #おまじない
	

in_f <- "sample3.txt"                  #入力ファイル名を指定してin_fに格納
param1 <- "euclidean"                  #距離(dist)を指定
param2 <- "average"                    #方法(method)を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
data.dist <- dist(data, method=param1) #遺伝子間の距離を計算した結果をdata.distに格納
out <- hclust(data.dist, method=param2)#階層的クラスタリングを実行した結果をoutに格納
plot(out)                              #樹形図（デンドログラム）の表示
	

in_f <- "sample3.txt"                  #入力ファイル名を指定してin_fに格納
param <- "average"                     #方法(method)を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#本番
data.dist <- as.dist(1 - cor(t(data), method="spearman")) #遺伝子間の距離を計算した結果をdata.distに格納
out <- hclust(data.dist, method=param) #階層的クラスタリングを実行した結果をoutに格納
plot(out)                              #樹形図（デンドログラム）の表示
	

in_f <- "data_Singh_RMA_3274.txt"      #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param1 <- 10                           #K個のクラスターに分割のKを指定
param2 <- "average"                    #方法(method)を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#まずはクラスタリングを実行
data.dist <- as.dist(1 - cor(t(data))) #遺伝子間の距離を計算し、結果をdata.distに格納
out <- hclust(data.dist, method=param2)#階層的クラスタリングを実行し、結果をoutに格納

#param1で指定した数に分割
cluster <- cutree(out, k=param1)       #param1個のクラスターに分割し、結果をclusterに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, cluster)#入力データの右側にclusterの情報を追加してtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

#以下は（こんなこともできますという）おまけ
#クラスターごとにソートした結果を得たい場合：
out_f2 <- "hoge2.txt"                  #出力ファイル名を指定してout_f2に格納
tmp2 <- tmp[order(cluster),]           #cluster列の値でソートした結果をtmp2に格納
write.table(tmp2, out_f2, sep="\t", append=F, quote=F, row.names=F)#tmp2の中身を指定したファイル名で保存
	

in_f <- "sample3.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param1 <- 3                            #K個のクラスターに分割のKを指定
param2 <- "euclidean"                  #距離(dist)を指定
param3 <- "average"                    #方法(method)を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#まずはクラスタリングを実行
data.dist <- dist(data, method=param2) #遺伝子間の距離を計算し、結果をdata.distに格納
out <- hclust(data.dist, method=param3)#階層的クラスタリングを実行し、結果をoutに格納

#param1で指定した数に分割
cluster <- cutree(out, k=param1)       #param1個のクラスターに分割し、結果をclusterに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, cluster)#入力データの右側にclusterの情報を追加してtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

#以下は（こんなこともできますという）おまけ
#クラスターごとにソートした結果を得たい場合：
out_f2 <- "hoge2.txt"                  #出力ファイル名を指定してout_f2に格納
tmp2 <- tmp[order(cluster),]           #cluster列の値でソートした結果をtmp2に格納
write.table(tmp2, out_f2, sep="\t", append=F, quote=F, row.names=F)#tmp2の中身を指定したファイル名で保存
	

in_f <- "data_GSE7623_rma.txt"         #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param1 <- 2:10                         #K個のクラスターに分割のKを指定(ここでは2-10を一度に指定←こんなこともできます)
param2 <- "average"                    #方法(method)を指定

#入力ファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
colnames(data) <- c(                   #サンプル名を変更している
"BAT_fed1", "BAT_fed2", "BAT_fed3", "BAT_fed4", #サンプル名を変更している
"BAT_fasted1", "BAT_fasted2", "BAT_fasted3", "BAT_fasted4", #サンプル名を変更している
"WAT_fed1", "WAT_fed2", "WAT_fed3", "WAT_fed4", #サンプル名を変更している
"WAT_fasted1", "WAT_fasted2", "WAT_fasted3", "WAT_fasted4", #サンプル名を変更している
"LIV_fed1", "LIV_fed2", "LIV_fed3", "LIV_fed4", #サンプル名を変更している
"LIV_fasted1", "LIV_fasted2", "LIV_fasted3", "LIV_fasted4")#サンプル名を変更している

#まずはクラスタリングを実行
data.dist <- as.dist(1 - cor(data))    #遺伝子間の距離を計算し、結果をdata.distに格納
out <- hclust(data.dist, method=param2)#階層的クラスタリングを実行し、結果をoutに格納

#param1で指定した数に分割
cluster <- cutree(out, k=param1)       #param1個のクラスターに分割し、結果をclusterに格納

#ファイルに保存
tmp <- cbind(colnames(data), cluster)  #サンプル名の右側にk個に分割した(例ではk=2,3,...,10)場合のclusterの情報を追加してtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納

#必要なパッケージをロード
library(som)                           #パッケージの読み込み

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#フィルタリング
sample2.f <- filtering(data, lt=10, ut=100, mmr=2.9, mmd=42)#解析に用いる遺伝子のフィルタリング(発現強度が10以下のものを10に、100以上のものを100に、「発現強度max/min < 2.9」または「(max-min) < 42」の行(遺伝子)を解析から除く)
sample2.f.n <- normalize(sample2.f, byrow=TRUE)#行(row)方向の正規化(平均=0,分散=1); (列方向にしたいときはcolにする)

#本番
foo <- som(t(sample2.f.n), xdim=3, ydim=5)#SOMの実行(「x軸方向3分割×y軸方向5分割」にする場合)
plot(foo)                              #結果のプロット
	

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納

#必要なパッケージをロード
library(som)                           #パッケージの読み込み

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み

#フィルタリング
sample2.f <- filtering(data, lt=10, ut=100, mmr=2.9, mmd=42)#解析に用いる遺伝子のフィルタリング(発現強度が10以下のものを10に、100以上のものを100に、「発現強度max/min < 2.9」または「(max-min) < 42」の行(遺伝子)を解析から除く)
sample2.f.n <- normalize(sample2.f, byrow=TRUE)#行(row)方向の正規化(平均=0,分散=1); (列方向にしたいときはcolにする)

#本番
foo <- som(sample2.f.n, xdim=3, ydim=5)#SOMの実行(「x軸方向3分割×y軸方向5分割」にする場合)
plot(foo)                              #結果のプロット
	

in_f <- "sample14.txt"                 #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param_G1 <- 3                          #G1群のサンプル数を指定
param_G2 <- 3                          #G2群のサンプル数を指定

#必要なパッケージをロード
library(OCplus)                        #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
out1 <- fdr1d(data, data.cl, verb=FALSE)#入力データ中の発現変動していない遺伝子(non-DEGs)の割合を調べた結果をout1に格納
out2 <- EOC(data, data.cl)             #入力データ中の発現変動していない遺伝子(non-DEGs)の割合を調べた結果をout2に格納
p0(out1)                               #得られたout1の中から目的のnon-DEGsの割合の数値を表示
p0(out2)                               #得られたout2の中から目的のnon-DEGsの割合の数値を表示
	

in_f <- "data_rma_2_BAT.txt"           #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 4                          #G1群のサンプル数を指定
param_G2 <- 4                          #G2群のサンプル数を指定

#必要なパッケージをロード
library(RankProd)                      #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(0, param_G1), rep(1, param_G2))#G1群を0、G2群を1としたベクトルdata.clを作成

#本番
out <- RankProducts(data, data.cl, logged=T, na.rm=F, plot=F, rand=123)#RankProductsの実行
stat_RP <- apply(out$RPs, 1, min)      #総合RP統計量を計算(20090527追加)
rank_RP <- as.matrix(rank(stat_RP, ties.method = "min"))#総合順位を計算(20090527追加)
stat_fc <- -out$AveFC                  #総合log(FC)統計量を計算
rank_fc <- rank(-abs(stat_fc))         #総合log(FC)統計量の順位を計算

#ファイルに保存
colnames(out$pfp) <- c("FDR(G1 < G2)","FDR(G1 > G2)")#列名を付与
colnames(out$RPs) <- c("stat(G1 < G2)","stat(G1 > G2)")#列名を付与
colnames(out$RPrank) <- c("rank(G1 < G2)","rank(G1 > G2)")#列名を付与
colnames(stat_fc) <- "log2(G2/G1)"     #列名を付与
tmp <- cbind(rownames(data), data, out$pfp, out$RPs, out$RPrank, stat_fc, rank_fc, stat_RP, rank_RP)#入力データの右側にFDR値, 各統計量, 各順位, 総合log(FC)統計量, その順位, 総合RP統計量, その順位の情報を付加したtmpを用意。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_Singh_RMA_3274.txt"      #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 50                         #G1群のサンプル数を指定
param_G2 <- 52                         #G2群のサンプル数を指定

#WAD統計量を計算するための関数を定義
WAD <- function(data=NULL, data.cl=NULL){#WAD統計量を計算するための関数
    x <- data                          #WAD統計量を計算するための関数
    cl <- data.cl                      #WAD統計量を計算するための関数
    mean1 <- rowMeans(as.matrix(x[, cl==1]))#WAD統計量を計算するための関数
    mean2 <- rowMeans(as.matrix(x[, cl==2]))#WAD統計量を計算するための関数
    x_ave <- (mean1 + mean2)/2         #WAD統計量を計算するための関数
    weight <- (x_ave - min(x_ave))/(max(x_ave) - min(x_ave))#WAD統計量を計算するための関数
    statistic <- (mean2 - mean1)*weight#WAD統計量を計算するための関数
    return(statistic)                  #WAD統計量を計算するための関数
}                                                                #WAD統計量を計算するための関数

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
stat_wad <- WAD(data=data, data.cl=data.cl)#WADを実行し、WAD統計量を計算した結果をstat_wadに格納
rank_wad <- rank(-abs(stat_wad))       #WAD統計量の順位を計算した結果をrank_wadに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, stat_wad, rank_wad)#入力データの右側に、「WAD統計量」と「その順位」を結合した結果をtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 6                          #G1群のサンプル数を指定
param_G2 <- 5                          #G2群のサンプル数を指定

#WAD統計量を計算するための関数を定義
WAD <- function(data=NULL, data.cl=NULL){#WAD統計量を計算するための関数
    x <- data                          #WAD統計量を計算するための関数
    cl <- data.cl                      #WAD統計量を計算するための関数
    mean1 <- rowMeans(as.matrix(x[, cl==1]))#WAD統計量を計算するための関数
    mean2 <- rowMeans(as.matrix(x[, cl==2]))#WAD統計量を計算するための関数
    x_ave <- (mean1 + mean2)/2         #WAD統計量を計算するための関数
    weight <- (x_ave - min(x_ave))/(max(x_ave) - min(x_ave))#WAD統計量を計算するための関数
    statistic <- (mean2 - mean1)*weight#WAD統計量を計算するための関数
    return(statistic)                  #WAD統計量を計算するための関数
}                                                                #WAD統計量を計算するための関数

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#発現データのシグナル強度が1未満のものを1にした後にlog2変換
data[data < 1] <- 1                    #1未満のシグナル強度のものを1とする
data <- log(data, 2)                   #log2スケーリング

#本番
stat_wad <- WAD(data=data, data.cl=data.cl)#WADを実行し、WAD統計量を計算した結果をstat_wadに格納
rank_wad <- rank(-abs(stat_wad))       #WAD統計量の順位を計算した結果をrank_wadに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, stat_wad, rank_wad)#入力データの右側に、「WAD統計量」と「その順位」を結合した結果をtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 6                          #G1群のサンプル数を指定
param_G2 <- 5                          #G2群のサンプル数を指定

#必要なパッケージをロード
library(TCC)                           #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
out <- WAD(data, data.cl, logged=F, floor=1)#WADを実行した結果をoutに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, out)#入力データの右側にDEG検出結果を結合したものをtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_Singh_RMA_3274.txt"      #入力ファイル名を指定してin_fに格納
param_G1 <- 50                         #G1群のサンプル数を指定
param_G2 <- 52                         #G2群のサンプル数を指定

#必要なパッケージをロード
library(varSelRF)                      #varSelRFパッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
rf.vs1 <- varSelRF(t(data), factor(data.cl))#RFをデフォルトのパラメータを用いて実行
rf.vs1$selected.vars                   #最終的に選ばれた遺伝子を表示
	

in_f <- "data_Singh_RMA_3274.txt"      #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param_G1 <- 50                         #G1群のサンプル数を指定
param_G2 <- 52                         #G2群のサンプル数を指定

#必要なパッケージをロード
library(st)                            #shrinkage t statisticを計算するためのパッケージの読み込み
library(samr)                          #SAMを計算するためのパッケージの読み込み
library(limma)                         #empirical Bayesを計算するためのパッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
stat_st <- shrinkt.stat(t(data), data.cl)#Shrinkage t統計量を計算し結果をstat_stに格納
rank_st <- rank(-abs(stat_st))         #Shrinkage t統計量の順位を計算し結果をrankt_stに格納
stat_efron <- efront.stat(t(data), data.cl)#Efron's t統計量を計算し結果をstat_efronに格納
rank_efron <- rank(-abs(stat_efron))   #Efron's t統計量の順位を計算し結果をrank_efronに格納
stat_ebayes <- modt.stat(t(data), data.cl)#Empirical Bayes t統計量を計算し結果をstat_ebayesに格納
rank_ebayes <- rank(-abs(stat_ebayes)) #Empirical Bayes t統計量の順位を計算し結果をrank_ebayesに格納
stat_sam <- sam.stat(t(data), data.cl) #SAM's t統計量を計算し結果をstat_samに格納
rank_sam <- rank(-abs(stat_sam))       #SAM's t統計量の順位を計算し結果をrank_samに格納

#ファイルに保存
tmp <- cbind(rownames(data), data, stat_st, rank_st, stat_efron, rank_efron, stat_ebayes, rank_ebayes, stat_sam, rank_sam)#「それぞれの計算方法で得られた統計量」と「その順位」をShrinkage, Efron, empirical Bayes, SAMの順に右のカラムに結合した結果をtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_Singh_RMA_3274.txt"      #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param_G1 <- 50                         #G1群のサンプル数を指定
param_G2 <- 52                         #G2群のサンプル数を指定
param3 <- "pareto"                     #ランキング法を指定

#必要なパッケージなどをロード
#source("http://gap.stat.sinica.edu.tw/Software/mvo.R")#layer ranking algorithmのRスクリプトの読み込み
source("http://www.iu.a.u-tokyo.ac.jp/~kadota/R/mvo.R")#layer ranking algorithmのRスクリプトの読み込み
library(samr)                          #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
data.tmp = list(x=data, y=data.cl, geneid=rownames(data), genenames=rownames(data), logged2=TRUE)#SAMを実行するsamr関数の入力フォーマットに合わせたdata.tmpという入力データを作成
out <- samr(data.tmp, resp.type="Two class unpaired", nperms=20)#samr関数を実行し、結果をoutに格納
stat_sam <- out$tt                     #SAM統計量をstat_samに格納
rank_sam <- rank(-abs(stat_sam))       #SAM統計量の順位をrank_samに格納
stat_fc <- log2(out$foldchange)        #log(FC)統計量をstat_fcに格納
rank_fc <- rank(-abs(stat_fc))         #log(FC)統計量の順位をrank_fcに格納
ranks <- cbind(stat_sam, stat_fc)      #SAM統計量とFC統計量をまとめたものをranksに格納。
ranks.out <- mvo(ranks,ignore=c(T,T),opposite=c(F,F), empty=F, method=param3)#param3で指定した総合ランキングを実行

#ファイルに保存
tmp <- cbind(rownames(data), data, stat_sam, stat_fc, rank_sam, rank_fc, ranks.out)#入力データの右側に「SAM統計量」「log(FC)」「SAM統計量の順位」「log(FC)の順位」「layer ranking」情報を追加してtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_Singh_RMA_3274.txt"      #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param_G1 <- 50                         #G1群のサンプル数を指定
param_G2 <- 52                         #G2群のサンプル数を指定
param3 <- "pareto"                     #ランキング法を指定

#必要な関数などをロード
source("http://www.iu.a.u-tokyo.ac.jp/~kadota/R/mvo.R")#layer ranking algorithmのRスクリプトの読み込み
source("http://www.iu.a.u-tokyo.ac.jp/~kadota/R/R_functions.R")#Welch t-testを行うWelch_ttest関数を含むファイルをあらかじめ読み込む

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
out <- t(apply(data, 1, Welch_ttest, data.cl))#各（行）遺伝子についてt検定を行った結果のt統計量とp値をoutに格納
stat_t <- out[,1]                      #t統計量をstat_tに格納
rank_t <- rank(-abs(stat_t))           #t統計量の順位をrank_tに格納
stat_fc <- apply(data, 1, AD, data.cl) #log(FC)統計量(AD統計量)をstat_fcに格納
rank_fc <- rank(-abs(stat_fc))         #log(FC)統計量の順位をrank_fcに格納
ranks <- cbind(stat_t, stat_fc)        #t統計量とFC統計量をまとめたものをranksに格納。
ranks.out <- mvo(ranks,ignore=c(T,T),opposite=c(F,F), empty=F, method=param3)#param3で指定した総合ランキングを実行

#ファイルに保存
tmp <- cbind(rownames(data), data, stat_t, stat_fc, rank_t, rank_fc, ranks.out)#入力データの右側に「t統計量」「log(FC)」「t統計量の順位」「log(FC)の順位」「layer ranking」情報を追加してtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_Singh_RMA_3274.txt"      #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param_G1 <- 50                         #G1群のサンプル数を指定
param_G2 <- 52                         #G2群のサンプル数を指定
param3 <- "pareto"                     #ランキング法を指定

#必要な関数などをロード
source("http://www.iu.a.u-tokyo.ac.jp/~kadota/R/mvo.R")#layer ranking algorithmのRスクリプトの読み込み
source("http://www.iu.a.u-tokyo.ac.jp/~kadota/R/R_functions.R")#Welch t-testを行うWelch_ttest関数を含むファイルをあらかじめ読み込む

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
out <- t(apply(data, 1, Welch_ttest, data.cl))#各（行）遺伝子についてt検定を行った結果のt統計量とp値をoutに格納
stat_t <- out[,1]                      #t統計量をstat_tに格納
rank_t <- rank(-abs(stat_t))           #t統計量の順位をrank_tに格納
stat_fc <- apply(data, 1, FC, data.cl) #log(FC)統計量をstat_fcに格納
rank_fc <- rank(-abs(stat_fc))         #log(FC)統計量の順位をrank_fcに格納
ranks <- cbind(stat_t, stat_fc)        #t統計量とFC統計量をまとめたものをranksに格納。
ranks.out <- mvo(ranks,ignore=c(T,T),opposite=c(F,F), empty=F, method=param3)#param3で指定した総合ランキングを実行

#ファイルに保存
tmp <- cbind(rownames(data), data, stat_t, stat_fc, rank_t, rank_fc, ranks.out)#入力データの右側に「t統計量」「log(FC)」「t統計量の順位」「log(FC)の順位」「layer ranking」情報を追加してtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "sample14.txt"                 #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param_G1 <- 3                          #G1群のサンプル数を指定
param_G2 <- 3                          #G2群のサンプル数を指定

#必要なパッケージをロード
library(OCplus)                        #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
tmpout <- fdr1d(data, data.cl, verb=FALSE)#入力データ中の発現変動していない遺伝子(non-DEGs)の割合を調べた結果をtmpoutに格納
p0(tmpout)                             #得られたtmpoutの中から目的のnon-DEGsの割合の数値を表示
out <- fdr2d(data, data.cl, p0=p0(tmpout), verb=FALSE)#「標準誤差のlog」と「t統計量」の二つの統計量の値を使って総合FDRを計算

#ファイルに保存
tmp <- cbind(rownames(data), data, out)#入力データの右側に「t統計量(tstat)」「標準誤差のlog(logse)」「総合FDR」情報を追加してtmpに格納。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2_LIV.txt"           #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 4                          #G1群のサンプル数を指定
param_G2 <- 4                          #G2群のサンプル数を指定

#必要なパッケージなどをロード
library(limma)                         #パッケージの読み込み
source("http://eh3.uc.edu/r/ibmtR.R")  #IBMTのRスクリプトの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
#design <- model.matrix(~ as.factor(data.cl))#デザイン行列を作成した結果をdesignに格納
design <- model.matrix(~data.cl)       #デザイン行列を作成した結果をdesignに格納
fit <- lmFit(data, design)             #モデル構築(ばらつきの程度を見積もっている)
fit$Amean<-rowMeans(data)              #おまじない
fit <- IBMT(fit,2)                     #IBMTプログラムの実行
p.value <- fit$IBMT.p                  #p値をp.valueに格納
q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納
ranking <- rank(p.value)               #p.valueでランキングした結果をrankingに格納
sum(q.value < 0.05)                    #FDR < 0.05を満たす遺伝子数を表示

#ファイルに保存
tmp <- cbind(rownames(data), data, p.value, q.value, ranking)#入力データの右側にp.value、q.value、rankingを結合した結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_GSE30533_rma.txt"        #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 5                          #G1群のサンプル数を指定
param_G2 <- 5                          #G2群のサンプル数を指定

#必要なパッケージをロード
library(limma)                         #パッケージの読み込み
source("http://eh3.uc.edu/r/ibmtR.R")  #IBMTのRスクリプトの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
#design <- model.matrix(~ as.factor(data.cl))#デザイン行列を作成した結果をdesignに格納
design <- model.matrix(~data.cl)       #デザイン行列を作成した結果をdesignに格納
fit <- lmFit(data, design)             #モデル構築(ばらつきの程度を見積もっている)
fit$Amean<-rowMeans(data)              #おまじない
fit <- IBMT(fit,2)                     #IBMTプログラムの実行
p.value <- fit$IBMT.p                  #p値をp.valueに格納
q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納
ranking <- rank(p.value)               #p.valueでランキングした結果をrankingに格納
sum(q.value < 0.05)                    #FDR < 0.05を満たす遺伝子数を表示

#ファイルに保存
tmp <- cbind(rownames(data), data, p.value, q.value, ranking)#入力データの右側にp.value、q.value、rankingを結合した結果をtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2_BAT.txt"           #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 4                          #G1群のサンプル数を指定
param_G2 <- 4                          #G2群のサンプル数を指定
param_perm <- 100                      #並べ替え回数を指定(数値が大きいほどより正確だがその分だけ時間がかかる。実際の解析ではサンプル数にもよるが最低でも1000以上を推奨)


#必要なパッケージをロード
library(RankProd)                      #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
out <- RP(data, data.cl, num.perm = param_perm, logged = TRUE, na.rm = FALSE, plot = FALSE, rand = 123)#Rank Product (RP)の実行
stat_RP <- apply(out$RPs, 1, min)      #総合RP統計量を計算(20090527追加)
rank_RP <- as.matrix(rank(stat_RP, ties.method = "min"))#総合順位を計算(20090527追加)
stat_fc <- -out$AveFC                  #総合log(FC)統計量を計算
rank_fc <- rank(-abs(stat_fc))         #総合log(FC)統計量の順位を計算

#ファイルに保存
colnames(out$pfp) <- c("FDR(G1 < G2)","FDR(G1 > G2)")#列名を付与
colnames(out$RPs) <- c("stat(G1 < G2)","stat(G1 > G2)")#列名を付与
colnames(out$RPrank) <- c("rank(G1 < G2)","rank(G1 > G2)")#列名を付与
colnames(stat_fc) <- "log2(G2/G1)"     #列名を付与
tmp <- cbind(rownames(data), data, out$pfp, out$RPs, out$RPrank, stat_fc, rank_fc, stat_RP, rank_RP)#入力データの右側にFDR値, 各統計量, 各順位, 総合log(FC)統計量, その順位, 総合RP統計量, その順位の情報を付加したtmpを用意。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2_BAT.txt"           #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 4                          #G1群のサンプル数を指定
param_G2 <- 4                          #G2群のサンプル数を指定
param_perm <- 100                      #並べ替え回数を指定(数値が大きいほどより正確だがその分だけ時間がかかる。実際の解析ではサンプル数にもよるが最低でも1000以上を推奨)
param_FDR <- 0.05                      #false discovery rate (FDR)閾値を指定

#必要なパッケージをロード
library(RankProd)                      #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
out <- RP(data, data.cl, num.perm = param_perm, logged = TRUE, na.rm = FALSE, plot = FALSE, rand = 123)#Rank Product (RP)の実行。
hoge_upG2 <- rownames(data)[out$pfp[1] < param_FDR]#「G2群 > G1群」の中で指定したFDR閾値を満たすIDを抽出
hoge_upG1 <- rownames(data)[out$pfp[2] < param_FDR]#「G2群 < G1群」の中で指定したFDR閾値を満たすIDを抽出
tmp <- union(hoge_upG1, hoge_upG2)       #（両方で共通して出現しているIDがごくまれにあるので念のため）和集合をとっている

#ファイルに保存
writeLines(tmp, out_f)                 #tmpの中身を指定したファイル名で保存
	

in_f <- "data_Singh_RMA_3274.txt"      #入力ファイル名を指定してin_fに格納
out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 50                         #G1群のサンプル数を指定
param_G2 <- 52                         #G2群のサンプル数を指定
param_perm <- 20                       #並べ替え回数を指定(数値が大きいほどより正確だがその分だけ時間がかかる。実際の解析ではサンプル数にもよるが最低でも1000以上を推奨)

#必要なパッケージをロード
library(RankProd)                      #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
out <- RP(data, data.cl, num.perm = param_perm, logged = TRUE, na.rm = FALSE, plot = FALSE, rand = 123)#Rank Product (RP)の実行。
stat_RP <- apply(out$RPs, 1, min)      #総合RP統計量を計算(20090527追加)
rank_RP <- as.matrix(rank(stat_RP, ties.method = "min"))#総合順位を計算(20090527追加)
stat_fc <- -out$AveFC                  #総合log(FC)統計量を計算
rank_fc <- rank(-abs(stat_fc))         #総合log(FC)統計量の順位を計算

#ファイルに保存
colnames(out$pfp) <- c("FDR(G1 < G2)","FDR(G1 > G2)")#列名を付与
colnames(out$RPs) <- c("stat(G1 < G2)","stat(G1 > G2)")#列名を付与
colnames(out$RPrank) <- c("rank(G1 < G2)","rank(G1 > G2)")#列名を付与
colnames(stat_fc) <- "log2(G2/G1)"     #列名を付与
tmp <- cbind(rownames(data), data, out$pfp, out$RPs, out$RPrank, stat_fc, rank_fc, stat_RP, rank_RP)#入力データの右側にFDR値, 各統計量, 各順位, 総合log(FC)統計量, その順位, 総合RP統計量, その順位の情報を付加したtmpを用意。
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

#以下は（こんなこともできますという）おまけ
plotRP(out, cutoff = 0.05)             #FDRを5%に設定したときのG2群(この場合Tumourサンプル)で高発現（2つの図の上）および低発現（2つの図の下）のdifferentially expressed genesが赤色で示される。この場合だとそれぞれ500個程度あることが分かる。
topGene(out, cutoff = 0.0001, gene.names = rownames(data))#FDR0.01%を満たす遺伝子をリストアップ。（いっぱいあることが分かります。）尚、オプションのgene.namesはrownames(data)とちゃんと指定してあげないと遺伝子名のところがシリアル番号？！(gene.indexカラムの数字と同じ)になってしまいます。
topGene(out, num.gene = 10, gene.names = rownames(data))#上位10遺伝子をリストアップしたいとき
summary(out)                           #outからどのような情報を抽出できるか調べる。
	

in_f <- "data_rma_2_LIV.txt"           #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 4                          #G1群のサンプル数を指定
param_G2 <- 4                          #G2群のサンプル数を指定

#必要なパッケージをロード
library(limma)                         #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
#design <- model.matrix(~ as.factor(data.cl))#デザイン行列を作成した結果をdesignに格納
design <- model.matrix(~data.cl)       #デザイン行列を作成した結果をdesignに格納
fit <- lmFit(data, design)             #モデル構築(ばらつきの程度を見積もっている)
out <- eBayes(fit)                     #検定(経験ベイズ)
p.value <- out$p.value[,ncol(design)]  #p値をp.valueに格納
q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納
ranking <- rank(p.value)               #p.valueでランキングした結果をrankingに格納
sum(q.value < 0.05)                    #FDR < 0.05を満たす遺伝子数を表示

#ファイルに保存
tmp <- cbind(rownames(data), data, p.value, q.value, ranking)#入力データの右側にDEG検出結果を結合したものをtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_GSE30533_rma.txt"        #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 5                          #G1群のサンプル数を指定
param_G2 <- 5                          #G2群のサンプル数を指定

#必要なパッケージをロード
library(limma)                         #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
#design <- model.matrix(~ as.factor(data.cl))#デザイン行列を作成した結果をdesignに格納
design <- model.matrix(~data.cl)       #デザイン行列を作成した結果をdesignに格納
fit <- lmFit(data, design)             #モデル構築(ばらつきの程度を見積もっている)
out <- eBayes(fit)                     #検定(経験ベイズ)
p.value <- out$p.value[,ncol(design)]  #p値をp.valueに格納
q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納
ranking <- rank(p.value)               #p.valueでランキングした結果をrankingに格納
sum(q.value < 0.05)                    #FDR < 0.05を満たす遺伝子数を表示

#ファイルに保存
tmp <- cbind(rownames(data), data, p.value, q.value, ranking)#入力データの右側にDEG検出結果を結合したものをtmpに格納
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存
	

in_f <- "data_rma_2_LIV.txt"           #入力ファイル名を指定してin_fに格納
out_f1 <- "hoge3.txt"                  #出力ファイル名を指定してout_f1に格納
out_f2 <- "hoge3.png"                  #出力ファイル名を指定してout_f2に格納
param_G1 <- 4                          #G1群のサンプル数を指定
param_G2 <- 4                          #G2群のサンプル数を指定
param_FDR <- 0.05                      #false discovery rate (FDR)閾値を指定
param_fig <- c(400, 380)               #ファイル出力時の横幅と縦幅を指定(単位はピクセル)

#必要なパッケージをロード
library(limma)                         #パッケージの読み込み

#入力ファイルの読み込みとラベル情報の作成
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data.cl <- c(rep(1, param_G1), rep(2, param_G2))#G1群を1、G2群を2としたベクトルdata.clを作成

#本番
#design <- model.matrix(~ as.factor(data.cl))#デザイン行列を作成した結果をdesignに格納
design <- model.matrix(~data.cl)       #デザイン行列を作成した結果をdesignに格納
fit <- lmFit(data, design)             #モデル構築(ばらつきの程度を見積もっている)
out <- eBayes(fit)                     #検定(経験ベイズ)
p.value <- out$p.value[,ncol(design)]  #p値をp.valueに格納
q.value <- p.adjust(p.value, method="BH")#q値をq.valueに格納
ranking <- rank(p.value)               #p.valueでランキングした結果をrankingに格納
sum(q.value < 0.05)                    #FDR < 0.05を満たす遺伝子数を表示

#ファイルに保存(テキストファイル)
tmp <- cbind(rownames(data), data, p.value, q.value, ranking)#入力データの右側にDEG検出結果を結合したものをtmpに格納
write.table(tmp, out_f1, sep="\t", append=F, quote=F, row.names=F)#tmpの中身を指定したファイル名で保存

#ファイルに保存(M-A plot)
mean_G1 <- apply(as.matrix(data[,data.cl==1]), 1, mean)#遺伝子ごとにG1群の平均を計算した結果をmean_G1に格納
mean_G2 <- apply(as.matrix(data[,data.cl==2]), 1, mean)#遺伝子ごとにG2群の平均を計算した結果をmean_G2に格納
M <- mean_G2 - mean_G1                 #M-A plotのM値(y軸の値)に相当するものをMに格納
A <- (mean_G1 + mean_G2)/2             #M-A plotのA値(x軸の値)に相当するものをAに格納
png(out_f2, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(A, M, xlab="A = (G2 + G1)/2", ylab="M = G2 - G1", cex=.1, pch=20)#M-A plotを描画
grid(col="gray", lty="dotted")         #指定したパラメータでグリッドを表示
obj <- as.logical(q.value < param_FDR) #条件を満たすかどうかを判定した結果をobjに格納(DEGがTRUE、non-DEGがFALSE)
points(A[obj], M[obj], col="magenta", cex=0.1, pch=20)#objがTRUEとなる要素のみ指定した色で描画
legend("topright", c(paste("DEG(FDR<", param_FDR, ")", sep=""), "non-DEG"),#凡例を作成している
       col=c("magenta", "black"), pch=20)#凡例を作成している
dev.off()                              #おまじない