基因芯片数据分析（四）：获取差异表达基因
然而对于大多数生物学工作者而言，学习和使用一种或者多种统计分析手段并不一定非常容易，这需要付出时间和努力。Bioconductor的很多软件包很好的避免了人们为学习统计分析手段而付出的时间。其中最为优秀的软件包就是LIMMA软件包了。
使用limma来分析差异表达的基因，主要分几步走：
- 读取数据
- 预处理数据
- 构建实验设计矩阵
- 使用线性模型估计差异表达的倍数
- 使用贝叶斯平滑标准差
- 试用不同的参数来输出差异表达基因结果。
因为前面几篇文章已经介绍了读取数据以及预处理的相关知识，这里我们直接使用Dilution数据来进行示例。
往期文章
数据预处理
library(affydata)
data(Dilution)
library(limma)
library(gcrma)
##使用gcrma算法来预处理数据
eset <- gcrma(Dilution)实验设计
构建实验设计矩阵，这一步非常关键。通常来讲，这分为两步：
- 给出每个实验样品的相关实验条件
- 设置比对的样品。
需要注意，在给定实验条件时，每一行都要对应到eset中的每一列中的样品。
pData(eset) #我们来看一看样品的排列，这一步的目的在于拿到样品的顺序。
f <- factor(c("liver20", "liver20", "liver10", "liver10"))
design <- model.matrix(~ 0 + f)
design
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 |
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 |
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
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转换，我们后续文章会介绍！
本文分享自 MedBioInfoCloud 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体分享计划 ，欢迎热爱写作的你一起参与！