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

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

param <- "GSE30533"                    #入手したいIDを指定
library(ArrayExpress) 
hoge <- getAE(param, type="raw", extract=F)
	

getwd()
list.files()
	

list.files()
	

list.files(path=getwd())
list.files(path="C:/Users/kadota/Desktop")
list.files(path = "C:/Users/kadota/Desktop")
	

getwd()
list.files(path = "C:/Users/kadota/Desktop/E-GEOD-30533.raw.1")

getwd()
list.files(recursive=TRUE)
	

?list.files
?getAE
	

getwd()
	

out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
library(affy)                          #パッケージの読み込み
hoge <- ReadAffy()                     #*.CELファイルの読み込み
eset <- mas5(hoge)                     #MAS5を実行し、結果をesetに保存
write.exprs(eset, file=out_f)          #結果をout_fで指定したファイル名で保存
	

list.files()
dim(exprs(eset))
	

exprs(eset)[,2]                        #2列目のデータを表示
exprs(eset)[,4:7]                      #4列目から7列目のデータを表示
exprs(eset)[,c(1,6,9)]                 #1,6,9列目のデータを表示
exprs(eset)[1:3,1:2]                   #1,2列目のデータの最初の3行を表示
head(exprs(eset)[,1:2], n=3)           #1,2列目のデータの最初の3行を表示(上と同じ結果)
head(exprs(eset)[,1:2])                #1,2列目のデータの最初の6行を表示(デフォルトは6)
	

colnames(exprs(eset))
colnames(exprs(eset)[,4:7])
colnames(exprs(eset)[,c(1,6,9)])

ccc <- c("GSM757155_-Fe_short_27.CEL", "GSM757160_control_28.CEL")
exprs(eset)[2:4,ccc]                   #cccで指定した列の2-4行目のデータを表示
rpr <- "1367453_at"
exprs(eset)[rpr,ccc]                   #ccc列のrprで指定したprobeset IDデータを表示
param <- "1369078_at"
exprs(eset)[param, 1]                  #1列目のparamで指定したprobeset IDデータを表示

new_colname <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
new_colname
colnames(exprs(eset)) <- new_colname
exprs(eset)[2:4,ccc]                   #cccで指定した列の2-4行目のデータを表示
ccc
exprs(eset)[2:4,]                      #2-4行目のデータを表示
	

library(affy)                          #パッケージの読み込み
hoge <- ReadAffy()                     #*.CELファイルの読み込み

data <- exprs(eset)                    #exprs(eset)をdataに代入
boxplot(data)                          # 1)まずはそのままboxplot
boxplot(data, log="y")                 # 2)y軸をlog10で表示
boxplot(data, log="y", las=3)          # 3)サンプル名を縦書きで表示
boxplot(data, log="y", las=3,          # 4)ラベル情報も追加
  xlab="samples", ylab="signal intensity")# 4)ラベル情報も追加

summary(data[,1:3])
	

summary(data[,1])
mean(data[,1], trim=0.02)
apply(data, 2, mean, trim=0.02)
	

out_f1 <- "hoge2.txt"                  #出力ファイル名を指定してout_f1に格納
out_f2 <- "hoge2.png"                  #出力ファイル名を指定してout_f2に格納
param_fig <- c(500, 300)               #ファイル出力時の横幅と縦幅を指定
###  MAS5前処理法を実行し、ファイルに保存  ###
library(affy)                          #パッケージの読み込み
hoge <- ReadAffy()                     #*.CELファイルの読み込み
eset <- mas5(hoge, normalize=F)        #MAS5を実行し、結果をesetに保存
write.exprs(eset, file=out_f1)         #結果をout_f1で指定したファイル名で保存。
###  箱ひげ図を描画し、ファイルに保存  ###
data <- exprs(eset)                    #exprs(eset)をdataに代入
colnames(data) <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
png(out_f2, width=param_fig[1], height=param_fig[2])
boxplot(data, log="y", las=3,          #boxplot
        xlab="samples", ylab="signal intensity")#boxplot
dev.off()                              #おまじない
	

list.files()
apply(data, 2, mean, trim=0.02)

mean(data[,1], trim=0.02)
NF <- 500/mean(data[,1], trim=0.02)
NF
	

summary(data[,1])
normalized <- data[,1]*NF
summary(normalized)
	

library(affy)                          #パッケージの読み込み
hoge <- ReadAffy()                     #*.CELファイルの読み込み
hoge
sampleNames(hoge)
	

hoge[1]
sampleNames(hoge[1])
probesignal <- hoge[,1]
sampleNames(probesignal)

eset <- mas5(probesignal, normalize=F)
summary(exprs(eset))
	

param <- "1369078_at"
exprs(eset)[param,]

PM <- pm(probesignal, param)
PM
MM <- mm(probesignal, param)
MM
	

source("http://www.bioconductor.org/biocLite.R")#書籍と違っているところ
biocLite("rat2302probe")               #書籍と違っているところ

out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param1 <- "rat2302probe"               #パッケージ名を指定
param2 <- "1369078_at"                 #プローブセットID情報を指定
### 必要なパッケージをロード ###
library(Biostrings)                    #パッケージの読み込み
library(param1, character.only=T)      #param1で指定したパッケージの読み込み
### 前処理(統一的なオブジェクト名に変更しているだけ) ###
hoge <- eval(parse(text=param1))       #オブジェクト名の変更(hogeとして取り扱う)
### 本番(フィルタリング) ###
obj <- is.element(as.character(hoge$Probe.Set.Name), param2)#条件判定
out <- hoge[obj,]                      #objがTRUEとなる行のみ抽出した結果をoutに格納
### 後処理(FASTA形式に変換) ###
fasta <- DNAStringSet(out$sequence)    #DNAStringSetオブジェクトに変換
des <- paste(out$Probe.Set.Name,out$Probe.Interrogation.Position,sep=".")
names(fasta) <- des                    #description行に相当する記述を追加している
### ファイルに保存 ###
writeXStringSet(fasta, file=out_f, format="fasta", width=50)
	

out$sequence

param3 <- 13                           #置換したい塩基の位置を指定
### 関数の作成 ###
s_chartr <- function(str, p) {
  t <- substring(str, p, p)            #置換したい位置の塩基を抽出
  t_c <- chartr("CGAT", "GCTA", t)     #置換後の塩基を作成
  substring(str, p, p) <- t_c          #置換
  return(str)                          #置換後のデータを返す
}
	

MMseq <- s_chartr(out$sequence, param3)#関数を実行
MMseq
	

out_f <- "fig2-3.png"                  #出力ファイル名を指定してout_fに格納
param_fig <- c(500, 300)               #ファイル出力時の横幅と縦幅を指定
param <- "1369078_at"                  #probeset IDを指定
### シグナル強度情報取得 ###
library(affy)
hoge <- ReadAffy()
probesignal <- hoge[,1]
sampleNames(probesignal)
PM <- pm(probesignal, param)
MM <- mm(probesignal, param)
### ファイルに保存 ###
png(out_f, width=param_fig[1], height=param_fig[2])
min <- min(c(PM, MM))                  #y軸の範囲情報（最小値）を取得
max <- max(c(PM, MM))                  #y軸の範囲情報（最大値）を取得
ylim <- c(min, max)                    #y軸の範囲情報をylimに格納
plot(PM, ylim=ylim, ylab="Intensity", pch=20, type="o", col="black")
par(new=T)                             #重ね書き指定
plot(MM, ylim=ylim, ylab="",pch=20, type="o", col="gray")
title(param)                           #タイトルを追加
legend("topright", c("PM", "MM"),col=c("black", "gray"), pch=20)#凡例を追加
dev.off()                              #おまじない
	

(671/87 + 627/297 + 969/234 + 478/197 + 770/81 + 473/97 + 258/194 + 193/286 + 194/83 + 197/58 + 232/92)/11
	

mean(PM/MM)

mean(log2(PM/MM))
mean(log(PM/MM, base=2))
mean(log2(PM) - log2(MM))
	

(4 + 1/4)/2
(log2(4) + log2(1/4))/2

log2(mean(PM/MM))                      # 間違い
mean(log2(PM/MM))                      # 正解

median(log2(PM/MM))                    # 中央値
mean(log2(PM/MM), trim=0.1)            # トリム平均
sort(log2(PM/MM))                      # log比のソート結果を表示
sort(log2(PM/MM))[2:10]                # 自力でトリム後のベクトルを作成
mean(sort(log2(PM/MM))[2:10])          # 別手段
	

x <- log2(PM/MM)
x
library(affy)
tukey.biweight(x)
	

M <- median(x)
S <- mad(x, constant=1)
M
S

abs(x - M)

median(abs(x - M))
	

mean(sort(log2(PM/MM))[2:10])          # 別手段
set.seed(1000)                         # 再現性のある乱数を発生させるためのおまじない
x <- rnorm(10000000)                   # 標準正規分布に従う乱数を発生
sd(x)                                  # SD
mad(x, constant=1)                     # MAD
mad(x)                                 # MAD
mad(x, constant=1.4826)                # MAD
	

