0.背景知識
參考地址：https://bioconductor.org/help/course-materials/2005/BioC2005/labs/lab01/beta7/
如果你在GEO數(shù)據(jù)庫挖過數(shù)據(jù)，有沒有碰到一些錯誤的、空的表達矩陣？有的作者上傳的表達矩陣是標準化過的，帶有負值，通常不拿這種標準化過的數(shù)據(jù)去做后續(xù)的分析。我們生信技能樹前面已經(jīng)分享過CEL格式的芯片原始數(shù)據(jù)處理方法。
最近復(fù)現(xiàn)文章時，發(fā)現(xiàn)了一些.gpr格式的x芯片原始數(shù)據(jù)，查了一下，發(fā)現(xiàn)是雙色芯片處理產(chǎn)生的文件，是用Genepix軟件得到的，比較古老的東西?？偨Y(jié)一下gpr格式的原始數(shù)據(jù)怎樣處理。
1.R包和文件準備
limma的userguide文檔里提到了gpr文件處理的代碼，沒有找到相應(yīng)的數(shù)據(jù)。我們使用的數(shù)據(jù)下載自：https://bioconductor.org/help/course-materials/2005/BioC2005/labs/lab01/Data/integrinbeta7.zip。
讀取gpr文件的R包不多，除了limma還有PAA、marray。limma夠用了，如果對另外兩個有興趣的話也可自己探索。
輸入數(shù)據(jù)是：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文件（可選）
這個芯片中設(shè)置了幾種除probe外的其他類型的位點,對應(yīng)的文件也在下載的文件夾里。使用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. 差異分析
試驗設(shè)計矩陣的構(gòu)建,有一點點差別，是直接從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
有了試驗設(shè)計矩陣，差異分析就很簡單，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)
可以統(tǒng)計一下差異基因的數(shù)量
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
差異基因數(shù)量個位數(shù)，可以考慮調(diào)整logFC的閾值，和常規(guī)的分析木有差別
