目前網(wǎng)頁上常見的是使用Affy的包對Affymetrix的芯片分析教程台舱，而忽略了對其他平臺芯片分析绒北。limma包是一個針對芯片數(shù)據(jù)進行分析的R包话侧，支持不同平臺的芯片數(shù)據(jù)筹陵。
數(shù)據(jù)獲取：http://www.ebi.ac.uk/arrayexpress/experiments/E-GEOD-33005.
研究內(nèi)容：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")
#控制探針質(zhì)量去除低表達的探針：計算每張芯片中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
根據(jù)上述分析發(fā)現(xiàn)，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
進一步導(dǎo)出比較差異的數(shù)據(jù)：
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