out_f <- "fig2-4.png"                  # 出力ファイル名を指定してout_fに格納
param_fig <- c(400, 350)               #ファイル出力時の横幅と縦幅を指定
png(out_f, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
x <- log2(PM/MM)
plot(x, ylim=c(M-3*S, M+3*S), ylab="xj",# xjをプロット
    xlab="j", pch=19, col="black")     # xjをプロット
abline(h=M, col="black", lwd=1)        # Mの線を引く
arrows(8, M, 8, M+1*S, code=3, length=0.13)# j=8のところで上向きの矢印を描く
arrows(8, M, 8, M-1*S, code=3, length=0.13)# j=8のところで上向きの矢印を描く
abline(h=M+1*S, col="black", lwd=1, lty=2)# M+1*MAD
abline(h=M-1*S, col="black", lwd=1, lty=2)# M-1*MAD
legend("topright", "M", col="black", lty=1)# 図の凡例を追加
dev.off()                              # おまじない
	

x <- log2(PM/MM)                       # 書籍と違っているところ(上のおまけを実行しなかったらxオブジェクトには標準正規分布乱数が入ったままなので...)

c <- 5
e <- 0.0001
u <- (x - M)/(c*S + e)
u
	

w <- 1 - abs(u)
matome <- cbind(x, u, w)
colnames(matome) <- c("x", "u", "w")
matome
	

out_f <- "fig2-5.png"                  # 出力ファイル名を指定してout_fに格納
param_fig <- c(400, 350)               #ファイル出力時の横幅と縦幅を指定
png(out_f, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
plot(u, ylim=c(-1, 1), ylab="uj",      # ujをプロット
    xlab="j", pch=19, col="black")     # ujをプロット
for(i in 1:10){
    text(x=i, y=u[i], sprintf("%3.2f", w[i]), pos=4, cex=1.2)
}
abline(h=0, col="black", lwd=1)        # u = 0の線を引く
legend("topright", "M", col="black", lty=1)# 図の凡例を追加
dev.off()                              # おまじない
	

sum(w*x)/sum(w)

pnorm(1.96)*2 - 1
pnorm(2.58)*2 - 1
pnorm(3.37)*2 - 1

c <- 2
e <- 0.0001
u <- (x - M)/(c*S + e)
w <- 1 - abs(u)
matome <- cbind(x, u, w)
colnames(matome) <- c("x", "u", "w")
matome
	

u[abs(u) > 1] <- 1
w <- 1 - abs(u)
matome <- cbind(x, u, w)
colnames(matome) <- c("x", "u", "w")
matome
sum(w*x)/sum(w)
	

c <- 5
e <- 0.0001
u <- (x - M)/(c*S + e)
u[abs(u) > 1] <- 1
w <- (1 - (u)^2)^2
matome <- cbind(x, u, w)
colnames(matome) <- c("x", "u", "w")
matome
sum(w*x)/sum(w)
tukey.biweight(x, c=5, epsilon=0.0001)

w <- (1 - (abs(u))^3)^3
matome <- cbind(x, u, w)
colnames(matome) <- c("x", "u", "w")
matome
	

sum(w*x)/sum(w)

M <- tukey.biweight(x, c=5, epsilon=0.0001)
S <- mad(x, center=M, constant=1)
u <- (x - M)/(c*S + e)
u[abs(u) > 1] <- 1
w <- (1 - (u)^2)^2
sum(w*x)/sum(w)

library(RobLoxBioC)
estimate(roblox(x, eps.lower=0.0, eps.upper=0.05, k=1))[1]
estimate(roblox(x, eps.lower=0.0, eps.upper=0.05, k=3))[1]
	

estimate(roblox(x, eps.lower=0.0, eps.upper=0.05, k=10))[1]

SB <- tukey.biweight(log2(PM/MM))
IM <- MM
obj <- ((MM >= PM) & SB > 0.03)
IM[obj] <- 2^(-SB)*PM[obj]
matome <- cbind(PM, MM, IM)
colnames(matome) <- c("PM", "MM", "IM")
matome
	

PV <- log2(PM - IM)
SigLogVal <- tukey.biweight(PV)
2^SigLogVal
	

out_f <- "hoge5.txt"                   #出力ファイル名を指定してout_fに格納
library(affy)                          #パッケージの読み込み
hoge <- ReadAffy()                     #*.CELファイルの読み込み
eset <- rma(hoge)                      #RMAを実行し、結果をesetに保存
write.exprs(eset, file=out_f)          #結果をout_fで指定したファイル名で保存
	

dim(exprs(eset))
sampleNames(hoge)

param <- 1:4
sampleNames(hoge[,param])
eset2 <- rma(hoge[,param])
dim(exprs(eset2))

head(exprs(eset)[1:3, 1:2])
head(exprs(eset2)[1:3, 1:2])
	

param <- 1:5
eset3 <- rma(hoge[,param])
param <- 5:9
sampleNames(hoge[,param])
eset4 <- rma(hoge[,param])
head(exprs(eset3)[1:3,4:5])
head(exprs(eset4)[1:3,1:2])

cor(exprs(eset3)[,5], exprs(eset4)[,1], method="spearman")
	

library(affy)
hoge <- ReadAffy()
newname <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
sampleNames(hoge) <- newname
hist(hoge)                             # 図2-6作成部分

normalized <- normalize(hoge)
hist(normalized)                       # 図2-7作成部分

dim(pm(hoge))                          # PMシグナル強度の行数と列数を表示
	

head(pm(hoge), n=3)                    # PMシグナル強度の最初の3行分を表示
hoge_sorted <- apply(pm(hoge),2,sort)  # 列ごとにソート
head(hoge_sorted, n=7)                 # ソート後の行列の最初の7行分を表示
reference <- rowMeans(hoge_sorted)     # 行ごとの平均値を算出
head(reference)
hist(log2(reference), probability=T)   # 図2-8作成部分
	

summary(pm(normalized)[,1:3])

normalized_sorted <- apply(pm(normalized),2,sort)# 列ごとにソート
head(normalized_sorted, n=7)           # ソート後の行列の最初の7行分を表示
	

param <- "1369078_at"                  # プローブセットIDを指定
PM <- pm(hoge, param)                  # PMシグナル情報取得
PM

apply(PM, 1, sd)                       # アレイ間のばらつき
mean(apply(PM, 1, sd))                 # アレイ間のばらつきの平均

apply(PM, 2, sd)                       # プローブ間のばらつき
mean(apply(PM, 2, sd))                 # プローブ間のばらつきの平均
	

library(frma)
eset <- frma(hoge)
	

param1 <- "GSE30533"                   # GEO IDを指定
param2 <- "hoge_"                      # 出力ファイル名の最初の部分を指定
library(GEOquery)                      # パッケージの読み込み
data <- getGEO(param1)                 # 指定したGEO IDのデータを取得
hoge <- sapply(data, annotation)       # 用いられたアレイ情報(GPL ID)をhogeに格納
for(i in 1:length(hoge)){              # hogeの要素数(用いたアレイ数)分ループを回す
    out_f <- paste(param2, hoge[i], ".txt", sep="")# ファイル名を作成
    out <- data[[i]]@featureData@data  # アノテーション情報抽出
    write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F) 
}
	

out_f <- "hoge6.txt"                   # 出力ファイル名を指定してout_fに格納
param <- "GPL1355"                     # 入手したいGEO IDを指定
library(GEOquery)                      # パッケージの読み込み
data <- getGEO(param)                  # 指定したGEO IDのデータを取得
out <- data@dataTable@table            # アノテーション情報抽出結果をoutに格納
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)
	

in_f <- "hoge1.txt"
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
colnames(data) <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
data.dist <- as.dist(1 - cor(data, method = "spearman"))
out <- hclust(data.dist, method = "average")
plot(out)                              # 図3-1作成部分
	

head(data[,c(8, 10)])
	

cor(data[,8], data[,10], method="spearman")
1 - cor(data[,8], data[,10], method="spearman")

cor(data[,8], data[,10], method="pearson")
cor(rank(data[,8]), rank(data[,10]), method="pearson")
	

in_f <- "hoge1.txt"
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
colnames(data) <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
par(mfrow=c(1,2))
### 最短距離法(single-linkage clustering) ###
data.dist <- as.dist(1 - cor(data, method="spearman"))
out <- hclust(data.dist, method = "single")
plot(out)                              # 図3-2左作成部分
### 最長距離法(complete-linkage clustering) ###
data.dist <- as.dist(1 - cor(data, method="spearman"))
out <- hclust(data.dist, method = "complete")
plot(out)                              # 図3-2右作成部分
	

summary(data[,8])
	

head(data[,8])
mean(data[,8])
sum(data[,8] < mean(data[,8]))
25486/31099
qqnorm(data[,8])
	

par(mfrow=c(1,2))
### 対数変換前の"G2_3"データのQ-Qプロット ###
qqnorm(data[,8])                       # 図3-3左作成部分
### 正規分布に従う乱数を10,000個発生させて描いたQ-Qプロット ###
set.seed(1000)
hoge <- rnorm(10000)
qqnorm(hoge)                           # 図3-3右作成部分
	

ks.test(hoge, "pnorm", mean=mean(hoge), sd=sd(hoge))
	

shapiro.test(hoge)
shapiro.test(rnorm(5000))

hoge <- log2(data[,8])
ks.test(hoge, "pnorm", mean=mean(hoge), sd=sd(hoge))
	

par(mfrow=c(1,2))
qqnorm(hoge)                           # 図3-4左作成部分
hist(hoge)                             # 図3-4右作成部分
	

in_f <- "hoge1.txt"
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
colnames(data) <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
par(mfrow=c(1,2))
### 対数変換前のデータ ###
data.dist <- as.dist(1 - cor(data, method = "pearson"))
out <- hclust(data.dist, method = "average")
plot(out)                              # 図3-5左作成部分
### 対数変換後のデータ ###
data.dist <- as.dist(1 - cor(log2(data), method = "pearson"))
out <- hclust(data.dist, method = "average")
plot(out)                              # 図3-5右作成部分
	

in_f <- "hoge1.txt"
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
colnames(data) <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
par(mfrow=c(1,2))
### マンハッタン距離 ###
data.dist <- dist(t(log2(data)), method = "manhattan")
out <- hclust(data.dist, method = "average")
plot(out)                              # 図3-6左作成部分
### ユークリッド距離 ###
data.dist <- dist(t(log2(data)), method = "euclidean")
out <- hclust(data.dist, method = "average")
plot(out)                              # 図3-6右作成部分
	

out_f1 <- "GSE30533_rma.txt"
out_f2 <- "GSE30533_rob.txt"
library(affy)
hoge <- ReadAffy()
### RMA正規化の実行および(対数変換後のデータなのでそのまま)ファイル保存 ###
eset <- rma(hoge)
write.exprs(eset, file=out_f1)
data_rma <- exprs(eset)
colnames(data_rma) <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
### RobLoxBioC正規化の実行、対数変換、およびファイル保存 ###
library(RobLoxBioC)
eset <- robloxbioc(hoge)
exprs(eset) <- log2(exprs(eset))
write.exprs(eset, file=out_f2)
data_rob <- exprs(eset)
colnames(data_rob) <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
### サンプル間クラスタリング(左：RMAデータ、右：RobLoxBioCデータ) ###
par(mfrow=c(1,2))
data.dist <- as.dist(1 - cor(data_rma, method = "spearman"))
out <- hclust(data.dist, method = "average")
plot(out)                              # 図3-7左作成部分
data.dist <- as.dist(1 - cor(data_rob, method = "spearman"))
out <- hclust(data.dist, method = "average")
plot(out)                              # 図3-7右作成部分
	

param_G1 <- 4                          #G1群のサンプル数を指定
param_G2 <- 4                          #G2群のサンプル数を指定
param_Ngene <- 10000                   #遺伝子数を指定
### ランダムデータの生成 ###
set.seed(1000)
hoge <- rnorm(param_Ngene*(param_G1+param_G2))
data <- matrix(hoge, nrow=param_Ngene, byrow=T)
rownames(data) <- paste("gene", 1:param_Ngene, sep="_")
colnames(data) <- c(paste("G1_", 1:param_G1, sep=""),
                    paste("G2_", 1:param_G2, sep=""))
	

dim(data)
summary(data)
apply(data, 2, sd)
	

param_range <- c(-4.7, 4.7)
param_label <- c("G1_1", "G1_2")
### 散布図の作成 ###
plot(data[,1], data[,2], pch=20, cex=0.1,# 図3-8作成部分
     xlim=param_range, ylim=param_range,# 図3-8作成部分
     xlab=param_label[1], ylab=param_label[2])# 図3-8作成部分
grid(col="gray", lty="dotted")         # 図3-8作成部分
	

cor.test(data[,1],data[,2])
	

par(mfrow=c(1,2))
data.dist <- as.dist(1 - cor(data, method = "pearson"))
out <- hclust(data.dist, method = "average")
plot(out)                              # 図3-9左作成部分
data.dist <- as.dist(1 - cor(data, method = "spearman"))
out <- hclust(data.dist, method = "average")
plot(out)                              # 図3-9右作成部分
	

param_G1 <- 4                          #G1群のサンプル数を指定
param_G2 <- 4                          #G2群のサンプル数を指定
library(genefilter)
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
### 本番(t検定) ###
out <- rowttests(data, as.factor(data.cl))
	

head(out, n=5)
	

data[1,]
rank(data[1,])

t.test(data[1,data.cl==1], data[1,data.cl==2], var.equal=T)
	

data[1,data.cl==1]
mean(data[1,data.cl==1])
mean(data[1,data.cl==2])
mean(data[1,data.cl==1]) - mean(data[1,data.cl==2])
?rowttests

p.value <- out$p.value
head(p.value)
sum(p.value < 0.05)
sum(p.value < 0.01)
	

min(p.value)
obj <- p.value == min(p.value)
out[obj,]
data[obj,]
hist(p.value)                          # 図3-10作成部分

nrow(data)
	

threshold <- 0.05
Bonferroni_thres <- threshold/nrow(data)
Bonferroni_thres
sum(p.value < Bonferroni_thres)
	

threshold <- c(0.001, 0.01, 0.03, 0.05, 0.10)
res <- NULL
for(i in 1:length(threshold)){
  observed <- sum(p.value < threshold[i])
  expected <- nrow(data)*threshold[i]
  FDR <- expected/observed
  res <- rbind(res, c(threshold[i], observed, expected, FDR))
}
colnames(res) <- c("threshold", "observed", "expected", "FDR")
	

res
head(data, n=3)
data[1:1000,1:4] <- data[1:1000, 1:4] + 3
head(data, n=3)
	

