专栏首页生物信息云基因芯片数据分析(四):获取差异表达基因

基因芯片数据分析(四):获取差异表达基因

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

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

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

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

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

往期文章

基因芯片数据分析(一):芯片数据初探

基因芯片数据分析(二):读取芯片数据

基因芯片数据分析(三):数据质控

数据预处理

library(affydata)
data(Dilution)
library(limma)
library(gcrma)
##使用gcrma算法来预处理数据
eset <- gcrma(Dilution)

实验设计

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

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

需要注意,在给定实验条件时,每一行都要对应到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(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-12-02

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • TCGA数据库:ATAC-Seq数据的下载整理及其可视化

    每一条染色单体由单个线性DNA分子组成。细胞核中的DNA是经过高度有序的包装,否则就是一团乱麻,不利于DNA复制和表达调控。这种有序的状态才能保证基因组的复制和...

    DoubleHelix
  • 不会编程的同学看过来:新的GEO数据在线分析工具——BART

    本地版:https://bitbucket.org/Luisa_amaral/bart

    DoubleHelix
  • PCR中的引物设计

    Nucleotide数据库的搜索结果还是比较智能的,搜索结果前面推荐的就是我们的要的,我们做基因检测,做的PCR是定量PCR,设计引物是以基因的转录本为模板设计...

    DoubleHelix
  • Fiddler抓包11-HTTPS证书Actions无法导出问题

    前言 在点Actions时候出现Export Failed:The root certificate could not be located.最近有很多小伙伴...

    上海-悠悠
  • C语言(各种基本定义)

    数组指针即“指向某个数组的指针”,指针数组即“存放了一堆指针的数组”,函数指针即“指向某个函数的指针”,这些与其说是编程语法,不如说是小学语文。

    用户2617681
  • 聊天机器人開發好消息!!DIALOGFLOW與微訊的天作之合!!

    虽然DIALOGFLOW暂未能够与微讯(WECHAT)或企业微讯(ENTERPRISE WECHAT)进行任何技制上的连接INTERGRATION),确实限制了...

    Spanz
  • pip 升级出错的解决办法

    伪君子
  • Pycharm无法安装第三方库,错误代码Non-zero exit code (1) 的解决方案之pip升级

    问题场景:在pycharm进行安装某些库,install失败,提示需要升级pip ,报错界面问题如下错误代码Non-zero exit code

    测试小兵
  • Using many Decision Trees – random forests使用多棵决策树--随机森林

    In this recipe, we'll use random forests for classification tasks. random forest...

    到不了的都叫做远方
  • [译] MobX 背后的基础原理

    不久之前 Bertalan Miklos 写了一篇很好的博文,比较了 MobX 和基于 proxy 的 NX-framework。这篇博文不仅证明了 proxy...

    江米小枣

扫码关注云+社区

领取腾讯云代金券