source("http://bioconductor.org/biocLite.R")

biocLite("ArrayExpress")

library(ArrayExpress)

AEset<-ArrayExpress("E-GEOD-470")%GEO中哮喘相关的基因表达谱数据（编号GSE470）

save(AEset,file="AEset.RData")

load("AEset.RData")

class(AEset)%三种查看方式

AEset

pData(AEset)

colnames(pData(AEset))

fac<-colnames(pData(AEset))[grep("Factor",colnames(pData(AEset)))]%查找colnames(pData(AEset))带中字符“Facotor”

biocLite("arrayQualityMetrics")

library(arrayQualityMetrics)

arrayQualityMetrics(expressionset=AEset,outdir="QAraw",force=FALSE,do.logtransform=TRUE,intgroup=fac)

library(affy)

rAEset<-rma(AEset)

arrayQualityMetrics(expressionset=rAEset,outdir="QAnorm",force=FALSE,do.logtransform=TRUE,intgroup=fac)

cAEset<-rma(AEset[,-2])%删除第二列

%根据假设检验的原理，进行基因表达差异分析

llibrary(siggenes)

library(multtest)

sc<-c(0,1,0,1,0,0,1,1,0,1,1)

sub<-exprs(cAEset[,c(1:6,7:11)])

gn<-geneNames(AEset)

sam.out<-sam(sub,sc,rand=123,gene.names=gn);

sam.out

plot(sam.out)%通过plot显示差异表达的图形% Obtain the SAM plot for Delta = 2>plot(sam.out, 2)

sum.sam.out<-summary(sam.out,3)% Get information about the genes called significant using Delta = 3.

sum.sam.out@row.sig.genes%Obtain the rows of golub containing the genes called differentially expressed

sum.sam.out@mat.sig%obtain the matrix containing the d-values, q-values etc. of the differentially expressed genes

sam(sub, sc, method="wilc.stat", rand=123)

list.siggenes(sam.out,3)%for Delta = 3

plot(sam.out,3)

biocLite("estrogen")

library(estrogen)

biocLite("genefilter")

library(genefilter)

heatmap(sub[siggn,],col=gentlecol(256))%对于检测到的差异表达基因，用heatmap图来直观显示其表达情况

library(vsn)

rsd<-rowSds(exprs(cAEset))%对基因表达数据进行标准分析，将标准差进行降序排列，选取前25位的基因进行双向聚类，得到heatmap

sel<-order(rsd,decreasing=TRUE)[1:50]

heatmap(exprs(cAEset)[sel,],col=gentlecol(256))