out <- rowttests(data, as.factor(data.cl))
p.value <- out$p.value
hist(p.value)                          # 図3-11作成部分

threshold <- 0.05
Bonferroni_thres <- threshold/nrow(data)
sum(p.value < Bonferroni_thres)
	

threshold <- c(0.001, 0.01, 0.03, 0.05, 0.10, 0.15)
res <- NULL
for(i in 1:length(threshold)){
  observed <- sum(p.value < threshold[i])
  expected <- nrow(data)*threshold[i]
  FDR <- expected/observed
  res <- rbind(res, c(threshold[i], observed, expected, FDR))
}
colnames(res) <- c("threshold", "observed", "expected", "FDR")
	

res
	

res[,2] - res[,3]

q.value <- p.adjust(p.value, method="BH")
sum(q.value < 0.1)
sum(q.value < 0.04405286)
sum(q.value < 0.12970169)
	

par(mfrow=c(1,2))
data.dist <- as.dist(1 - cor(data, method = "pearson"))
hoge <- hclust(data.dist, method = "average")
plot(hoge)                             # 図3-12左作成部分
data.dist <- as.dist(1 - cor(data, method = "spearman"))
hoge <- hclust(data.dist, method = "average")
plot(hoge)                             # 図3-12右作成部分
	

sum(abs(out$dm) > 1)

data.cl
colnames(data)
set.seed(1)
posi <- sample(8)
posi
	

data.cl[posi]
out <- rowttests(data, as.factor(data.cl[posi]))
sum(abs(out$dm) > 1)

posi <- sample(8)
posi
out <- rowttests(data, as.factor(data.cl[posi]))
sum(abs(out$dm) > 1)
posi <- sample(8)
posi
out <- rowttests(data, as.factor(data.cl[posi]))
sum(abs(out$dm) > 1)
	

in_f <- "hoge1.txt"
param_posi <- c(2, 4)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
colnames(data) <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
### 散布図の作成 ###
par(mfrow=c(1,2))
unlogged <- data[,param_posi]
plot(unlogged, pch=20, cex=0.1)        # 図3-13左作成部分
grid(col="gray", lty="dotted")         # 図3-13左作成部分
logged <- log2(data[,param_posi])
plot(logged, pch=20, cex=0.1)          # 図3-13右作成部分
grid(col="gray", lty="dotted")         # 図3-13右作成部分
	

head(logged)
M <- logged[,2] - logged[,1]
A <- (logged[,2] + logged[,1])/2
plot(A, M)                             # 図3-14作成部分
grid(col="gray", lty="dotted")         # 図3-14作成部分
	

obj <- (abs(M) >= 1)
points(A[obj], M[obj], col="gray")
sum(obj)
sum(obj)/length(obj)
	

points(13, 4, col="black", cex=2, pch=19)
points(3, 4, col="black", cex=2, pch=15)
	

obj <- (A < 5)
sum(obj)
summary(M[obj])

head(sort(abs(M[obj]), decreasing=T))
sort(abs(M[obj]), decreasing=T)[82]
plot(A, M)                             # 図3-17作成部分
grid(col="gray", lty="dotted")         # 図3-17作成部分(グリッド線を追加)
points(A[obj], M[obj], col="gray")     # 図3-17作成部分(A < 5に相当する部分を灰色にしている)
obj2 <- abs(M) > 4.536463              # 図3-17作成部分(Mの絶対値が4.536463より大きい遺伝子の位置情報を取得)
hoge <- as.logical(obj*obj2)           # 図3-17作成部分(|M| > 4.536463かつA < 5を満たすものの位置情報を取得)
points(A[hoge], M[hoge], col="black", pch=20)# 図3-17作成部分(黒丸をつけている)
	

obj <- (A > 10)
sum(obj)
sort(abs(M[obj]), decreasing=T)[45]
	

in_f <- "GSE30533_rob.txt"
param_posi <- c(2, 4)
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
### 散布図の作成 ###
data <- data[,param_posi]
M <- data[,2] - data[,1]
A <- (data[,2] + data[,1])/2
plot(A, M)                             # 図3-18作成部分
grid(col="gray", lty="dotted")         # 図3-18作成部分(グリッド線を追加)
	

in_f <- "hoge1.txt"
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
par(mfrow=c(1,2))
### G1群の平均-分散プロット(対数変換前) ###
hoge <- data[,1:5]
MEAN <- apply(hoge, 1, mean)
VARIANCE <- apply(hoge, 1, var)
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1,# 図3-19左作成部分
     xlim=c(1e-02,1e+08), ylim=c(1e-02,1e+08), col="black")# 図3-19左作成部分
grid(col="gray", lty="dotted")         # 図3-19左作成部分(グリッド線を追加)
abline(a=0, b=1, col="gray")           # 図3-19左作成部分(y=xの灰色直線を追加)
### G1群の平均-分散プロット(対数変換後) ###
hoge <- log2(data[,1:5])
MEAN <- apply(hoge, 1, mean)
VARIANCE <- apply(hoge, 1, var)
plot(MEAN, VARIANCE, pch=20, cex=.1, xlim=c(-2,16), ylim=c(-2,16), col="black")# 図3-19右作成部分
grid(col="gray", lty="dotted")         # 図3-19右作成部分(グリッド線を追加)
abline(a=0, b=1, col="gray")           # 図3-19右作成部分(y=xの灰色直線を追加)
	

in_f1 <- "GSE30533_rma.txt"
in_f2 <- "GSE30533_rob.txt"
par(mfrow=c(1,2))
### G1群の平均-分散プロット(RMAデータ) ###
data <- read.table(in_f1, header=TRUE, row.names=1, sep="\t", quote="")
hoge <- 2^data[,1:5]
MEAN <- apply(hoge, 1, mean)
VARIANCE <- apply(hoge, 1, var)
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1, xlim=c(1e-02,1e+07), ylim=c(1e-02,1e+07), col="black")
grid(col="gray", lty="dotted")         # 図3-20左作成部分(グリッド線を追加)
abline(a=0, b=1, col="gray")           # 図3-20左作成部分(y=xの灰色直線を追加)
### G1群の平均-分散プロット(RobLoxBioCデータ) ###
data <- read.table(in_f2, header=TRUE, row.names=1, sep="\t", quote="")
hoge <- 2^data[,1:5]
MEAN <- apply(hoge, 1, mean)
VARIANCE <- apply(hoge, 1, var)
plot(MEAN, VARIANCE, log="xy", pch=20, cex=.1, xlim=c(1e-02,1e+07), ylim=c(1e-02,1e+07), col="black")
grid(col="gray", lty="dotted")         # 図3-20右作成部分(グリッド線を追加)
abline(a=0, b=1, col="gray")           # 図3-20右作成部分(y=xの灰色直線を追加)
	

in_f <- "hoge1.txt"                    #入力ファイル名を指定してin_fに格納
out_f <- "GSE30533_MAS_DEG.txt"        #出力ファイル名を指定してout_fに格納
param_G1 <- 5                          #G1群のサンプル数を指定
param_G2 <- 5                          #G2群のサンプル数を指定
### ファイルの読み込み、log2変換、列名変更 ###
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data <- log2(data)
colnames(data) <- c(paste("G1_", 1:5, sep=""), paste("G2_", 1:5, sep=""))
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
### DEG検出 ###
library(limma)
design <- model.matrix(~ as.factor(data.cl))
fit <- lmFit(data, design)
out <- eBayes(fit)
p.value <- out$p.value[,ncol(design)]
q.value <- p.adjust(p.value, method="BH")
ranking <- rank(p.value)
### ファイルに保存 ###
tmp <- cbind(rownames(data), data, p.value, q.value, ranking)
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)
	

topTable(out, coef=colnames(design)[ncol(design)], adjust="BH", number=8)

head(sort(p.value))
head(sort(q.value))
sum(q.value < 0.1)
sum(q.value < 0.05)
	

length(q.value)
	

mean_G1 <- apply(as.matrix(data[,data.cl==1]), 1, mean)
mean_G2 <- apply(as.matrix(data[,data.cl==2]), 1, mean)
M <- mean_G2 - mean_G1
A <- (mean_G1 + mean_G2)/2
plot(A, M)                             # 図4-1作成部分
grid(col="gray", lty="dotted")         # 図4-1作成部分(グリッド線を追加)
obj <- q.value < 0.1
points(A[obj], M[obj], col="gray", pch=15, cex=1.5)
	

M[obj]
A[obj]
	

in_f <- "GSE30533_rob.txt"             #入力ファイル名を指定してin_fに格納
param_G1 <- 5                          #G1群のサンプル数を指定
param_G2 <- 5                          #G2群のサンプル数を指定
param_FDR <- 0.1                       #DEG検出時のFDR閾値を指定
### ファイルの読み込み ###
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data.cl <- c(rep(1, param_G1), rep(2, param_G2))
### DEG検出 ###
library(limma)
design <- model.matrix(~ as.factor(data.cl))
fit <- lmFit(data, design)
out <- eBayes(fit)
hoge <- topTable(out,coef=colnames(design)[ncol(design)], 
                 adjust="BH",number= nrow(data))
### M-A plot ###
plot(hoge$AveExpr, hoge$logFC, xlab="AveExpr", ylab="logFC")
grid(col="gray", lty="dotted")
obj <- hoge$adj.P.Val < param_FDR
points(hoge$AveExpr[obj], hoge$logFC[obj], col="gray", pch=15, cex=1.5)
	

sum(obj)
head(hoge)
	

AD <- hoge$logFC 
A <- hoge$AveExpr
w <- (A - min(A))/(max(A) - min(A))
WAD <- AD * w
names(WAD) <- rownames(hoge)
ranking <- rank(-abs(WAD))
	

head(WAD)
head(ranking)
w[4]
w[2]
AD[4]
AD[2]
AD[4]*w[4]
AD[2]*w[2]
	

### M-A plot ###
plot(A, AD, xlab="AveExpr", ylab="logFC")# 図4-2作成部分
grid(col="gray", lty="dotted")         # 図4-2作成部分(グリッド線を追加)
### limmaランキング結果上位10個を塗りつぶし灰色四角で表示 ###
obj <- rank(hoge$P.Value) <= 10
points(A[obj], AD[obj], col="gray", pch=15, cex=1.5)
### WADランキング結果上位10個を塗りつぶし黒丸で表示 ###
obj <- ranking <= 10
points(A[obj], AD[obj], col="black", pch=16, cex=1.2)
### Legend ###
legend("bottomright", c("limma","WAD"), col=c("gray","black"), pch=c(15,16))
	

G1 <- c(10.1, 7.5, 12.5, 9.3, 6.1)
G2 <- c(10.8, 8.1, 13.0, 9.8, 6.8)
t.test(G1, G2, paired=TRUE)
t.test(G1, G2, paired=FALSE)
	

in_f <- "GSE30533_rob.txt"             #入力ファイル名を指定してin_fに格納
out_f <- "GSE30533_rob_pair_DEG.txt"   #出力ファイル名を指定してout_fに格納
param <- 5                             #個体数またはsib-ships数を指定
### ファイルの読み込み ###
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
pair <- as.factor(rep(1:param, 2))
cl <- as.factor(c(rep(1, param), rep(2, param)))
### DEG検出 ###
library(limma)
design <- model.matrix(~ pair + cl)
fit <- lmFit(data, design)
out <- eBayes(fit)
p.value <- out$p.value[,ncol(design)]
q.value <- p.adjust(p.value, method="BH")
ranking <- rank(p.value)
### ファイルに保存 ###
tmp <- cbind(rownames(data), data, p.value, q.value, ranking)
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)
	

topTable(out, coef=colnames(design)[ncol(design)], adjust="BH", number=4)
	

pair
cl
design
dim(design)
colnames(design)
colnames(design)[ncol(design)]
head(out$p.value, n=3)
	

