0.背景知识
参考地址：https://bioconductor.org/help/course-materials/2005/BioC2005/labs/lab01/beta7/
如果你在GEO数据库挖过数据，有没有碰到一些错误的、空的表达矩阵？有的作者上传的表达矩阵是标准化过的，带有负值，通常不拿这种标准化过的数据去做后续的分析。我们生信技能树前面已经分享过CEL格式的芯片原始数据处理方法。
最近复现文章时，发现了一些.gpr格式的x芯片原始数据，查了一下，发现是双色芯片处理产生的文件，是用Genepix软件得到的，比较古老的东西。总结一下gpr格式的原始数据怎样处理。
1.R包和文件准备
limma的userguide文档里提到了gpr文件处理的代码，没有找到相应的数据。我们使用的数据下载自：https://bioconductor.org/help/course-materials/2005/BioC2005/labs/lab01/Data/integrinbeta7.zip。
读取gpr文件的R包不多，除了limma还有PAA、marray。limma够用了，如果对另外两个有兴趣的话也可自己探索。
输入数据是：6个gpr文件，一个target文件。
rm(list = ls())
library(limma)
dir()
## [1] "6Hs.166.gpr" "6Hs.168.gpr" "6Hs.187.1.gpr" "6Hs.194.gpr"
## [5] "6Hs.195.1.gpr" "6Hs.243.1.gpr" "gpr.html" "gpr.Rmd"
## [9] "gpr_example.R" "gprdata.R" "main.Rproj" "SpotTypes.txt"
## [13] "TargetBeta7.txt"
targets <- rio::import("TargetBeta7.txt")
targets
## FileNames Subject ID # Cy3 Cy5 Hyb buffer Hyb Temp (deg C)
## 1 6Hs.195.1.gpr 1 b7 - b7 + Ambion Hyb Slide 55
## 2 6Hs.168.gpr 3 b7 + b7 - Ambion Hyb Slide 55
## 3 6Hs.166.gpr 4 b7 + b7 - Ambion Hyb Slide 55
## 4 6Hs.187.1.gpr 6 b7 - b7 + Ambion Hyb Slide 55
## 5 6Hs.194.gpr 8 b7 - b7 + Ambion Hyb Slide 55
## 6 6Hs.243.1.gpr 11 b7 + b7 - Ambion Hyb Slide 55
## Hyb Time (h) Date of Blood Draw Amplification Slide Type Date of Scan
## 1 40 2002.10.11 R2 aRNA Aminosilane 2003.07.25
## 2 40 2003.01.16 R2 aRNA Aminosilane 2003.08.07
## 3 40 2003.01.16 R2 aRNA Aminosilane 2003.08.07
## 4 40 2002.09.16 R2 aRNA Aminosilane 2003.07.18
## 5 40 2002.09.18 R2 aRNA Aminosilane 2003.07.25
## 6 40 2003.01.13 R2 aRNA Aminosilane 2003.08.06
f <- function(x) as.numeric(x$Flags > -75)
RG <- read.maimages(targets, source="genepix", wt.fun=f)
## Read 6Hs.195.1.gpr
## Read 6Hs.168.gpr
## Read 6Hs.166.gpr
## Read 6Hs.187.1.gpr
## Read 6Hs.194.gpr
## Read 6Hs.243.1.gpr
target文件在limma帮助文档中也有说明，有三列必须的：一列是gpr文件名，另两列是Cy3和Cy5,表示在每个gpr文件中Cy3和Cy5两种染料标记了哪组样本。 ### 2.读取spottypes文件（可选）
这个芯片中设置了几种除probe外的其他类型的位点,对应的文件也在下载的文件夹里。使用plotMA可以查看每个gpr文件中spot类型的分布。
spottypes <- readSpotTypes()
spottypes
## SpotType Name ID col cex
## 1 Probe * * black 0.2
## 2 Empty Empty|EMPTY Empty|EMPTY yellow 1.0
## 3 Buffer Buffer* Buffer orange 1.0
## 4 Negative Randomized negative control H2NC* magenta 1.0
## 5 Arabidopsis *A. thaliana* AT00* green 1.0
## 6 SPC *Stratagene Positive Control AT00* blue 1.0
## 7 BetaActin ACTB - Actin, beta H2* red 1.0
RG$genes$Status <- controlStatus(spottypes, RG)
## Matching patterns for: Name ID
## Found 23184 Probe
## Found 1328 Empty
## Found 3 Buffer
## Found 192 Negative
## Found 30 Arabidopsis
## Found 3 SPC
## Found 34 BetaActin
## Setting attributes: values col cex
plotMA(RG)
plotMA(RG,array=6)
3.背景校正和标准化
RGne <- backgroundCorrect(RG, method="normexp", offset=25)
## Array 1 corrected
## Array 2 corrected
## Array 3 corrected
## Array 4 corrected
## Array 5 corrected
## Array 6 corrected
## Array 1 corrected
## Array 2 corrected
## Array 3 corrected
## Array 4 corrected
## Array 5 corrected
## Array 6 corrected
MA <- normalizeWithinArrays(RGne)
此时得到的MA文件就类似于表达矩阵的前体了，可以直接拿来做差异分析
4. 差异分析
试验设计矩阵的构建,有一点点差别，是直接从targets里面得到的分组。
design <- modelMatrix(targets, ref="b7 -")
## Found unique target names:
## b7 - b7 +
design
## b7 +
## 1 1
## 2 -1
## 3 -1
## 4 1
## 5 1
## 6 -1
design2 <- cbind(Dye=1, design)
colnames(design2)[2] = "Beta7"
design2
## Dye Beta7
## 1 1 1
## 2 1 -1
## 3 1 -1
## 4 1 1
## 5 1 1
## 6 1 -1
有了试验设计矩阵，差异分析就很简单，lmFit、eBayes、topTable三部搞定：
fit <- lmFit(MA, design2)
## Warning: Partial NA coefficients for 3 probe(s)
fit <- eBayes(fit)
deg <- topTable(fit, coef="Beta7", adjust="fdr",number = Inf)
deg = deg[,c((ncol(deg)-5):ncol(deg),1:(ncol(deg)-6))]
deg = na.omit(deg)
可以统计一下差异基因的数量
deg$change = ifelse(deg$logFC > 1 & deg$adj.P.Val<0.05,
"up",
ifelse(deg$logFC < -1 & deg$adj.P.Val<0.05,
"down",
"stable"))
table(deg$change)
##
## down stable up
## 4 21844 3
差异基因数量个位数，可以考虑调整logFC的阈值，和常规的分析木有差别
