R实例：limma包分析Agilent4x44KArrays

目前网页上常见的是使用Affy的包对Affymetrix的芯片分析教程，而忽略了对其他平台芯片分析。limma包是一个针对芯片数据进行分析的R包，支持不同平台的芯片数据。

数据获取：http://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-33005.
研究内容：Effects of corn oil gavage on gene expression profiles of thymus in Wistar Rats

Download.png

下载上述两个文件，其中sample and data relationship format (SDRF) 文件之后的用处等同于limma包中的targets文件（targets文件见下图）

Targets文件(in Userguide).png

#读取SDRF文件
SDRF <- read.delim("E-GEOD-33005.sdrf.txt",check.names=FALSE,stringsAsFactors=FALSE)
#查看SDRF中有哪些属性，其中Array Data File是我们需要的
names(SDRF)
[1] "Source Name"
[2] "Comment [Sample_description]"
[3] "Comment [Sample_source_name]"
[4] "Characteristics[Organism]"
[5] "Term Source REF"
[6] "Term Accession Number"
[7] "Characteristics[age]"
[8] "Characteristics[gender]"
[9] "Term Source REF"
[10] "Term Accession Number"
[11] "Characteristics[strain]"
[12] "Term Source REF"
[13] "Term Accession Number"
[14] "Characteristics[tissue]"
[15] "Term Source REF"
[16] "Term Accession Number"
[17] "Characteristics[treatment]"
[18] "Protocol REF"
[19] "Term Source REF"
[20] "Protocol REF"
[21] "Term Source REF"
[22] "Protocol REF"
[23] "Term Source REF"
[24] "Extract Name"
[25] "Material Type"
[26] "Protocol REF"
[27] "Term Source REF"
[28] "Labeled Extract Name"
[29] "Label"
[30] "Protocol REF"
[31] "Term Source REF"
[32] "Hybridization Name"
[33] "Array Design REF"
[34] "Term Source REF"
[35] "Protocol REF"
[36] "Term Source REF"
[37] "Protocol REF"
[38] "Term Source REF"
[39] "Protocol REF"
[40] "Term Source REF"
[41] "Array Data File"
[42] "Comment [ArrayExpress FTP file]"
[43] "Protocol REF"
[44] "Term Source REF"
[45] "Normalization Name"
[46] "Derived Array Data File"
[47] "Comment [Derived ArrayExpress FTP file]"
[48] "FactorValue [TREATMENT]"
setwd("~/Downloads/agilent/E-GEOD-33005.raw.1")
#read.maimages中的source是指检测芯片荧光值所用的仪器平台而非芯片平台，可以在https://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-33005/protocols/连接中查看
x <- read.maimages(SDRF[,"Array Data File"],source="agilent",green.Only=TRUE)
Read GSM819076_US10283824_252828210181_S01_GE1_107_Sep09_1_4.txt
Read GSM819075_US10283824_252828210181_S01_GE1_107_Sep09_1_3.txt
Read GSM819074_US10283824_252828210181_S01_GE1_107_Sep09_1_2.txt
Read GSM819073_US10283824_252828210180_S01_GE1_107_Sep09_1_4.txt
Read GSM819072_US10283824_252828210180_S01_GE1_107_Sep09_1_3.txt
Read GSM819071_US10283824_252828210180_S01_GE1_107_Sep09_1_2.txt
Read GSM819070_US10283824_252828210180_S01_GE1_107_Sep09_1_1.txt
Read GSM819069_US10283824_252828210179_S01_GE1_107_Sep09_1_4.txt
Read GSM819068_US10283824_252828210179_S01_GE1_107_Sep09_1_3.txt
Read GSM819067_US10283824_252828210179_S01_GE1_107_Sep09_1_2.txt
Read GSM819066_US10283824_252828210179_S01_GE1_107_Sep09_1_1.txt
Read GSM819065_US10283824_252828210178_S01_GE1_107_Sep09_1_4.txt
Read GSM819064_US10283824_252828210178_S01_GE1_107_Sep09_1_3.txt
Read GSM819063_US10283824_252828210178_S01_GE1_107_Sep09_1_2.txt
Read GSM819062_US10283824_252828210178_S01_GE1_107_Sep09_1_1.txt
Read GSM819061_US10283824_252828210177_S01_GE1_107_Sep09_1_4.txt
Read GSM819060_US10283824_252828210177_S01_GE1_107_Sep09_1_3.txt
Read GSM819059_US10283824_252828210177_S01_GE1_107_Sep09_1_2.txt
Read GSM819058_US10283824_252828210177_S01_GE1_107_Sep09_1_1.txt