in_f <- "GSE30533_rob.txt"             #入力ファイル名を指定してin_fに格納
### ファイルの読み込み ###
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
factor1 <- c(1,1,1,1,1,2,2,2,2,2)
factor2 <- c(1,1,2,2,2,1,1,1,2,2)
merge <- as.factor(paste(factor1, factor2, sep="."))
### DEG検出 ###
library(limma)
design <- model.matrix(~ 0 + merge)
fit <- lmFit(data, design)
	

merge
design
	

contrast <- makeContrasts(
    analysis1 = merge1.2 - merge2.2,
    analysis2 = merge1.1 - merge2.1,
    analysis3 = (merge1.2 - merge2.2) - (merge1.1 - merge2.1),
    levels = design)
fit2 <- contrasts.fit(fit, contrast)
out <- eBayes(fit2)
	

head(out$p.value, n=3)
	

topTable(out, coef="analysis1", adjust="BH", number=3)
topTable(out, coef=1, adjust="BH", number=3)
topTable(out, coef=2, adjust="BH", number=3)
topTable(out, coef=3, adjust="BH", number=3)

vennDiagram(decideTests(out, adjust.method="BH", p.value=0.7)) 
	

in_f <- "GSE30533_rob.txt"             #入力ファイル名を指定してin_fに格納
out_f <- "GSE30533_rob_multif_DEG.txt" #出力ファイル名を指定してout_fに格納
### ファイルの読み込み ###
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
factor1 <- c(rep("G1", 5), rep("G2", 5))
factor2 <- c("F","F","M","M","M", "F","F","F","M","M")
merge <- as.factor(paste(factor1, factor2, sep="."))
### DEG検出 ###
library(limma)
design <- model.matrix(~ 0 + merge)
colnames(design) <- levels(merge)
fit <- lmFit(data, design)
contrast <- makeContrasts(
    G1vsG2inM = G1.M - G2.M,
    G1vsG2inF = G1.F - G2.F,
    Diff = (G1.M - G2.M) - (G1.F - G2.F),
    levels = design)
fit2 <- contrasts.fit(fit, contrast)
out <- eBayes(fit2)
### ファイルに保存 ###
tmp <- cbind(rownames(data), data, out$p.value)
	

topTable(out, coef="G1vsG2inM", adjust="BH", number=3)
topTable(out, coef=1, adjust="BH", number=3)
topTable(out, coef=2, adjust="BH", number=3)
topTable(out, coef=3, adjust="BH", number=3)

vennDiagram(decideTests(out, adjust.method="BH", p.value=0.7)) 
	

in_f <- "GSE30533_rob.txt"             #入力ファイル名を指定してin_fに格納
out_f <- "GSE30533_rob_3groups_DEG.txt"#出力ファイル名を指定してout_fに格納
param_G1 <- 4                          #G1群のサンプル数を指定
param_G2 <- 3                          #G2群のサンプル数を指定
param_G3 <- 3                          #G3群のサンプル数を指定
### ファイルの読み込み ###
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")
data.cl <- c(rep("G1", param_G1), rep("G2", param_G2) , rep("G3", param_G3))
### DEG検出 ###
library(limma)
design <- model.matrix(~ 0 + as.factor(data.cl))
colnames(design) <- levels(as.factor(data.cl))
fit <- lmFit(data, design)
contrast <- makeContrasts(
    G1vsG2 = G1 - G2, G1vsG3 = G1 - G3, G2vsG3 = G2 - G3, levels = design)
fit2 <- contrasts.fit(fit, contrast)
out <- eBayes(fit2)
### ファイルに保存 ###
tmp <- cbind(rownames(data), data, out$p.value)
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)
	

dim(out$p.value)
head(out$p.value, n=4)
param_FDR <- 0.5
sum(p.adjust(out$p.value[,1], method="BH") < param_FDR)
sum(p.adjust(out$p.value[,2], method="BH") < param_FDR)
sum(p.adjust(out$p.value[,3], method="BH") < param_FDR)

vennDiagram(decideTests(out, adjust.method="BH", p.value=0.5))

q.value <- apply(out$p.value, MARGIN=2, p.adjust, method="BH")
head(q.value, n=4)
	

sum(q.value[,3] < 0.5)
	

library(TCC)
data(hypoData_ts)
dim(hypoData_ts)
?hypoData_ts

hypoData_ts[1,]
	

data <- hypoData_ts
par(mfrow=c(2,4))                      #2行×4列からなる図を作成し、上の行(row)から順番に埋めていくという宣言
for(i in 1:nrow(data)){
    plot(data[i,], ylab="Expression", xlab=rownames(data)[i], ylim=c(0,10))
}
	

out_f <- "fig4-5.png"                  #出力ファイル名を指定してout_fに格納
param_fig <- c(600, 350)               #ファイル出力時の横幅と縦幅を指定
library(TCC)
data(hypoData_ts)
data <- hypoData_ts
### ファイルに保存 ###
png(out_f, width=param_fig[1], height=param_fig[2])
par(mfrow=c(2,4))
for(i in 1:nrow(data)){
    plot(data[i,], ylab="Expression", xlab=rownames(data)[i], ylim=c(0,10),
        cex=2.0, cex.lab=1.7, cex.axis=1.5, pch=1)
}
dev.off()
	

out_f <- "result_ROKU.txt"             #出力ファイル名を指定してout_fに格納
### ファイルの読み込み ###
library(TCC)
data(hypoData_ts)
data <- hypoData_ts
### ROKUの実行 ###
out <- ROKU(data, upper.limit=0.25)
### ファイルに保存 ###
tmp <- cbind(rownames(data), out$outliers, out$modH, out$rank)
write.table(tmp, out_f, sep="\t", append=F, quote=F, row.names=F)
	

dim(out$outlier)
out$outlier
	

out$H
sort(out$H)
	

out$Tbw

x <- data[6,]
Tbw <- out$Tbw[6]
y <- abs(x - Tbw)
x
Tbw
y
	

out_f <- "fig4-6.png"                  #出力ファイル名を指定してout_fに格納
param_fig <- c(600, 350)               #ファイル出力時の横幅と縦幅を指定
library(TCC)
data(hypoData_ts)
data <- hypoData_ts
### ROKUの実行とデータ変換 ###
out <- ROKU(data, upper.limit=0.25)
y <- abs(data - out$Tbw)
### ファイルに保存 ###
png(out_f, width=param_fig[1], height=param_fig[2])
par(mfrow=c(2,4))
for(i in 1:nrow(data)){
    ### データ変換前 ###
    plot(data[i,], ylab="Expression", xlab=rownames(data)[i], ylim=c(0,10),#黒丸をプロット
        cex=2.0, cex.lab=1.7, cex.axis=1.5, col="black", pch=1)#黒丸をプロット(pch=1は塗りつぶさない丸に相当、cex=2.0は丸の大きさを標準の2.0倍に、cex.lab=1.7は"Expression"に相当する軸のラベルの大きさを標準の1.7倍に、など)
    ### Tukey biweight (Tbw)の水平線 ###
    abline(h=out$Tbw[i], lty="dashed", col="black")
    ### データ変換後 ###
    points(y[i,], cex=2.0, cex.lab=1.7, cex.axis=1.5, col="gray", pch=17)#pch=17が塗りつぶし三角に相当
}
dev.off()
	

out$modH
sort(out$modH)
	

param <- "GSE30533"                    #入手したいIDを指定

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

#本番(データ取得)
hoge <- getPlatforms(param)            #アレイ（プラットフォーム）情報を取得
eset <- getDataset(param, hoge[[1]], norm="ORIGINAL", genes=F)#paramで指定したIDのsubmitterが登録した正規化後のデータ(norm="ORIGINAL")を取得
dim(exprs(eset))                       #行数と列数を表示
head(exprs(eset))                      #最初の数行を表示
	

param <- "GSE30533"                    #入手したいIDを指定

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

#本番(データ取得)
hoge <- getPlatforms(param)            #アレイ（プラットフォーム）情報を取得
eset <- getDataset(param, hoge[[1]], norm="FRMA", genes=F)#paramで指定したIDのfRMA正規化後のデータを取得
dim(exprs(eset))                       #行数と列数を表示
	

param <- "GSE30533"                    #入手したいIDを指定

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

#本番(データ取得)
hoge <- getPlatforms(param)            #アレイ（プラットフォーム）情報を取得
eset <- getDataset(param, hoge[[1]], norm="ORIGINAL", genes=T)#paramで指定したIDのsubmitterが登録した正規化後のデータ(log2 RMA)を取得
dim(exprs(eset))                       #行数と列数を表示
head(exprs(eset))                      #最初の数行を表示
	

out_f <- "data_rma.txt"                #出力ファイル名を指定してout_fに格納
param <- "E-MEXP-1422"                 #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのデータを取得した結果をhogeに格納

#本番
eset <- rma(hoge)                      #RMAを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- getAE(param, type="raw", extract=F)#paramで指定したIDの生データをダウンロード("extract=T"にすると圧縮ファイルを自動的に解凍してくれますが大量のファイル数になります...)
	

param <- "GSE7623"                     #入手したいIDを指定

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

#前処理(データ取得)
hoge <- getAE(param, type="raw", extract=F)#paramで指定したIDの生データをダウンロード("extract=T"にすると圧縮ファイルを自動的に解凍してくれますが大量のファイル数になります...)
	

param <- "GSE7866"                     #入手したいIDを指定

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

#前処理(データ取得)
hoge <- getAE(param, type="raw", extract=F)#paramで指定したIDの生データをダウンロード("extract=T"にすると圧縮ファイルを自動的に解凍してくれますが大量のファイル数になります...)
	

param <- "GSE10417"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- getAE(param, type="raw", extract=F)#paramで指定したIDの生データをダウンロード("extract=T"にすると圧縮ファイルを自動的に解凍してくれますが大量のファイル数になります...)
	

param <- "GSE781"                      #入手したいIDを指定

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

#前処理(データ取得)
hoge <- getAE(param, type="raw", extract=F)#paramで指定したIDの生データをダウンロード("extract=T"にすると圧縮ファイルを自動的に解凍してくれますが大量のファイル数になります...)
	

param <- "GSE7623"                     #入手したいGEO IDを指定

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

#本番
getGEOSuppFiles(param)                 #指定したGEO IDのサプリメントファイル（CELファイルに相当）をダウンロード
	

param1 <- "GSE7623"                    #入手したいGEO IDを指定
param2 <- "hoge2_"                     #出力ファイル名の最初の部分を指定

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

#前処理
eset <- getGEO(param1)                 #指定したGEO IDのデータを取得した結果をesetに格納
sapply(eset, annotation)               #用いられたアレイ情報(GPL ID)を表示

#本番
hoge <- sapply(eset, annotation)       #用いられたアレイ情報(GPL ID)をhogeに格納
for(i in 1:length(hoge)){              #hogeの要素数(用いられたアレイ数)分だけループを回す
    out_f <- paste(param2, hoge[i], ".txt", sep="")#出力ファイル名を作成した結果をout_fに格納
    write.exprs(eset[[i]], file=out_f) #outの中身を指定したファイル名で保存。
}
	

param1 <- "GSE781"                     #入手したいGEO IDを指定
param2 <- "hoge3_"                     #出力ファイル名の最初の部分を指定

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

#前処理
eset <- getGEO(param1)                 #指定したGEO IDのデータを取得した結果をesetに格納
sapply(eset, annotation)               #用いられたアレイ情報(GPL ID)を表示

#本番
hoge <- sapply(eset, annotation)       #用いられたアレイ情報(GPL ID)をhogeに格納
for(i in 1:length(hoge)){              #hogeの要素数(用いられたアレイ数)分だけループを回す
    out_f <- paste(param2, hoge[i], ".txt", sep="")#出力ファイル名を作成した結果をout_fに格納
    write.exprs(eset[[i]], file=out_f) #outの中身を指定したファイル名で保存。
}
	

param1 <- "GSE1322"                    #入手したいGEO IDを指定
param2 <- "hoge4_"                     #出力ファイル名の最初の部分を指定

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

