芯片表达谱分析（二）

本文取 GSE45636 数据演示单色荧光芯片数据分析。GSE45636 采用 Affymetrix Human Genome U133 Plus 2.0 Array 平台在 Bioconductor - 3.13 AnnotationData Packages 网站可查到其芯片注释包为 hgu133plus2.db.

第一步导入需要的 R 包，hgu133plus2.db 是芯片注释包，oligo 是处理 Affymetrix 芯片数据的包，limma 是进行差异基因分析的包，而 tidyverse 方便数据框操作，以及提供管道符。

library(hgu133plus2.db)
library(oligo)
library(limma)
library(tidyverse)

准备好需要读取的数据路径列表。用 list.celfiles 函数自动获取目录下的 CEL 文件，就不需要手打文件名了，避免错误。

wd <- file.path("/Example/GSE45636")
sampleName <- paste("GSM", 1111159:1111164, sep = "")
celFiles <- list.celfiles(wd, listGzipped = TRUE, full.names = TRUE)

> sampleName
[1] "GSM1111159" "GSM1111160" "GSM1111161" "GSM1111162" "GSM1111163"
[6] "GSM1111164"

> celFiles
[1] "/Example/GSE45636/GSM1111159_1681_035_HT31_1-_h01_AK_050712.CEL.gz"
[2] "/Example/GSE45636/GSM1111160_1682_035_HT31_1+_h01_AK_050712.CEL.gz"
[3] "/Example/GSE45636/GSM1111161_1683_035_HT31_2-_h01_AK_050712.CEL.gz"
[4] "/Example/GSE45636/GSM1111162_1684_035_HT31_2+_h01_AK_050712.CEL.gz"
[5] "/Example/GSE45636/GSM1111163_1685_035_HT31_3-_h01_AK_050712.CEL.gz"
[6] "/Example/GSE45636/GSM1111164_1686_035_HT31_3+_h01_AK_050712.CEL.gz"

用 read.celfiles 函数读取原始数据（CEL 格式文件）成 ExpressionSet 对象，然后用 oligo 的 boxplot 函数查看表达值的箱线图。如果准备有 phenoData 和 featureData 数据，可以此时用相应参数整合至 ExpressionSet 对象。
ExpressionSet 详细解释在《ExpressionSet 数据结构》文章。

rawData <- read.celfiles(filenames = celFiles)
oligo::boxplot(rawData)

不同样品信号强度存在差异，不适合进行比较

读取以后要进行背景矫正、样品间矫正、总聚（summarization）等，否则因为技术误差，表达值的比较将是不准确的。总聚是因为芯片设计转录本是由多个探针构成的。使用 oligo 包，这些步骤都被包装到 rma 函数，使用 RMA(Robust Multichip Average) 算法处理。要注意的是 RMA 处理后数据为 log2 后的，后续可以直接用于 limma 包，不要再手动去 log2 处理。

rmaData <- rma(rawData)
oligo::boxplot(rmaData)

样品总体的数值分布已经类似，中位数基本一致

然后可以移除一些低表达基因，因为芯片数据许多探针信号会处在背景信号范围，那些探针不仅低表达，且往往不是差异表达的。
生成每个探针的表达中位值频率图。

probeMedian <- rowMeans(exprs(rmaData))
hist(probeMedian, 100)

探针信号频数

在这里，可以移除左边的一部分探针数据，阈值根据图像自己选择，我这里选择 4.8 为界。这里一共 6 样本，我过滤设置至少 2 样本表达值不低于 4.8 的探针保留。过滤阈值的选择是自由的，要记得过滤的目的是为了移除一些不表达或低表达探针。

keepProbe <- rowSums(exprs(rmaData) >= 4.8) >= 2
rmaData2 <- subset(rmaData, keepProbe)

下一步是注释，也就是添加 featureData 信息，因为我们在导入数据时没有添加。
先查看注释包有哪些注释信息，自己需要哪些信息。
注释 R 包的使用见《AnnotationDbi 使用》文章。

> keytypes(hgu133plus2.db)
 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS" "ENTREZID"    
 [7] "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