Source.png

#进行背景校正
y <- backgroundCorrect(x,method="normexp")
Array 1 corrected
Array 2 corrected
Array 3 corrected
Array 4 corrected
Array 5 corrected
Array 6 corrected
Array 7 corrected
Array 8 corrected
Array 9 corrected
Array 10 corrected
Array 11 corrected
Array 12 corrected
Array 13 corrected
Array 14 corrected
Array 15 corrected
Array 16 corrected
Array 17 corrected
Array 18 corrected
Array 19 corrected
#进行归一化
y <- normalizeBetweenArrays(y,method="quantile")
#控制探针质量去除低表达的探针：计算每张芯片中95%的 negative control probes，保留在至少四张芯片中比 negative control probes强度高10%以上的探针
neg95 <- apply(y$E[y$genes$COntrolType==-1,],2,function(x) quantile(x,p=0.95))
cutoff <- matrix(1.1*neg95,nrow(y),ncol(y),byrow=TRUE)
isexpr <- rowSums(y$E > cutoff) >= 4
table(isexpr)
isexpr
FALSE TRUE
11500 32754
#Regular probes are code as “0” in the ControlType column
y0 <- y[y$genes$COntrolType==0 & isexpr,]


对芯片进行完前期处理后，进一步进行基因差异的分析：

#获取不同芯片的实验处理
Treatment <- SDRF[,"Characteristics[treatment]"]
levels <- c("10 ml/kg saline","2 ml/kg corn oil","5 ml/kg corn oil","10 ml/kg corn oil")
Treatment <- factor(Treatment,levels=levels)
#制作treatment 的矩阵
design <- model.matrix(~Treatment)
fit <- lmFit(y0,design)
#贝叶斯校正
fit <- eBayes(fit,trend=TRUE)
plotSA(fit, main="Probe-level")
summary(decideTests(fit[,-1]))
Treatment2 ml/kg corn oil Treatment5 ml/kg corn oil Treatment10 ml/kg corn oil
Down 0 0 911
NotSig 32723 32723 30063
Up 0 0 1749


plotSA-Probe.png

根据上述分析发现，only the 10 ml/kg treatment 与the saline control
有所不同。

同样的分析，但是是针对每个基因而不是探针：

yave <- avereps(y0,ID=y0$genes[,"SystematicName"])
fit <- lmFit(yave,design)
fit <- eBayes(fit,trend=TRUE)
plotSA(fit, main="Gene-level")
summary(decideTests(fit[,-1]))
Treatment2 ml/kg corn oil Treatment5 ml/kg corn oil Treatment10 ml/kg corn oil
-1 0 0 407
0 19112 19112 17878
1 0 0 827


plotSA-Gene.png

进一步导出比较差异的数据：

tempOutput = topTable(fit, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG,"limma_notrend.results.csv",quote = F)
head(nrDEG)
Row Col ControlType ProbeName SystematicName logFC AveExpr t
NM_021264 177 40 0 A_64_P085690 NM_021264 17.35790 17.36775 519.4180
NM_212510 84 56 0 A_64_P074077 NM_212510 17.00460 17.00126 488.1011
NM_012876 187 12 0 A_64_P049933 NM_012876 17.12797 17.11662 486.5568
ENSRNOT00000030778 87 50 0 A_64_P296417 ENSRNOT00000030778 16.98037 16.97319 477.2043
XM_230930 363 34 0 A_44_P103028 XM_230930 17.09420 17.05256 468.9523
NM_001009239 92 25 0 A_64_P057441 NM_001009239 16.66621 16.68367 466.1435
P.Value adj.P.Val B
NM_021264 1.534832e-39 2.826651e-35 63.47183
NM_212510 4.779014e-39 2.826651e-35 63.32308
NM_012876 5.063788e-39 2.826651e-35 63.31506
ENSRNOT00000030778 7.218211e-39 2.826651e-35 63.26496
XM_230930 9.926515e-39 2.826651e-35 63.21847
NM_001009239 1.107768e-38 2.826651e-35 63.20214


保存的CSV文件如下：

geneexpression.png