#前処理
eset <- getGEO(param1)                 #指定したGEO IDのデータを取得した結果をesetに格納
sapply(eset, annotation)               #用いられたアレイ情報(GPL ID)を表示

#本番
hoge <- sapply(eset, annotation)       #用いられたアレイ情報(GPL ID)をhogeに格納
for(i in 1:length(hoge)){              #hogeの要素数(用いられたアレイ数)分だけループを回す
    out_f <- paste(param2, hoge[i], ".txt", sep="")#出力ファイル名を作成した結果をout_fに格納
    write.exprs(eset[[i]], file=out_f) #outの中身を指定したファイル名で保存。
}
		

param1 <- "GSE30533"                   #入手したいGEO IDを指定
param2 <- "hoge5_"                     #出力ファイル名の最初の部分を指定

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

#前処理
eset <- getGEO(param1)                 #指定したGEO IDのデータを取得した結果をesetに格納
sapply(eset, annotation)               #用いられたアレイ情報(GPL ID)を表示

#本番
hoge <- sapply(eset, annotation)       #用いられたアレイ情報(GPL ID)をhogeに格納
for(i in 1:length(hoge)){              #hogeの要素数(用いられたアレイ数)分だけループを回す
    out_f <- paste(param2, hoge[i], ".txt", sep="")#出力ファイル名を作成した結果をout_fに格納
    write.exprs(eset[[i]], file=out_f) #outの中身を指定したファイル名で保存。
}
	

out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GPL1355"                     #入手したいGEO IDを指定

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

#前処理
data <- getGEO(param)                  #指定したGEO IDのデータを取得

#本番
out <- data@dataTable@table            #アノテーション情報抽出結果をoutに格納
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
	

param1 <- "GSE7623"                    #入手したいGEO IDを指定
param2 <- "hoge2_"                     #出力ファイル名の最初の部分を指定

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

#前処理
data <- getGEO(param1)                 #指定したGEO IDのデータを取得
sapply(data, annotation)               #用いられたアレイ情報(GPL ID)を表示

#本番
hoge <- sapply(data, annotation)       #用いられたアレイ情報(GPL ID)をhogeに格納
for(i in 1:length(hoge)){              #hogeの要素数(用いられたアレイ数)分だけループを回す
    out_f <- paste(param2, hoge[i], ".txt", sep="")#出力ファイル名を作成した結果をout_fに格納
    out <- data[[i]]@featureData@data  #アノテーション情報抽出結果をoutに格納
    write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
}
	

param1 <- "GSE7623"                    #入手したいGEO IDを指定
param2 <- "hoge3_"                     #出力ファイル名の最初の部分を指定
param_posi <- c(1, 11)                 #outオブジェクト中のID列とgene symbol列の位置情報を指定

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

#前処理
data <- getGEO(param1)                 #指定したGEO IDのデータを取得
sapply(data, annotation)               #用いられたアレイ情報(GPL ID)を表示

#本番
hoge <- sapply(data, annotation)       #用いられたアレイ情報(GPL ID)をhogeに格納
for(i in 1:length(hoge)){              #hogeの要素数(用いられたアレイ数)分だけループを回す
    out_f <- paste(param2, hoge[i], ".txt", sep="")#出力ファイル名を作成した結果をout_fに格納
    out <- data[[i]]@featureData@data  #アノテーション情報抽出結果をoutに格納
    colnames(out)                      #確認してるだけです
    out <- out[,param_posi]            #サブセットを抽出
    write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
}
	

param1 <- "GSE781"                     #入手したいGEO IDを指定
param2 <- "hoge4_"                     #出力ファイル名の最初の部分を指定

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

#前処理
data <- getGEO(param1)                 #指定したGEO IDのデータを取得
sapply(data, annotation)               #用いられたアレイ情報(GPL ID)を表示

#本番
hoge <- sapply(data, annotation)       #用いられたアレイ情報(GPL ID)をhogeに格納
for(i in 1:length(hoge)){              #hogeの要素数(用いられたアレイ数)分だけループを回す
    out_f <- paste(param2, hoge[i], ".txt", sep="")#出力ファイル名を作成した結果をout_fに格納
    out <- data[[i]]@featureData@data  #アノテーション情報抽出結果をoutに格納
    write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
}
	

out_f <- "hoge5.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GPL17586"                    #入手したいGEO IDを指定

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

#前処理
data <- getGEO(param)                  #指定したGEO IDのデータを取得

#本番
out <- data@dataTable@table            #アノテーション情報抽出結果をoutに格納
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
	

param1 <- "GSE1322"                    #入手したいGEO IDを指定
param2 <- "hoge6_"                     #出力ファイル名の最初の部分を指定
param_posi <- c(1, 10)                 #outオブジェクト中のID列とgene symbol列の位置情報を指定

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

#前処理
data <- getGEO(param1)                 #指定したGEO IDのデータを取得
sapply(data, annotation)               #用いられたアレイ情報(GPL ID)を表示

#本番
hoge <- sapply(data, annotation)       #用いられたアレイ情報(GPL ID)をhogeに格納
for(i in 1:length(hoge)){              #hogeの要素数(用いられたアレイ数)分だけループを回す
    out_f <- paste(param2, hoge[i], ".txt", sep="")#出力ファイル名を作成した結果をout_fに格納
    out <- data[[i]]@featureData@data  #アノテーション情報抽出結果をoutに格納
    colnames(out)                      #確認してるだけです
    out <- out[,param_posi]            #サブセットを抽出
    write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
}
	

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

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

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

#前処理(抽出可能な情報のリストアップおよびキーを表示)
hoge <- eval(parse(text=param1))       #オブジェクト名の変更（hogeとして取り扱いたいだけです）
cols(hoge)                             #hoge中に含まれている情報を表示
keytypes(hoge)                         #hogeから実際に使用するキーを表示
param_key <- keys(hoge, keytype="PROBEID")#param1で指定したものをキーにしている

#本番
out <- select(hoge, keys=param_key, keytype="PROBEID", columns=param2)#アノテーション情報抽出結果をoutに格納
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
	

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

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

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

#前処理(抽出可能な情報のリストアップおよびキーを表示)
hoge <- eval(parse(text=param1))       #オブジェクト名の変更（hogeとして取り扱いたいだけです）
cols(hoge)                             #hoge中に含まれている情報を表示
keytypes(hoge)                         #hogeから実際に使用するキーを表示
param_key <- keys(hoge, keytype="PROBEID")#param1で指定したものをキーにしている

#本番
out <- select(hoge, keys=param_key, keytype="PROBEID", columns=param2)#アノテーション情報抽出結果をoutに格納
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
	

out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param1 <- "illuminaMousev2.db"         #パッケージ名を指定
param2 <- c("REFSEQ", "SYMBOL", "GENENAME", "CHRLOC", "ENSEMBLTRANS")#欲しいアノテーション情報を指定（記述可能なリストは以下のkeytypes(dbname)から取得可能）

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

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

#前処理(抽出可能な情報のリストアップおよびキーを表示)
hoge <- eval(parse(text=param1))       #オブジェクト名の変更（hogeとして取り扱いたいだけです）
cols(hoge)                             #hoge中に含まれている情報を表示
keytypes(hoge)                         #hogeから実際に使用するキーを表示
param_key <- keys(hoge, keytype="PROBEID")#param1で指定したものをキーにしている

#本番
out <- select(hoge, keys=param_key, keytype="PROBEID", columns=param2)#アノテーション情報抽出結果をoutに格納
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
	

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

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

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

#前処理(抽出可能な情報のリストアップおよびキーを表示)
hoge <- eval(parse(text=param1))       #オブジェクト名の変更（hogeとして取り扱いたいだけです）
cols(hoge)                             #hoge中に含まれている情報を表示
keytypes(hoge)                         #hogeから実際に使用するキーを表示
param_key <- keys(hoge, keytype="PROBEID")#param1で指定したものをキーにしている

#本番
out <- select(hoge, keys=param_key, keytype="PROBEID", columns=param2)#アノテーション情報抽出結果をoutに格納
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#outの中身を指定したファイル名で保存
	

out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- "hgu133plus2probe"            #パッケージ名を指定

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

#前処理(統一的なオブジェクト名に変更しているだけ)
hoge <- eval(parse(text=param))        #オブジェクト名の変更（hogeとして取り扱いたいだけです）

#本番（ファイルに保存）
write.table(hoge, out_f, sep="\t", append=F, quote=F, row.names=F)#hogeの中身を指定したファイル名で保存
		

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param1 <- "hgu133plus2probe"           #パッケージ名を指定
param2 <- "1552263_at"                 #プローブセットID情報を指定

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

#前処理(統一的なオブジェクト名に変更しているだけ)
hoge <- eval(parse(text=param1))       #オブジェクト名の変更（hogeとして取り扱いたいだけです）

#本番(フィルタリング)
obj <- is.element(as.character(hoge$Probe.Set.Name), param2)#条件を満たすかどうかを判定した結果をobjに格納
out <- hoge[obj,]                      #objがTRUEとなる行のみ抽出した結果をoutに格納

#ファイルに保存
write.table(out, out_f, sep="\t", append=F, quote=F, row.names=F)#hogeの中身を指定したファイル名で保存
		

out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param1 <- "hgu133plus2probe"           #パッケージ名を指定
param2 <- "1552263_at"                 #プローブセットID情報を指定

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

#前処理(統一的なオブジェクト名に変更しているだけ)
hoge <- eval(parse(text=param1))       #オブジェクト名の変更（hogeとして取り扱いたいだけです）

#本番(フィルタリング)
obj <- is.element(as.character(hoge$Probe.Set.Name), param2)#条件を満たすかどうかを判定した結果をobjに格納
out <- hoge[obj,]                      #objがTRUEとなる行のみ抽出した結果をoutに格納

#後処理(FASTA形式に変換)
fasta <- DNAStringSet(out$sequence)    #塩基配列情報をDNAStringSetオブジェクトとしてfastaに格納
description <- paste(out$Probe.Set.Name, out$Probe.Interrogation.Position, sep=".")#description部分を作成
names(fasta) <- description            #description行に相当する記述を追加している

#ファイルに保存
writeXStringSet(fasta, file=out_f, format="fasta", width=50)#fastaの中身を指定したファイル名で保存
		

out_f <- "hoge4.txt"                   #出力ファイル名を指定してout_fに格納
param1 <- "rat2302probe"               #パッケージ名を指定
param2 <- "1369078_at"                 #プローブセットID情報を指定

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

#前処理(統一的なオブジェクト名に変更しているだけ)
hoge <- eval(parse(text=param1))       #オブジェクト名の変更（hogeとして取り扱いたいだけです）

#本番(フィルタリング)
obj <- is.element(as.character(hoge$Probe.Set.Name), param2)#条件を満たすかどうかを判定した結果をobjに格納
out <- hoge[obj,]                      #objがTRUEとなる行のみ抽出した結果をoutに格納

#後処理(FASTA形式に変換)
fasta <- DNAStringSet(out$sequence)    #塩基配列情報をDNAStringSetオブジェクトとしてfastaに格納
description <- paste(out$Probe.Set.Name, out$Probe.Interrogation.Position, sep=".")#description部分を作成
names(fasta) <- description            #description行に相当する記述を追加している

#ファイルに保存
writeXStringSet(fasta, file=out_f, format="fasta", width=50)#fastaの中身を指定したファイル名で保存
		

out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param_dataset <- "rnorvegicus_gene_ensembl"#データセット名を指定(「listDatasets(useMart("ensembl"))」にリストアップされているものを指定可能)
param_attribute <- "affy_rat230_2"     #配列を取得したい属性名を指定(「listAttributes(mart)」にリストアップされているものを指定可能)
param_seqtype <- "cdna"                #配列のタイプを指定(cdna, peptide, 3utr, 5utr, genomicのいずれかが指定可能です)

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

