R实例：limma包分析Agilent 4x44K Arrays

目前网页上常见的是使用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