实验思路：

程序代码及运行结果：

library(affy)%把affy包载入R中

library(tcltk)%把tcltk包载入R中

filters <-matrix(c("CEL file", ".[Cc][Ee][Ll]", "All", ".*"), ncol = 2, byrow = T)%生成2*2字符矩阵，用来定义选择文件类型

cel.files <-tk_choose.files(caption = "Select CELs", multi = TRUE,filters = filters, index = 1)%跳出选择文件窗口并选择文件

data.raw <-ReadAffy(filenames = cel.files)%读入所有文件的文件名

n.cel <-length(cel.files)%显示出所选文件的数量

n.cel

#[1] 8

sampleNames(data.raw)%查看文件名

#[1] "GSM49953.CEL" "GSM49954.CEL" "GSM49955.CEL"

#[4]"GSM49956.CEL" "GSM49957.CEL" "GSM50017.CEL"

#[7] "GSM50018.CEL" "GSM50019.CEL"

sampleNames(data.raw)  <-paste("S",1:n.cel, sep='')%为了直观将文件名按顺序依次命名为S1，S2...Sn

pData(data.raw)$treatment <-rep(c("0h", "1h", "24h", "7d"), each=2)%生成0h，1h，24h，7d四个值依次重复两次所组成的数列，数列命名为treatment

pData(data.raw)%指针pData函数读取文件

#       sample treatment

#   S1      1        0h

#   S2      2        0h

#   S3      3        1h

#   S4      4        1h

#   S5      5       24h

#   S6      6       24h

#   S7      7        7d

#   S8      8        7d

%1、计算基因表达量

eset.rma <-rma(data.raw)%调用RMA算法对数据进行预处理

#Background correcting

#Normalizing

#Calculating Expression

eset.mas5 <-mas5(data.raw)%调用mas5算法对数据进行预处理

#background correction: mas

#PM/MM correction : mas

#expression values: mas

#background correcting...done.

#22283 ids to be processed

#|                    |

#|####################|

emat.rma.log2 <-exprs(eset.rma)%rma的eset结果是经过对数变换的

emat.mas5.nologs <-exprs(eset.mas5)%mas5的eset结果是原始信号强度

class(emat.rma.log2)%结果是matrix格式

#[1] "matrix"

emat.rma.nologs <-2^emat.rma.log2%虽然表达量是用对数变换的信号值表示的，但是有些计算过程要用到未经变换的原始值，应该把它们都计算出来

emat.mas5.log2 <-log2(emat.mas5.nologs)%对mas5的原始信号强度进行对数处理

rm(eset.mas5)%这里仅使用rma的结果做演示。计算平均表达量和差异表达倍数（和0h对照比)rm函数能永久地从workspace中删除mas5的eset结果

rm(emat.mas5.nologs)%删除未经变换的mas5的eset结果

rm(emat.mas5.log2)%删除经过对数变换的mas5的eset结果

results.rma <-data.frame((emat.rma.log2[,c(1,3,5,7)] + emat.rma.log2[,c(2,4,6,8)])/2)%计算平均值，并转换为数据框格式#计算表达量差异倍数

results.rma$fc.1h <- results.rma[,2]-results.rma[,1]

results.rma$fc.24h <- results.rma[,3]-results.rma[,1]

results.rma$fc.7d <- results.rma[,4]-results.rma[,1]

subset.logic <- results.rma$fc.7d>0

subset.data <- results.rma[subset.logic,]%要注意的是逻辑向量的长度要和相应维度的数据长度一致，逻辑向量中为TRUE的就保留，FALSE的就丢弃。

length(subset.logic); nrow(results.rma)%显示逻辑向量的长度

#[1] 22283

#[1] 22283

head(subset.logic)%显示逻辑向量

#[1] FALSE  TRUE  TRUE  TRUE  TRUE FALSE

%%2、选取表达基因

%选取“表达”基因的方法常见的有两种，一是使用genefilter软件包，另外一种是调用affy包的mas5calls()函数。使用 genefilter需要设定筛选阈值，不同的人可能有不同的标准。mas5calls方法使用探针水平数据（AffyBatch类型数据）进行处理，一般使用没经过预处理的芯片数据通用性强些，其他参数用默认就可以。

data.mas5calls <-mas5calls(data.raw)%用mas5calls函数对文件数据进行处理

#Getting probe level data...

#Computing p-values

#Making P/M/A Calls

eset.mas5calls <-exprs(data.mas5calls)%继续用exprs计算“表达”量，得到的数据只有三个值P/M/A。

head(eset.mas5calls)%显示处理后的表达值

#            S1  S2  S3  S4  S5  S6  S7  S8

#1007_s_at "P" "P" "P" "P" "P" "P" "P" "P"

#1053_at   "P" "P" "P" "P" "P" "P" "P" "P"

#117_at    "A" "A" "A" "A" "A" "A" "A" "A"

#121_at    "P" "P" "P" "P" "P" "P" "P" "P"

#1255_g_at "A" "A" "A" "A" "A" "A" "A" "A"

#1294_at   "P" "P" "P" "P" "P" "M" "P" "P"

AP <-apply(eset.mas5calls, 1, function(x) any(x=="P"))%把至少一个芯片中有表达的基因选出来

present.probes <-names(AP[AP])%？

paste(length(present.probes),"/",length(AP))

#[1] "12757 / 22283"

rm(data.mas5calls)%删除data.mas5calls

rm(eset.mas5calls)%删除data.mas5calls

results.present <- results.rma[present.probes,]%present.probes是名称向量，用它进行数据子集提取。

%%3、选取差异表达基因

%生物学数据分析时的"差异"应该有两个意思，一是统计学上的差异，另外一个是生物学上的差异。一个基因在两个条件下的表达量分别有3个测量值：99,100,101 和 102,103,104。统计上两种条件下的基因表达数值是有差异的，后者比前者表达量要大。但生物学上有意义吗？未必。按平均值计算表达变化上升了3%，能产生什么样的生物学效应？这得看是什么基因了。所以差异表达基因的选取一般设置至少两个阈值：基因表达变化量和统计显著性量度（p值、q值等）。

1）t-测验

apply(abs(results.present[,5:7]), 2, max)%经常使用的筛选阈值是表达量变化超过2倍，即|log2(fc)|>=log(2)%上面语句的意思是对数据 abs(results.present[,5:7]) ，按列（第二维），使用统计函数max（计算最大值）

sum(abs(results.present[,"fc.7d"])>=log2(2))

results.st <- results.present[abs(results.present$fc.7d)>=log2(2),]%t测验

sel.genes <-row.names(results.st)

p.value <-apply(emat.rma.log2[sel.genes,], 1, function(x){t.test(x[1:2], x[7:8])$p.value})

results.st$p.value <- p.value

names(results.st)