#前処理（IDリスト情報を取得）
mart <- useMart("ensembl",dataset=param_dataset)#データベース名("ensembl")とparam_datasetで指定したデータセット名を与えてmartに格納
hoge <- getBM(attributes=param_attribute, filters="", values="", mart=mart)#指定したIDリスト取得結果をhogeに格納
head(hoge)                             #確認してるだけです
length(hoge)                           #ID数を表示

#本番（配列取得）
out <- getSequence(id=hoge, type=param_attribute, seqType=param_seqtype, mart=mart)#指定したパラメータで配列取得した結果をoutに格納
dim(out)                               #配列数がID数に比べて若干少ないのは複数のプローブセットが同じ転写物由来なためです。

#ファイルに保存
exportFASTA(out, file=out_f)           #outの中身を指定したファイル名で保存
	

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param_dataset <- "rnorvegicus_gene_ensembl"#データセット名を指定(「listDatasets(useMart("ensembl"))」にリストアップされているものを指定可能)
param_attribute <- "affy_rat230_2"     #配列を取得したい属性名を指定(「listAttributes(mart)」にリストアップされているものを指定可能)
param_seqtype <- "5utr"                #配列のタイプを指定(cdna, peptide, 3utr, 5utr, genomicのいずれかが指定可能です)

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

#前処理（IDリスト情報を取得）
mart <- useMart("ensembl",dataset=param_dataset)#データベース名("ensembl")とparam_datasetで指定したデータセット名を与えてmartに格納
hoge <- getBM(attributes=param_attribute, filters="", values="", mart=mart)#指定したIDリスト取得結果をhogeに格納
head(hoge)                             #確認してるだけです
length(hoge)                           #ID数を表示

#本番（配列取得）
out <- getSequence(id=hoge, type=param_attribute, seqType=param_seqtype, mart=mart)#指定したパラメータで配列取得した結果をoutに格納
dim(out)                               #配列数を表示

#ファイルに保存
exportFASTA(out, file=out_f)           #outの中身を指定したファイル名で保存
	

out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param_dataset <- "rnorvegicus_gene_ensembl"#データセット名を指定(「listDatasets(useMart("ensembl"))」にリストアップされているものを指定可能)
param_attribute <- "affy_rat230_2"     #配列を取得したい属性名を指定(「listAttributes(mart)」にリストアップされているものを指定可能)
param_seqtype <- "cdna"                #配列のタイプを指定(cdna, peptide, 3utr, 5utr, genomicのいずれかが指定可能です)
param_value <- "1369078_at"           #プローブセットIDを指定

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

#前処理
mart <- useMart("ensembl",dataset=param_dataset)#データベース名("ensembl")とparam_datasetで指定したデータセット名を与えてmartに格納

#本番（配列取得）
out <- getSequence(id=param_value, type=param_attribute, seqType=param_seqtype, mart=mart)#指定したパラメータで配列取得した結果をoutに格納
dim(out)                               #配列数がID数に比べて若干少ないのは複数のプローブセットが同じ転写物由来なためです。

#ファイルに保存
exportFASTA(out, file=out_f)           #outの中身を指定したファイル名で保存
	

out_f <- "hoge4.txt"                   #出力ファイル名を指定してout_fに格納
param_dataset <- "rnorvegicus_gene_ensembl"#データセット名を指定(「listDatasets(useMart("ensembl"))」にリストアップされているものを指定可能)
param_attribute <- "ensembl_gene_id"   #配列を取得したい属性名を指定(「listAttributes(mart)」にリストアップされているものを指定可能)
param_seqtype <- "cdna"                #配列のタイプを指定(cdna, peptide, 3utr, 5utr, genomicのいずれかが指定可能です)

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

#前処理（IDリスト情報を取得）
mart <- useMart("ensembl",dataset=param_dataset)#データベース名("ensembl")とparam_datasetで指定したデータセット名を与えてmartに格納
hoge <- getBM(attributes=param_attribute, filters="", values="", mart=mart)#指定したIDリスト取得結果をhogeに格納
head(hoge)                             #確認してるだけです
length(hoge)                           #ID数を表示

#本番（配列取得）
out <- getSequence(id=hoge, type=param_attribute, seqType=param_seqtype, mart=mart)#指定したパラメータで配列取得した結果をoutに格納
dim(out)                               #配列数を表示

#ファイルに保存
exportFASTA(out, file=out_f)           #outの中身を指定したファイル名で保存
	

out_f <- "hoge5.txt"                   #出力ファイル名を指定してout_fに格納
param_dataset <- "hsapiens_gene_ensembl"#データセット名を指定(「listDatasets(useMart("ensembl"))」にリストアップされているものを指定可能)
param_attribute <- "ensembl_gene_id"   #配列を取得したい属性名を指定(「listAttributes(mart)」にリストアップされているものを指定可能)
param_seqtype <- "cdna"                #配列のタイプを指定(cdna, peptide, 3utr, 5utr, genomicのいずれかが指定可能です)

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

#前処理（IDリスト情報を取得）
mart <- useMart("ensembl",dataset=param_dataset)#データベース名("ensembl")とparam_datasetで指定したデータセット名を与えてmartに格納
hoge <- getBM(attributes=param_attribute, filters="", values="", mart=mart)#指定したIDリスト取得結果をhogeに格納
head(hoge)                             #確認してるだけです
length(hoge)                           #ID数を表示

#本番（配列取得）
out <- getSequence(id=hoge, type=param_attribute, seqType=param_seqtype, mart=mart)#指定したパラメータで配列取得した結果をoutに格納

#ファイルに保存
exportFASTA(out, file=out_f)           #outの中身を指定したファイル名で保存
	

out_f <- "hoge6.txt"                   #出力ファイル名を指定してout_fに格納
param_dataset <- "hsapiens_gene_ensembl"#データセット名を指定(「listDatasets(useMart("ensembl"))」にリストアップされているものを指定可能)
param_attribute <- "refseq_mrna"       #配列を取得したい属性名を指定(「listAttributes(mart)」にリストアップされているものを指定可能)
param_seqtype <- "cdna"                #配列のタイプを指定(cdna, peptide, 3utr, 5utr, genomicのいずれかが指定可能です)

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

#前処理（IDリスト情報を取得）
mart <- useMart("ensembl",dataset=param_dataset)#データベース名("ensembl")とparam_datasetで指定したデータセット名を与えてmartに格納
hoge <- getBM(attributes=param_attribute, filters="", values="", mart=mart)#指定したIDリスト取得結果をhogeに格納
head(hoge)                             #確認してるだけです
length(hoge)                           #ID数を表示

#本番（配列取得）
out <- getSequence(id=hoge, type=param_attribute, seqType=param_seqtype, mart=mart)#指定したパラメータで配列取得した結果をoutに格納

#ファイルに保存
exportFASTA(out, file=out_f)           #outの中身を指定したファイル名で保存
	

out_f <- "hoge7.txt"                   #出力ファイル名を指定してout_fに格納
param_dataset <- "hsapiens_gene_ensembl"#データセット名を指定(「listDatasets(useMart("ensembl"))」にリストアップされているものを指定可能)
param_attribute <- "refseq_mrna"       #配列を取得したい属性名を指定(「listAttributes(mart)」にリストアップされているものを指定可能)
param_seqtype <- "cdna"                #配列のタイプを指定(cdna, peptide, 3utr, 5utr, genomicのいずれかが指定可能です)
param_value <- "NM_138957"             #IDを指定

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

#前処理
mart <- useMart("ensembl",dataset=param_dataset)#データベース名("ensembl")とparam_datasetで指定したデータセット名を与えてmartに格納

#本番（配列取得）
out <- getSequence(id=param_value, type=param_attribute, seqType=param_seqtype, mart=mart)#指定したパラメータで配列取得した結果をoutに格納

#ファイルに保存
exportFASTA(out, file=out_f)           #outの中身を指定したファイル名で保存
	

param <- "NM_031237"                   #Genbank accession番号を指定

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

#本番（配列取得）
out <- getSEQ(param)                   #指定したパラメータで配列取得した結果をoutに格納
out                                    #確認してるだけです
	

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- "NM_031237"                   #Genbank accession番号を指定

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

#本番（配列取得）
out <- getSEQ(param)                   #指定したパラメータで配列取得した結果をoutに格納
out                                    #確認してるだけです

#ファイルに保存
write.fasta(out, names=param, file.out=out_f, nbchar=50)#outの中身を指定したファイル名で保存
	

param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納
hoge                                   #確認してるだけです
sampleNames(hoge)                      #読み込んだCELファイル名を表示

#本番(各種情報取得)
PM <- pm(hoge)                         #PMプローブシグナル強度情報を取得した結果をPMに格納
head(PM)                               #確認してるだけです
MM <- mm(hoge)                         #MMプローブシグナル強度情報を取得した結果をMMに格納
head(MM)                               #確認してるだけです

length(PM[,1])                         #このアレイ上に搭載されているPMプローブ数を表示
length(MM[,1])                         #このアレイ上に搭載されているMMプローブ数を表示
	

param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番(PM-MMシグナルの値がマイナスになるものがどれくらいの割合を占めるのか調べる)
PM <- pm(hoge)                         #PMプローブシグナル強度情報を取得した結果をPMに格納
MM <- mm(hoge)                         #MMプローブシグナル強度情報を取得した結果をMMに格納

out <- PM - MM                         #PM - MM計算結果をoutに格納
length(out)                            #行列outの全要素数を表示
sum(out > 0)                           #(PM - MM) > 0となる要素数を表示
sum(out <= 0)                          #(PM - MM) <= 0となる要素数を表示
sum(out <= 0) / length(out)            #(PM - MM) <= 0となった割合を表示
mean(MM >= PM)                         #(MM >= PM)となった割合を表示
mean(MM > PM)                          #(MM > PM)となった割合を表示
	

param1 <- "GSE30533"                   #入手したいGSE IDを指定
param2 <- "1369078_at"                 #プローブセットIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param1, save=F)   #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納
sampleNames(hoge)                      #確認してるだけです(サンプル名およびその順番を)
hoge <- hoge[,1]                       #一つめのサンプルの情報のみ抽出した結果をhogeに格納
sampleNames(hoge)                      #確認してるだけです

#本番(目的のprobeset IDのシグナル強度情報のみ抽出)
PM <- pm(hoge, param2)                 #PMプローブシグナル強度情報を取得した結果をPMに格納
MM <- mm(hoge, param2)                 #MMプローブシグナル強度情報を取得した結果をMMに格納
PM                                     #確認してるだけです
MM                                     #確認してるだけです
	

out_f <- "hoge4.png"                   #出力ファイル名を指定してout_fに格納
param1 <- "GSE30533"                   #入手したいGSE IDを指定
param2 <- "1369078_at"                 #プローブセットIDを指定
param_fig <- c(500, 300)               #ファイル出力時の横幅と縦幅を指定(単位はピクセル)

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

#前処理(データ取得)
hoge <- ArrayExpress(param1, save=F)   #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納
hoge <- hoge[,1]                       #一つめのサンプルの情報のみ抽出した結果をhogeに格納
sampleNames(hoge)                      #確認してるだけです

#本番(目的のprobeset IDのシグナル強度情報のみ抽出)
PM <- pm(hoge, param2)                 #PMプローブシグナル強度情報を取得した結果をPMに格納
MM <- mm(hoge, param2)                 #MMプローブシグナル強度情報を取得した結果をMMに格納

#ファイルに保存
png(out_f, pointsize=13, width=param_fig[1], height=param_fig[2])#出力ファイルの各種パラメータを指定
min <- min(c(PM, MM))                  #y軸の範囲を決めるべく、PMとMMの数値ベクトル中の最小値を取得している
max <- max(c(PM, MM))                  #y軸の範囲を決めるべく、PMとMMの数値ベクトル中の最大値を取得している
ylim <- c(min, max)                    #y軸の範囲情報をylimに格納
plot(PM, ylim=ylim, ylab="Intensity",  #PMシグナル強度をプロット
    pch=20, type="o", col="black")     #PMシグナル強度をプロット
