Affymetrix芯片分析：获取差异表达基因系列一

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

library(tcltk)%把tcltk包载入R中
filters <-matrix(c("CEL file", ".[Cc][Ee][Ll]", "All", "."), ncol = 2,byrow = T)%生成22字符矩阵，用来定义选择文件类型
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)%显示出所选文件的数量
sampleNames(data.raw)%查看文件名
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四个值依次重复两次所组成的数列
pData(data.raw)%指针pData函数读取文件

%1、计算基因表达量

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

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

%%用一个exprs函数就可以从eset数据中提取出表达量，得到的数据类型是矩阵。但是应该注意rma的eset结果是经过对数变换的，而mas5的eset结果是原始信号强度。虽然表达量是用对数变换的信号值表示的，但是有些计算过程要用到未经变换的原始值，应该把它们都计算出来：

emat.rma.log2 <-exprs(eset.rma)%rma的eset结果是经过对数变换的
emat.mas5.nologs <-exprs(eset.mas5)%mas5的eset结果是原始信号强度

class(emat.rma.log2)%结果是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结果

results.rma <-data.frame((emat.rma.log2[,c(1,3,5,7)] + emat.rma.log2[,c(2,4,6,8)])/2)%计算平均值，并转换为数据框格式#计算表达量差异倍数
results.rmafc.1h <\- results.rma[,2]-results.rma[,1] results.rmafc.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)%显示逻辑向量的长度
head(subset.logic)

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

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

eset.mas5calls <-exprs(data.mas5calls)%继续用exprs计算“表达”量，得到的数据只有三个值P/M/A。对于这三个值的具体解释可以用?mas5calls查看帮助。P为present，A为absent，M为marginal（临界值）。

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

AP <-apply(eset.mas5calls, 1, function(x) any(x=="P"))%把至少一个芯片中有表达的基因选出来
present.probes <-names(AP[AP])%？
paste(length(present.probes),"/",length(AP))
results.present <- results.rma[present.probes,]%present.probes是名称向量，用它进行数据子集提取。

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

3.1 简单t-测验
这种方法不用太多的统计学知识，生物专业的人很容易想到，而且确实有不少人在用。经常使用的筛选阈值是表达量变化超过2倍，即|log2(fc)|>=log(2)。先简单看看有没有：
apply(abs(results.present[,5:7]), 2, max)

fc.1h fc.24h fc.7d

5.309 6.688 6.844

apply是一个很有用的函数，它对数据按某个维度批量应用一个函数进行计算。第一个参数为向量或矩阵（或者是能转成向量或矩阵的数据，如数据框），第三个参数表示要使用的函数，第二个参数为应用的维度。上面语句的意思是对数据 abs(results.present[,5:7]) 按列（第二维）使用统计函数max（计算最大值）。表达变化超过2倍的基因共有842个：

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

[1] 842

results.st <- results.present[abs(results.present$fc.7d)>=log2(2),]
sel.genes <- row.names(results.st)

t测验，并选出p<0.05的差异表达基因：

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

[1] "S1" "S3" "S5" "S7" "fc.1h" "fc.24h" "fc.7d"

[8] "p.value"

results.st <- results.st[, c(1,4,7,8)]
results.st <- results.st[p.value<0.05,]
head(results.st, 2)
nrow(results.st)

3.2 SAM（Significance Analysis of Microarrays）
这种方法流行过一段时间，但由于FDR（错误检出率）控制太差，现在基本不用了。
要用也不复杂。但是注意SAM函数使用的emat表达数据是present.probes筛选出来的“表达”基因子集，如果你用没有经过筛选的数据，得到的结果会差别很大，不信可以自己试试（这点可能也是这种方法的毛病之一）。

library(BiocInstaller)
biocLite("samr")
library(samr)
samfit <- SAM(emat.rma.nologs[present.probes,c(1,2,7,8)], c(1,1,2,2), resp.type="Two class unpaired", genenames=present.probes)

SAM函数返回值一个列表结构，可以自己用?SAM看看。差异表达基因的数据在siggenes.table中，也是一个列表结构：
str(samfitsiggenes.table) 上调基因在siggenes.table的genes.up中，下调基因在genes.lo中。从上面的数据结构显示还可以看到差异表达基因的数量： ngenes.up和ngenes.lo。提取差异表达基因数据： results.sam <\- data.frame(rbind(samfitsiggenes.tablegenes.up,samfitsiggenes.table$genes.lo),
row.names=1, stringsAsFactors=FALSE)
for(i in 1:ncol(results.sam)) results.sam[,i] <- as.numeric(results.sam[,i])
head(results.sam, 2)

应用表达倍数进行筛选，有861个基因表达变化超过2倍（和前面简单t测验结果仅差1个，说明t测验还是可以的嘛！）：
results.sam <- results.sam[abs(log2(results.sam$Fold.Change))>=log2(2), ] ;

nrow(results.sam)

[1] 861

应用q值筛选，q<0.05只有10个，而q<0.1则有685个，选择筛选阈值也成了这种方法的一个问题：
#samr的q值表示方式为%，即5表示5%
nrow(results.sam[results.sam$q.val<5,])