[13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"     "ONTOLOGYALL"  "PATH"        
[19] "PFAM"         "PMID"         "PROBEID"      "PROSITE"      "REFSEQ"       "SYMBOL"      
[25] "UCSCKG"       "UNIGENE"      "UNIPROT" 

我这里只取探针 ID 和基因 SYMBOL 的注释。注意明确指定用 AnnotationDbi 包的 select 函数。可能存在命名冲突的函数，调用时指定 R 包来源是好习惯。

probes <- featureNames(rmaData2)
annotData <- AnnotationDbi::select(x = hgu133plus2.db, keys = probes, 
                                   columns = c("PROBEID", "SYMBOL"), keytype = "PROBEID") %>% 
  dplyr::filter(!is.na(SYMBOL))

检查一下是否有探针对应到多个基因的情况。

> count(annotData, PROBEID) %>% dplyr::filter(n > 1) %>% glimpse()
Rows: 2,031
Columns: 2
$ PROBEID <chr> "1007_s_at", "1294_at", "1552283_s_at", "1552401_a_at", "1552411_at", "1552412_a_a…
$ n       <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 4, 5, 2, 2, 2, 2, 2, 3, 2, 3, 2, 2, 2, 2, 2, 3, 2, 2, 2…

可见有 2031 个探针存在对应着多个基因的情况。取几个看看情况。

> dplyr::filter(annotData, PROBEID %in% c("1007_s_at", "1552401_a_at", "1552411_at"))
       PROBEID    SYMBOL
1    1007_s_at      DDR1
2    1007_s_at   MIR4640
3 1552401_a_at     BACH1
4 1552401_a_at GRIK1-AS2
5   1552411_at  DEFB106A
6   1552411_at  DEFB106B

这种情况，是无法区分对应到哪个基因是 “正确” 的，因此选择直接移除这些探针注释。

fData(rmaData2)$PROBEID <- probes
fData(rmaData2) <- fData(rmaData2) %>% left_join(annotData2, by="PROBEID")

此时，表达数据已经整理好，经 RMA 算法处理已经是 log2 后的，可以直接用 limma 包进行差异基因分析。
先准备分组信息和 design model 等。在这里，建议学习一下 R 语言的线性模型和公式（formula）等知识点。

sampleGroup <- factor(c(rep_len(c("CTR", "KD"), 3)))
dgm <- model.matrix(~ 0 + sampleGroup)
rownames(dgm) <- sampleName
colnames(dgm) <- c("CTR", "KD")
contrast <- makeContrasts(contrasts = c("KD-CTR"), levels = dgm)

结果如下。设计矩阵 dgm 行为样品，列为线性模型变量。

> sampleGroup
[1] CTR KD  CTR KD  CTR KD 
Levels: CTR KD
> dgm
           CTR KD
GSM1111159   1  0
GSM1111160   0  1
GSM1111161   1  0
GSM1111162   0  1
GSM1111163   1  0
GSM1111164   0  1
attr(,"assign")
[1] 1 1
attr(,"contrasts")
attr(,"contrasts")$sampleGroup
[1] "contr.treatment"

> contrast
      Contrasts
Levels KD-CTR
   CTR     -1
   KD       1

用 limma 完成差异基因分析。因为样品名之前读取数据时自动获取的，跟设计矩阵的行名不同，这里先进行统一。然后用 lmFit 应用线性模型，而 topTable 函数取得差异基因数据，它的 number 参数设置提取差异基因数据的数目，用 Inf 表示提取所有数据，无论是否差异。方便后面自己根据需要去筛选差异基因。

sampleNames(rmaData2) <- sampleName
fit1 <- lmFit(rmaData2, design = dgm)
fit2 <- contrasts.fit(fit1, contrasts = contrast)
fit3 <- eBayes(fit2)
allGene <- topTable(fit3, coef = 1, number = Inf)

看一下结果格式，可见之前添加的 featureData 会添加到结果里，输出的结果可以自己根据差异倍数和矫正后 P 值进行筛选。

> head(allGene)
                  PROBEID  SYMBOL     logFC   AveExpr         t      P.Value    adj.P.Val        B
1557129_a_at 1557129_a_at FAM111B -3.486052  7.386077 -31.15664 3.569462e-10 4.310448e-06 13.50997
205436_s_at   205436_s_at    H2AX -2.934340  9.705511 -29.77089 5.269543e-10 4.310448e-06 13.21950
206157_at       206157_at    PTX3  4.725968  7.463342  29.71943 5.348144e-10 4.310448e-06 13.20829
225956_at       225956_at  CREBRF  3.474989  8.444295  29.23024 6.164295e-10 4.310448e-06 13.10014
206461_x_at   206461_x_at    MT1H  2.567942 11.844680  29.17308 6.268389e-10 4.310448e-06 13.08731
36711_at         36711_at    MAFF  3.124939  9.646905  28.91136 6.770750e-10 4.310448e-06 13.02806

用火山图展示：

火山图

参考资料
limma: Linear Models for Microarray and RNA-Seq Data User’s Guide