par(new=T)                             #図の重ね書き指定
plot(MM, ylim=ylim, ylab="",           #MMシグナル強度をプロット
    pch=20, type="o", col="gray")      #MMシグナル強度をプロット
title(param2)                          #図のタイトルを追加
legend("topright", c("PM", "MM"),      #図の凡例を追加
    col=c("black", "gray"), pch=20)    #図の凡例を追加
dev.off()                              #おまじない
	

out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- "robust_weighted_average"     #summarization法を指定

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- frma(hoge, summarize=param)    #fRMAを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- "median_polish"               #summarization法を指定

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- frma(hoge, summarize=param)    #fRMAを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- frma(hoge, summarize=param)    #fRMAを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- robloxbioc(hoge)               #rmxを実行し、結果をesetに保存

#対数変換
summary(exprs(eset))                   #得られたesetの遺伝子発現行列のシグナル強度分布を表示
exprs(eset)[exprs(eset) < 1] <- 1      #対数変換（log2）できるようにシグナル強度が1未満のものを1にしておく
summary(exprs(eset))                   #上記処理後のシグナル強度分布を再び表示させて確認
exprs(eset) <- log(exprs(eset), 2)     #底を2として対数変換

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- robloxbioc(hoge)               #rmxを実行し、結果をesetに保存

#対数変換
summary(exprs(eset))                   #得られたesetの遺伝子発現行列のシグナル強度分布を表示
exprs(eset)[exprs(eset) < 1] <- 1      #対数変換（log2）できるようにシグナル強度が1未満のものを1にしておく
summary(exprs(eset))                   #上記処理後のシグナル強度分布を再び表示させて確認
exprs(eset) <- log(exprs(eset), 2)     #底を2として対数変換

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納
hoge                                   #確認してるだけです
sampleNames(hoge)                      #読み込んだCELファイル名を表示

#本番(最初の2サンプルの正規化)
posi <- 1:2
eset <- robloxbioc(hoge[posi])         #rmxを実行し、結果をesetに保存
head(exprs(eset))                      #確認してるだけです

#本番(最初の3サンプルの正規化)
posi <- 1:3
eset <- robloxbioc(hoge[posi])         #rmxを実行し、結果をesetに保存
head(exprs(eset))                      #確認してるだけです
	

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

#必要な関数などをロード
source("http://www.iu.a.u-tokyo.ac.jp/~kadota/R/dfw.R")#DFWを実行する関数を含むファイルをあらかじめ読み込む
library(affy)                          #パッケージの読み込み

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #CELファイルの読み込み

#本番
eset <- expresso(hoge, bgcorrect.method="none", normalize.method="quantiles", pmcorrect.method="pmonly", summary.method="dfw")#DFW法を実行し、結果をesetに保存。

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
		

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- q.farms(hoge)                  #qFARMSを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
		

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- l.farms(hoge)                  #lFARMSを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
		

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- mmgmos(hoge)                   #multi-mgMOSを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- mmgmos(hoge)                   #multi-mgMOSを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- gcrma(hoge)                    #GCRMAを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- gcrma(hoge)                    #GCRMAを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- justPlier(hoge)                #PLIERを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
		

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- justPlier(hoge, normalize=TRUE)#quantile normalizationつきのPLIERを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
		

out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- justPlier(hoge)                #PLIERを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
		

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- vsnrma(hoge)                   #VSNを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- vsnrma(hoge)                   #VSNを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- rma(hoge)                      #RMAを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- rma(hoge)                      #RMAを実行し、結果をesetに保存

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- mas5(hoge)                     #MASを実行し、結果をesetに保存

#対数変換
summary(exprs(eset))                   #得られたesetの遺伝子発現行列のシグナル強度分布を表示
exprs(eset)[exprs(eset) < 1] <- 1      #対数変換（log2）できるようにシグナル強度が1未満のものを1にしておく
summary(exprs(eset))                   #上記処理後のシグナル強度分布を再び表示させて確認
exprs(eset) <- log(exprs(eset), 2)     #底を2として対数変換

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- mas5(hoge)                     #MASを実行し、結果をesetに保存

#対数変換
summary(exprs(eset))                   #得られたesetの遺伝子発現行列のシグナル強度分布を表示
exprs(eset)[exprs(eset) < 1] <- 1      #対数変換（log2）できるようにシグナル強度が1未満のものを1にしておく
summary(exprs(eset))                   #上記処理後のシグナル強度分布を再び表示させて確認
exprs(eset) <- log(exprs(eset), 2)     #底を2として対数変換

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- mas5(hoge)                     #MASを実行し、結果をesetに保存

#対数変換
summary(exprs(eset))                   #得られたesetの遺伝子発現行列のシグナル強度分布を表示
exprs(eset) <- log(exprs(eset), 2)     #底を2として対数変換

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- expresso(hoge, normalize.method="invariantset", bg.correct=FALSE,#dChip(PM onlyモデル)を実行し、結果をesetに保存
                 pmcorrect.method="pmonly", summary.method="liwong")#dChip(PM onlyモデル)を実行し、結果をesetに保存

#対数変換
summary(exprs(eset))                   #得られたesetの遺伝子発現行列のシグナル強度分布を表示
exprs(eset)[exprs(eset) < 1] <- 1      #対数変換（log2）できるようにシグナル強度が1未満のものを1にしておく
summary(exprs(eset))                   #上記処理後のシグナル強度分布を再び表示させて確認
exprs(eset) <- log(exprs(eset), 2)     #底を2として対数変換

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

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

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

#データファイルの読み込み(*.CELファイル)
hoge <- ReadAffy()                     #*.CELファイルの読み込み

#本番
eset <- expresso(hoge, normalize.method="invariantset", bg.correct=FALSE,#dChip(PM-MMモデル)を実行し、結果をesetに保存
                 pmcorrect.method="subtractmm", summary.method="liwong")#dChip(PM-MMモデル)を実行し、結果をesetに保存

#対数変換
summary(exprs(eset))                   #得られたesetの遺伝子発現行列のシグナル強度分布を表示
exprs(eset)[exprs(eset) < 1] <- 1      #対数変換（log2）できるようにシグナル強度が1未満のものを1にしておく
summary(exprs(eset))                   #上記処理後のシグナル強度分布を再び表示させて確認
exprs(eset) <- log(exprs(eset), 2)     #底を2として対数変換

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納
param <- "GSE30533"                    #入手したいIDを指定

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

#前処理(データ取得)
hoge <- ArrayExpress(param, save=F)    #paramで指定したIDのCELファイルなどを取得した結果をhogeに格納

#本番
eset <- expresso(hoge, normalize.method="invariantset", bg.correct=FALSE,#dChip(PM-MMモデル)を実行し、結果をesetに保存
                 pmcorrect.method="subtractmm", summary.method="liwong")#dChip(PM-MMモデル)を実行し、結果をesetに保存

#対数変換
summary(exprs(eset))                   #得られたesetの遺伝子発現行列のシグナル強度分布を表示
exprs(eset)[exprs(eset) < 1] <- 1      #対数変換（log2）できるようにシグナル強度が1未満のものを1にしておく
summary(exprs(eset))                   #上記処理後のシグナル強度分布を再び表示させて確認
exprs(eset) <- log(exprs(eset), 2)     #底を2として対数変換

#ファイルに保存
write.exprs(eset, file=out_f)          #結果を指定したファイル名で保存
	

in_f <- "sample19.txt"                 #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param <- 10                            #各サンプルの正規化後の平均値を指定

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

#本番
tmp_mean <- apply(data_tmp, 2, mean, na.rm=TRUE)#各サンプル（列）の平均シグナル強度を計算した結果をtmp_meanに格納
data <- sweep(data_tmp, 2, param/tmp_mean, "*")#各列中の全てのシグナル値にparam/tmp_meanを掛け、その結果をdataに格納
data                                   #結果を表示
apply(data, 2, mean)                   #各列のmeanを表示させ、正常に動作しているか確認

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

in_f <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge.txt"                    #出力ファイル名を指定してout_fに格納
param <- 100                           #各サンプルの正規化後の平均値を指定

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

#本番
tmp_mean <- apply(data_tmp, 2, mean, na.rm=TRUE)#各サンプル（列）の平均シグナル強度を計算した結果をtmp_meanに格納
data <- sweep(data_tmp, 2, param/tmp_mean, "*")#各列中の全てのシグナル値にparam/tmp_meanを掛け、その結果をdataに格納
data                                   #結果を表示
apply(data, 2, mean)                   #各列のmeanを表示させ、正常に動作しているか確認

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

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

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

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

#本番
data.z <- normalize(data, byrow=FALSE) #列方向に正規化した結果をdata.zに格納
data.z                                 #結果を表示
apply(data.z, 2, mean)                 #各列のmeanを表示させ、正常に動作しているか確認

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

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

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

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

#本番
data.z <- genescale(data, axis=2, method="Z")#列方向に正規化した結果をdata.zに格納
data.z                                 #結果を表示
apply(data.z, 2, mean)                 #各列のmeanを表示させ、正常に動作しているか確認

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

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

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

data.z <- scale(data)                  #列方向に正規化した結果をdata.zに格納
data.z                                 #結果を表示
apply(data.z, 2, mean)                 #各列のmeanを表示させ、正常に動作しているか確認

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

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

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

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

#本番
data.z <- normalize(data, byrow=TRUE)  #行方向に正規化した結果をdata.zに格納
data.z                                 #結果を表示
apply(data.z, 1, mean)                 #各行のmeanを表示させ、正常に動作しているか確認

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

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

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

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

#本番
data.z <- genescale(data, axis=1, method="Z")#行方向に正規化した結果をdata.zに格納
data.z                                 #結果を表示
apply(data.z, 1, mean)                 #各行のmeanを表示させ、正常に動作しているか確認

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

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

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

#本番
tmp <- scale(t(data))                  #scale関数は列方向のscalingしかしてくれないので、t関数を使って転置行列に対してscale関数を実行し、結果をtmpに格納
data.z <- t(tmp)                       #前の行で行列を転置させていたので、もう一度転置しなおした結果をdata.zに格納
data.z                                 #結果を表示
apply(data.z, 1, mean)                 #各行のmeanを表示させ、正常に動作しているか確認

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

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

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#入力ファイルを読み込んでdataに格納

#本番
data.t <- t(data)                      #行列の転置（この後に使うscale関数が列ごとの操作を行うため、scale関数をそのまま使えるように予め行列を入れ替えておく必要があるため）
data.t.mean.mad <- scale(data.t, apply(data.t,2,mean), apply(data.t,2,mad,constant=1))#各列のmean=0, MAD=1になるようにスケーリングし、結果をdata.t.mean.madに格納
data.mean.mad <- t(data.t.mean.mad)    #scale関数の適用が終わったので、もう一度行列の転置を行って元に戻し、結果をdata.mean.madに格納
apply(data.mean.mad, 1, mean)          #各行のmeanを表示させ、正常に動作しているか確認

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

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

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

#本番
data.mean.mad <- scale(data, apply(data,2,mean), apply(data,2,mad,constant=1))#各列のmean=0, MAD=1になるようにスケーリングし、結果をdata.mean.madに格納
apply(data.mean.mad, 2, mean)          #各列のmeanを表示させ、正常に動作しているか確認

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

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

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

#本番
data.t <- t(data)                      #行列の転置（この後に使うscale関数が列ごとの操作を行うため、scale関数をそのまま使えるように予め行列を入れ替えておく必要があるため）
data.t.m.sd <- scale(data.t, apply(data.t,2,median), apply(data.t,2,sd))#各列のmedian=0, SD=1になるようにスケーリングし、結果をdata.t.m.sdに格納
data.m.sd <- t(data.t.m.sd)            #scale関数の適用が終わったので、もう一度行列の転置を行って元に戻し、結果をdata.m.sdに格納
apply(data.m.sd, 1, median)            #各行のmedianを表示させ、正常に動作しているか確認

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

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

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

