基因芯片数据分析（四）：获取差异表达基因

2019-12-13 10:33:48 浏览数 (1)
从基因芯片当中提取生物学的信息需要合理的统计学方法。人们已经为优化传统统计学方法在基因芯片方面的应用做出了多年的努力。但是直到现在，最主要的努力依然还是依据实验设计的差别，用统计学方法提取出差异表达的基因，然后再转回使用实验的方法去验证这个结果。

然而对于大多数生物学工作者而言，学习和使用一种或者多种统计分析手段并不一定非常容易，这需要付出时间和努力。Bioconductor的很多软件包很好的避免了人们为学习统计分析手段而付出的时间。其中最为优秀的软件包就是LIMMA软件包了。

使用limma来分析差异表达的基因，主要分几步走：

  1. 读取数据
  2. 预处理数据
  3. 构建实验设计矩阵
  4. 使用线性模型估计差异表达的倍数
  5. 使用贝叶斯平滑标准差
  6. 试用不同的参数来输出差异表达基因结果。

因为前面几篇文章已经介绍了读取数据以及预处理的相关知识，这里我们直接使用Dilution数据来进行示例。

往期文章

基因芯片数据分析（一）：芯片数据初探

基因芯片数据分析（二）：读取芯片数据

基因芯片数据分析（三）：数据质控

数据预处理

代码语言：javascript复制
library(affydata)
data(Dilution)
library(limma)
library(gcrma)
##使用gcrma算法来预处理数据
eset <- gcrma(Dilution)

实验设计

构建实验设计矩阵，这一步非常关键。通常来讲，这分为两步：

  1. 给出每个实验样品的相关实验条件
  2. 设置比对的样品。

需要注意，在给定实验条件时，每一行都要对应到eset中的每一列中的样品。

代码语言：javascript复制
pData(eset) #我们来看一看样品的排列，这一步的目的在于拿到样品的顺序。
代码语言：javascript复制
f <- factor(c("liver20", "liver20", "liver10", "liver10"))
design <- model.matrix(~ 0   f)
design
代码语言：javascript复制
colnames(design) <- levels(f)
contrast <- c("liver20-liver10")
cont.matrix <- makeContrasts(contrasts = contrast, levels=design)
cont.matrix

对于构建design matrix，我们可以考虑多种实验因素，比如我们有实验：

实验

condition

pair

batch

A

WT

1

1

B

WT

2

2

C

WT

3

1

D

KD

1

2

E

KD

2

1

F

KD

3

2

代码语言：javascript复制
df <- data.frame(Sample=c("A", "B", "C", "D", "E", "F"), 
                 condition=c("WT", 'WT', "WT", "KD", "KD", "KD"),
                 pair = c(1, 2, 3, 1, 2, 3),
                 batch = c(1, 2, 1, 2, 1, 2))
des <- model.matrix(~0   condition   pair   batch, data=df)
des

再如有实验是多个时间点之间的比较：

实验

condition

time

A

WT

0h

B

WT

0h

C

WT

3h

D

WT

3h

E

KD

0h

F

KD

0h

G

KD

3h

H

KD

3h

代码语言：javascript复制
df <- data.frame(Sample=c("A", "B", "C", "D", "E", "F", "G", "H"), 
                 condition=c("WT", 'WT', "WT", "WT", "KD", "KD", "KD", "KD"),
                 time = c(0, 0, 3, 3, 0, 0, 3, 3), stringsAsFactors = FALSE)
df$f <- factor(paste(df$condition, df$time, sep="."))
des <- model.matrix(~0   f, data=df)
colnames(des) <- levels(df$f)
des
代码语言：javascript复制
makeContrasts(contrasts = c("KD.0-WT.0", "KD.3-WT.3", "KD.3-KD.0", "WT.3-WT.0"), levels=des)

对于model.matrix中的formula，下表可以帮助大家理解。

公式

模型

意义

Y ∼ A

Y = β0 β1 A

一次线性关系, 并且允许在 Y 方向上有截矩

Y ∼ −1 A

Y = β1 A

过原点的一次线性关系

Y ∼ A I (A2 )

Y = β0 β1 A β2 A2

多项式模型; 其中 I( ) 函数可以引入普通的数学公式.

Y ∼ A B

Y = β0 β1 A β2 B

二元一次方程式.

Y ∼ A : B

Y = β0 β1 AB

与 A 及 B 的相互关系成一次线性关系

Y ∼ A * B

Y = β0 β1 A β2 B β3 AB

与A和B相关也与A及 B的相互 关系有关也可以写成 Y ∼ A B A:B.

Y ∼ (A B C )^2

Y = β0 β1 A β2 B β3 C β4 AB β5 AC β6 BC

与多个一次变量相关, 同时也与它们的n个元素的组合有关, 这里的 n就是 ( )^n 中的n.也可以写成Y ∼ A*B*C– A:B:C.

但是事实上，生物学工作者最多用到的是前两行，其余的几乎都不会涉及。

差异表达分析

fit <- lmFit(eset, design) fit1 <- contrasts.fit(fit, cont.matrix) fit2 <- eBayes(fit1) topTable(fit2, coef = contrast, n=5)

注释

对于ExpressionSet对象，如果对象中含有featureData的话，当使用topTable等输出结果时，就会自动生成注释信息。所以，基因芯片的注释，多半可以在构建ExpressionSet对象时就准备好。

具体ID转换，我们后续文章会介绍！

数据分析 express

0 人点赞