gpr格式的芯片原始数据处理（神技能）

导读

我在生信技能树和生信菜鸟团都系统性总结过常见GEO数据库的芯片处理技巧，但是也有一些小众芯片我懒得去研发它的流程。但是我们生信技能树R语言之统计可视化讲师小洁却可以拿一整年的时间为学员们专心解决各式各样的疑难杂症，前面她已经分享了：aglient芯片原始数据处理，相信大家都获益良多，今天再来一个小众芯片：gpr格式的芯片原始数据处理

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)

image.png

plotMA(RG,array=6)

image.png

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的阈值，和常规的分析木有差别！后续常规相关分析流程大家可以在B站学习视频哈！

当然了，如果你计算机基础知识非常差，再多的资料给到你也是没有意义的，所以推荐我们的GEO数据挖掘课程哈，全年滚动开班，直播互动教学以及答疑，下一期是7月6号开课，可以花时间了解一下：