#本番
data.m.sd <- scale(data, apply(data,2,median), apply(data,2,sd))#各列のmedian=0, SD=1になるようにスケーリングし、結果をdata.m.sdに格納
apply(data.m.sd, 2, median)            #各列のmedianを表示させ、正常に動作しているか確認

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

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

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

#本番
data.t <- t(data)                      #行列の転置（この後に使うscale関数が列ごとの操作を行うため、scale関数をそのまま使えるように予め行列を入れ替えておく必要があるため）
data.t.m.mad <- scale(data.t, apply(data.t,2,median), apply(data.t,2,mad,constant=1))#各列のmedian=0, MAD=1になるようにスケーリングし、結果をdata.t.m.madに格納
data.m.mad <- t(data.t.m.mad)          #scale関数の適用が終わったので、もう一度行列の転置を行って元に戻し、結果をdata.m.madに格納
apply(data.m.mad, 1, median)           #各行のmedianを表示させ、正常に動作しているか確認

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

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

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

#本番
data.m.mad <- scale(data, apply(data,2,median), apply(data,2,mad,constant=1))#各列のmedian=0, MAD=1になるようにスケーリングし、結果をdata.m.madに格納
apply(data.m.mad, 2, median)           #各列のmedianを表示させ、正常に動作しているか確認

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

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

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

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

#本番
data.t <- t(data)                      #行列の転置（この後に使うscale関数が列ごとの操作を行うため、scale関数をそのまま使えるように予め行列を入れ替えておく必要があるため）
data.t.tukey.mad <- scale(data.t,apply(data.t,2,tukey.biweight),apply(data.t,2,mad,constant=1))#各列のTukey's biweight=0, MAD=1になるようにスケーリングし、結果をdata.t.tukey.madに格納
data.tukey.mad <- t(data.t.tukey.mad)  #scale関数の適用が終わったので、もう一度行列の転置を行って元に戻し、結果をdata.tukey.madに格納
apply(data.tukey.mad, 1, tukey.biweight)#各行のTukey's biweightを表示させ、正常に動作しているか確認

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

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

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

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#入力ファイルを読み込んでdataに格納

data.tukey.mad <- scale(data,apply(data,2,tukey.biweight), apply(data,2,mad,constant=1))#各列のTukey's biweight=0, MAD=1になるようにスケーリングし、結果をdata.tukey.madに格納
apply(data.tukey.mad, 2, tukey.biweight)#各列のTukey's biweightを表示させ、正常に動作しているか確認

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

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

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

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

#本番
data.r <- genescale(data, axis=1, method="R")#スケーリングした結果をdata.rに格納
apply(data.r, 1, range)                #各行の最小値と最大値を表示させ、正常に動作しているか確認

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

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

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

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

#本番
data.r <- genescale(data, axis=2, method="R")#スケーリングした結果をdata.rに格納
apply(data.r, 2, range)                #各列の最小値と最大値を表示させ、正常に動作しているか確認

#ファイルに保存
tmp <- cbind(rownames(data.r), data.r) #遺伝子IDの列を行列data.rの左端に挿入し、結果を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で指定したファイルの読み込み

#本番
summary(data)                          #シグナル強度の分布を確認し、ダイナミックレンジが4桁程度あることを確認
data[data < 1] <- 1                    #シグナル強度が1未満のものを1にする
data.log <- log(data, base=2)          #log2-transformed dataをdata.logに格納
summary(data.log)                      #対数変換後のシグナル強度の分布を確認

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

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

#データファイルの読み込み
data <- read.table(in_f, header=TRUE, row.names=1, sep="\t", quote="")#in_fで指定したファイルの読み込み
data$IDENTIFIER <- NULL                #余分なIDENTIFIER列の消去
summary(data)                          #シグナル強度の分布を確認し、ダイナミックレンジが4桁程度あることを確認

#本番
data[data < 1] <- 1                    #シグナル強度が1未満のものを1にする
data.log <- log(data, base=2)          #log2-transformed dataをdata.logに格納
summary(data.log)                      #対数変換後のシグナル強度の分布を確認

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

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

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

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

#本番
data.q <- normalizeBetweenArrays(as.matrix(data))#Quantile正規化を実行し、結果をdata.qに格納

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

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

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

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

#本番
data.q <- normalizeBetweenArrays(as.matrix(data))#Quantile正規化を実行し、結果をdata.qに格納

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

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

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

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

#本番
data.q <- normalizeBetweenArrays(as.matrix(data))#Quantile正規化を実行し、結果をdata.qに格納

#ファイルに保存
tmp <- cbind(rownames(data.q), data.q) #遺伝子IDの列を行列data.qの左端に挿入し、結果を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に格納
param1 <- 5                            #組織数kを指定
param2 <- 50                           #シグナル強度Aを指定

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

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

#本番
f1 <- kOverA(param1, A=param2)         #「k組織以上でシグナル強度A以上を持つ遺伝子を抽出」という条件（filter）をf1に格納
ffun <- filterfun(f1)                  #フィルタリング用の関数(filtering function)を作成しffunに格納
obj <- genefilter(data, ffun)          #条件を満たすかどうかを判定した結果をobjに格納
obj                                    #確認してるだけです
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
data                                   #確認してるだけです

#ファイルに保存
tmp <- cbind(rownames(data), data)     #保存したい情報を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に格納
param1 <- 0.7                          #組織数の割合pを指定
param2 <- 20                           #シグナル強度Aを指定

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

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

#本番
f1 <- pOverA(p=param1, A=param2)       #「シグナル強度がAよりも大きい組織数の割合がpよりも大きい遺伝子を抽出」という条件（filter）をf1に格納
ffun <- filterfun(f1)                  #フィルタリング用の関数(filtering function)を作成しffunに格納
obj <- genefilter(data, ffun)          #条件を満たすかどうかを判定した結果をobjに格納
obj                                    #確認してるだけです
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
data                                   #確認してるだけです

#ファイルに保存
tmp <- cbind(rownames(data), data)     #保存したい情報を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 <- 92                            #シグナル強度Aを指定

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

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

#本番
f1 <- maxA(param)                      #「少なくとも一つの組織でシグナル強度がparam以上の遺伝子を抽出」という条件（filter）をf1に格納
ffun <- filterfun(f1)                  #フィルタリング用の関数(filtering function)を作成しffunに格納
obj <- genefilter(data, ffun)          #条件を満たすかどうかを判定した結果をobjに格納
obj                                    #確認してるだけです
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
data                                   #確認してるだけです

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

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

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

#本番
ave <- rowMeans(data)                  #行ごとの平均値を算出した結果をaveに格納
hoge <- quantile(ave, probs=param)     #paramで指定した下位の位置情報をhogeに格納
obj <- (ave >= hoge)                   #条件を満たすかどうかを判定した結果を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 <- "sample2.txt"                  #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納
param <- 2.9                           #シグナル強度最大と最小の発現比を指定

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

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

#本番
data <- filtering(data,lt=min(data),ut=max(data),mmr=param,mmd=0)#フィルタリングを実行し、結果をdataに格納
data                                   #確認してるだけです

#ファイルに保存
tmp <- cbind(rownames(data), data)     #保存したい情報を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に格納
param1 <- 10                           #シグナル強度の下限を指定
param2 <- 100                          #シグナル強度の上限を指定
param3 <- 2.9                          #シグナル強度最大と最小の発現比を指定
param4 <- 42                           #シグナル強度最大と最小の差を指定

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

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

#本番
data <- filtering(data, lt=param1, ut=param2, mmr=param3, mmd=param4)#フィルタリングを実行し、結果をdataに格納
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 <- 2.0                           #シグナル強度最大と最小の発現比を指定(対数変換前の値)

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

#本番
max <- apply(data, 1, max)             #行ごとのシグナル強度の最大値をmaxに格納
min <- apply(data, 1, min)             #行ごとのシグナル強度の最小値をminに格納
logratio <- max - min                  #フィルタリングを実行した結果をobjに格納
obj <- (logratio >= log2(param))       #条件を満たすかどうかを判定した結果をobjに格納
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
summary(logratio)                      #確認してるだけです
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 <- "hoge4.txt"                   #出力ファイル名を指定してout_fに格納
param <- 0.1                           #フィルタリングしたい下位x%を指定

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

#本番
max <- apply(data, 1, max)             #行ごとのシグナル強度の最大値をmaxに格納
min <- apply(data, 1, min)             #行ごとのシグナル強度の最小値をminに格納
logratio <- max - min                  #フィルタリングを実行した結果をobjに格納
hoge <- quantile(logratio, probs=param)#paramで指定した下位の位置情報をhogeに格納
obj <- (logratio >= hoge)              #条件を満たすかどうかを判定した結果をobjに格納
data <- data[obj,]                     #objがTRUEとなる要素のみ抽出した結果をdataに格納
summary(logratio)                      #確認してるだけです
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 <- "sample13_7vs7.txt"            #入力ファイル名を指定してin_fに格納
out_f <- "hoge1.txt"                   #出力ファイル名を指定してout_fに格納

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

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

#本番
obj <- apply(data, 1, allNA)           #条件を満たすかどうかを判定した結果をobjに格納(全ての要素がNAならFALSE,それ以外ならTRUE)
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 <- "sample13_7vs7.txt"            #入力ファイル名を指定してin_fに格納
out_f <- "hoge2.txt"                   #出力ファイル名を指定してout_fに格納

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

#本番
hoge <- t(apply(data, 1, is.na))       #各要素がNAであるかどうか(TRUE or FALSE)を返すis.na関数を各行に対して適用した結果をhogeに格納
hoge2 <- apply(hoge, 1, sum)           #TRUEの数(NAの数)を返すsum関数を各行に対して適用し、結果をhoge2に格納
obj <- hoge2 < ncol(data)              #条件を満たすかどうかを判定した結果を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 <- "sample13_7vs7.txt"            #入力ファイル名を指定してin_fに格納
out_f <- "hoge3.txt"                   #出力ファイル名を指定してout_fに格納

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

#本番
hoge <- t(apply(data, 1, is.na))       #各要素がNAであるかどうか(TRUE or FALSE)を返すis.na関数を各行に対して適用した結果をhogeに格納
hoge2 <- apply(hoge, 1, sum)           #TRUEの数(NAの数)を返すsum関数を各行に対して適用し、結果をhoge2に格納
obj <- (hoge2 == 0)                    #条件を満たすかどうかを判定した結果を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 <- "sample13_7vs7.txt"            #入力ファイル名を指定してin_fに格納
out_f <- "hoge4.txt"                   #出力ファイル名を指定してout_fに格納
param_G1 <- 7                          #G1群のサンプル数を指定
param_G2 <- 7                          #G2群のサンプル数を指定
param1 <- 5                            #G1群のNA数の閾値を指定
param2 <- 3                            #G2群のNA数の閾値を指定

#入力ファイルの読み込みとラベル情報の作成
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を作成

#前処理(NAの要素数を各行についてカウント)と本番
hoge <- t(apply(data[data.cl==1],1,is.na))#各要素がNAであるかどうかを返すis.na関数をG1群のみに対して適用した結果をhogeに格納
obj1 <- apply(hoge,1,sum) < param1     #条件を満たすかどうかを判定した結果をobj1に格納(G1群)
hoge <- t(apply(data[data.cl==2],1,is.na))#各要素がNAであるかどうかを返すis.na関数をG2群のみに対して適用した結果をhogeに格納
obj2 <- apply(hoge,1,sum) < param2     #条件を満たすかどうかを判定した結果をobj2に格納(G2群)
obj <- (obj1 & obj2)                   #条件を満たすかどうかを判定した結果を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の中身を指定したファイル名で保存