本文演示双色荧光芯片数据分析。选取 GSE9187 数据集 GSM232034-GSM232039 样品作为案例，包含对照组和敲低组各 3 样品。
导入需要的包并且准备好文件路径。
library(tidyverse)
library(limma)
sampleName <- c("GSM232034", "GSM232035", "GSM232036", "GSM232037",
"GSM232038", "GSM232039")
filePath <- file.path(paste(sampleName, ".txt.gz", sep = ""))
names(filePath) <- sampleName
> head(filePath, n = 2)
GSM232034 GSM232035
"GSM232034.txt.gz" "GSM232035.txt.gz"
用 limma 包的 read.maimages 函数读取文件。
rg1 <- read.maimages(files = filePath, source = "agilent", names = sampleName)
根据不同的芯片设置不同的 source 参数，总结如下表。
source 参数
读取后数据结构如下，其中 G 是绿色荧光，R 是红色荧光。
> head(rg1, n = 2)
An object of class "RGList"
$G
GSM232034 GSM232035 GSM232036 GSM232037 GSM232038 GSM232039
[1,] 17008 14792 13002 20347.5 13638.5 17413.5
[2,] 24 31 34 20.0 31.0 32.5
$Gb
GSM232034 GSM232035 GSM232036 GSM232037 GSM232038 GSM232039
[1,] 19 18 17 19 18 18
[2,] 19 19 18 19 18 18
$R
GSM232034 GSM232035 GSM232036 GSM232037 GSM232038 GSM232039
[1,] 56639 36160.5 30719 40252.5 18309.5 31936
[2,] 20 24.0 25 25.0 24.0 24
$Rb
GSM232034 GSM232035 GSM232036 GSM232037 GSM232038 GSM232039
[1,] 21 20 21 20 20 20
[2,] 20 20 21 20 20 20
$targets
FileName
GSM232034 GSM232034.txt.gz
GSM232035 GSM232035.txt.gz
GSM232036 GSM232036.txt.gz
GSM232037 GSM232037.txt.gz
GSM232038 GSM232038.txt.gz
GSM232039 GSM232039.txt.gz
$genes
Row Col Start Sequence ProbeUID ControlType ProbeName
1 1 1 0 0 1 GE_BrightCorner
2 1 2 0 1 1 DarkCorner
GeneName SystematicName Description
1 GE_BrightCorner GE_BrightCorner
2 DarkCorner DarkCorner
$source
[1] "agilent"
对每个样品生成 MA 图并检查。MA 图在芯片和 RNA-Seq 分析都常用，展示总体的表达均值和差异倍数关系。一般来说分析的模型都会假定不同组总体表达水平不变，大部分基因表达无差异，少部分基因表达有差异。
在双色芯片，假设 表示 i 基因的绿色荧光强度，
表示 i 基因的红色荧光强度，那么一般 MA 值计算公式为：
对于 "dye-swap" 实验设计，假设第一张芯片 Treatment/Control 为 Green/Red 那么第二张将为 Red/Green, 两张芯片计算 M 值的分子分母颜色相反，最终用两芯片 M 值平均值。
> plotMA(rg1, array = 1, main = sampleName[1])
MA 图
进行背景信号矫正，使用 offset 参数给所有信号值加上一个常数，从而低信号值在 log 转换后的 log 比值将趋近于零。方法 "normexp" 处理后将使所有矫正后信号为正值。
If method="normexp" a convolution of normal and exponential distributions is fitted to the fore-ground intensities using the background intensities as a covariate, and the expected signal given the observed foreground becomes the corrected intensity. This results in a smooth monotonic trans-formation of the background subtracted intensities such that all the corrected intensities are positive.
rg2 <- backgroundCorrect(rg1, method = "normexp", offset = 50)
进行芯片内标准化（Within-Array Normalization），这将调整芯片内部信号的技术误差。这步将计算 log2 后的 MA 值。安捷伦的芯片适合 "loess" 方法。
rg3 <- normalizeWithinArrays(rg2, method = "loess")
> head(rg3, n = 2)
An object of class "MAList"
$targets
FileName
GSM232034 GSM232034.txt.gz
GSM232035 GSM232035.txt.gz
GSM232036 GSM232036.txt.gz
GSM232037 GSM232037.txt.gz
GSM232038 GSM232038.txt.gz
GSM232039 GSM232039.txt.gz
$genes
Row Col Start Sequence ProbeUID ControlType ProbeName
1 1 1 0 0 1 GE_BrightCorner
2 1 2 0 1 1 DarkCorner
GeneName SystematicName Description
1 GE_BrightCorner GE_BrightCorner
2 DarkCorner DarkCorner
$source
[1] "agilent"
$M
GSM232034 GSM232035 GSM232036 GSM232037 GSM232038
[1,] -0.33832995 -0.6505946 -0.7186429 -0.4015876 -1.04712940
[2,] -0.09966227 0.0200928 -0.1233429 0.1504817 0.01073333
GSM232039
[1,] -0.6033589
[2,] -0.0532774
$A
GSM232034 GSM232035 GSM232036 GSM232037 GSM232038 GSM232039
[1,] 14.923482 14.500260 14.289791 14.806421 13.953195 14.527340
[2,] 5.786474 6.266751 6.106519 5.794016 6.481864 5.885135
再度检查 MA 图。
> plotMA(rg3, array = 1, main = sampleName[1])
MA 图
进行芯片间标准化，这一步要注意不能有 0 或者负数的值输入，在背景矫正时我们用 "normexp" 方法已经确保了数值为正。
rg4 <- normalizeBetweenArrays(rg3, method = "Aquantile")
比较标准化前后的信号密度图。
> plotDensities(rg3)
> plotDensities(rg4)
标准化前密度图
标准化后密度图
移除后续分析不需要的 control 探针，探针类型由 ControlType 值表示，下面表格总结了一些探针类型的值。其中阳性对照一般采用物种特异性的内源序列，保证有强烈的信号；而阴性对照用于衡量信号背景水平。
| 值 | 含义 |
|---|---|
| 0 | Control type none |
| 1 | Positive control |
| -1 | Negative control |
| -15000 | SNP |
| -20000 | Not probe |
| -30000 | Ignore |
先查看目前数据包含哪些探针类型。
> rg4$genes$ControlType %>% unique()
[1] 1 0 -1
移除 control 探针。
keepProbe <- rg4$genes$ControlType == 0
rg5 <- rg4[keepProbe,]
芯片一般会不少探针是分布在多个 spot 的，因此会有多条记录。
> rg5$genes$ProbeName %>% table() %>% sort(decreasing = TRUE) %>% head(n = 3)
.
A_23_P102000 A_23_P102471 A_23_P103361
10 10 10
用 avereps 函数为这些探针取平均值。
rg6 <- avereps(rg5, ID = rg5$genes$ProbeName)
之后进行差异基因（探针）分析，先构建实验设计矩阵，从 GEO 上可以看到样本都使用人种标准参考 RNA 并用绿色荧光（Cy3）标记。
将标记荧光信息添加到 targets 数据框。
rg6$targets$Cy3 <- "Ref"
rg6$targets$Cy5 <- c("CTR", "CTR", "CTR", "KD", "KD", "KD")
> rg6$targets
FileName Cy3 Cy5
GSM232034 GSM232034.txt.gz Ref CTR
GSM232035 GSM232035.txt.gz Ref CTR
GSM232036 GSM232036.txt.gz Ref CTR
GSM232037 GSM232037.txt.gz Ref KD
GSM232038 GSM232038.txt.gz Ref KD
GSM232039 GSM232039.txt.gz Ref KD
然后使用 modelMatrix 函数，其中 ref 参数指定参考的 RNA 样品，在这里是 Cy3 的人种参考 RNA.
rgDesign <- modelMatrix(rg6$targets, ref = "Ref")
> rgDesign
CTR KD
GSM232034 1 0
GSM232035 1 0
GSM232036 1 0
GSM232037 0 1
GSM232038 0 1
GSM232039 0 1
用 limma 分析得到差异的探针（基因）。
fit1 <- lmFit(rg6, design = rgDesign)
contrast <- makeContrasts(KD-CTR, levels = rgDesign)
fit2 <- contrasts.fit(fit1, contrast)
fit3 <- eBayes(fit2)
allGene <- topTable(fit3, coef = 1, number = Inf, genelist = rg6$genes)
结果概览
> glimpse(allGene, width = 50)
Rows: 41,000
Columns: 16
$ Row <int> 41, 50, 35, 257, 153, 370…
$ Col <int> 28, 24, 70, 28, 41, 47, 5…
$ Start <int> 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ Sequence <chr> "TTTCCTTCTGATGCAGACATGGTC…
$ ProbeUID <int> 3055, 3696, 2650, 19943, …
$ ControlType <int> 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ ProbeName <chr> "A_32_P72351", "A_23_P120…
$ GeneName <chr> "AK026140", "HMOX1", "KLR…
$ SystematicName <chr> "AK026140", "NM_002133", …
$ Description <chr> "Homo sapiens cDNA: FLJ22…
$ logFC <dbl> 3.564953, -2.431149, 2.33…
$ AveExpr <dbl> 10.206344, 10.882607, 9.2…
$ t <dbl> 55.37582, -52.03170, 51.9…
$ P.Value <dbl> 7.906573e-11, 1.243140e-1…
$ adj.P.Val <dbl> 1.726809e-06, 1.726809e-0…
$ B <dbl> 14.41471, 14.16047, 14.15…
火山图如下
火山图
参考资料
Agilent Feature Extraction 12.0 Reference Guide
limma: Linear Models for Microarray and RNA-Seq Data User’s Guide
网友评论