[1] 10

nrow(results.sam[results.sam$q.val<10,])

[1] 685

3 Wilcoxon's signed-rank test
这个方法发表在 Liu, W.-m. et al, Analysis of high density expression microarrays with signed-rank call algorithms. Bioinformatics, 2002, 18, 1593-1599。R软件包simpleaffy的detection.p.val函数有实现，可以通过pairwise.comparison函数调用：

library(simpleaffy)
#注意下面语句中的数据顺序
sa.fit <- pairwise.comparison(eset.rma, "treatment", c("7d", "0h"))
pairwise.comparison返回的数据为simpleaffy自定义的"PairComp"类型，提取数据要用它专门的函数：平均值用means函数获得，变化倍数（log2）用fc函数获得，t测验的p值用tt函数获得：
class(sa.fit)

[1] "PairComp"

attr(,"package")

[1] "simpleaffy"

results.sa <- data.frame(means(sa.fit), fc(sa.fit), tt(sa.fit))
#选择有表达的基因

results.sa <- results.sa[present.probes,]
head(results.sa, 2)

X7d X0h fc.sa.fit. tt.sa.fit.

244901_at 4.047 4.203 -0.1562 0.43982

244902_at 3.938 4.295 -0.3570 0.05824

colnames(results.sa) <- c("7d", "0h", "fold.change", "p.val")
head(results.sa, 2)

7d 0h fold.change p.val

244901_at 4.047 4.203 -0.1562 0.43982

244902_at 3.938 4.295 -0.3570 0.05824

应用表达倍数筛选得到表达倍数超过2倍的基因数量有862个，应用p值筛选后得到562个差异表达基因：
results.sa <- results.sa[abs(results.sa$fold.change)>=log2(2), ]; nrow(results.sa)

[1] 862

results.sa <- results.sa[results.sa$p.val<0.05,]; nrow(results.sa)

[1] 562

4 Moderated T statistic
这种方法在R软件包limma里面实现得最好。limma最初主要用于双色（双通道）芯片的处理，现在不仅支持单色芯片处理，新版还添加了对RNAseq数据的支持，很值得学习使用。安装方法同前面其他Bioconductor软件包的安装。载入limm软件包后可以用limmaUsersGuide()函数获取pdf格式的帮助文档。limma需要先产生一个design矩阵，用于描述RNA样品：

library(limma)
treatment <- factor(pData(eset.rma)$treatment)
design <- model.matrix(~ 0 + treatment)
colnames(design) <- c("C0h", "T1h", "T24h", "T7d")
design

可以看到：矩阵的每一行代表一张芯片，每一列代表一种RNA来源（或处理）。此外，你可能还需要另外一个矩阵，用来说明你要进行哪些样品间的对比分析：

contrast.matrix <- makeContrasts(T1h-C0h, T24h-C0h, T7d-C0h, levels=design)

下一步建立线性模型，并进行分组比较和p值校正：
fit <- lmFit(eset.rma[present.probes,], design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)

nrow(topTable(fit2, coef=3, adjust.method="fdr", lfc=1, number=30000))

[1] 842

nrow(topTable(fit2, coef=3, adjust.method="fdr", p.value=0.05, lfc=1, number=30000))

[1] 740

把topTable函数的返回结果存到其他变量就可以了，它是数据框类型数据，可以用write或write.csv函数保存到文件：

results.lim <- topTable(fit2, coef=3, adjust.method="fdr", p.value=0.05, lfc=1, number=30000)
class(results.lim)

head(results.lim)

为什么以上几种方法仅用表达倍数（2倍）筛选得到的数字不大一样？limma和直接计算的结果都是842个，而simpleaffy和SAM为862/861个。这是对eset信号值取对数和求平均值的先后导致的，limma先取对数再求平均值，而simpleaffy和SAM是先求平均值再取对数。
5 其他方法：
如Rank products方法，在R软件包RankProd里实现，方法文献为：Breitling R, Armengaud P, Amtmann A, et al. Rank products: a simple, yet powerful, new method to detect differentially regulated genes in replicated microarray experiments[J]. FEBS letters, 2004, 573(1): 83-92.

#所有表达基因的名称
write(present.probes, "genes.expressed.txt")
#处理7天的差异表达基因
write.csv(results.lim, "results.lim.7d.csv")
#emat.rma.log2
write.csv(emat.rma.log2[present.probes,], "emat.rma.log2.csv")
如果要全部结果：
results.lim.all <- topTable(fit2, coef=1:3, adjust.method="fdr", p.value=1, lfc=0, number=30000)
head(results.lim.all, 3)
results.lim.all <- results.lim.all[, 1:3]
colnames(results.lim.all) <- c('T1h', 'T24h', 'T7d')
head(results.lim.all, 3)

T1h T24h T7d

254818_at 0.34085 6.024 6.215

245998_at -0.13675 3.676 2.778

265119_at -0.02536 6.061 4.380

write.csv(results.lim.all, 'results.lim.all.csv')