查看: 1020|回复: 6

[软件使用] 如何进行差异表达分析?

[复制链接]

迅猛龙

Rank: 8Rank: 8

主题
257
注册时间
2020.6.16
在线时间
165 小时

发表于 2021.5.14 09:25:23 | 显示全部楼层 |阅读模式
本帖最后由 基迪奥-Jt桃 于 2021.5.14 09:24 编辑

做差异表达分析,需要一定的R语言基础,推荐大家参加基迪奥的第二期转录组培训班。2021年转录组生信培训班第二期(5月24-28日)火热报名中,培训内容由基础进阶,实操课程提供代码和实时答疑。零生信分析、编程基础也能上手分析,限时招生,赶快报名吧~

培训时间:2021年5月24日-5月28日
授课形式:腾讯会议线上授课

报名费:1800元/人
团购报名:1500元(≥2人)

参加培训班即赠送OmicShare季会员,生信绘图工具免费使用,更可免费观看海量视频课程!

报名方式

方式1:
扫描下方二维码,填写信息报名


方式2:
发送姓名、学校、电话到邮箱contact@genedenovo.com,主题注明“转录组培训班”

客服:020-39341079或18054271626 邓小姐

课程表


言归正传,在上一期的文章中(文章链接:基于edgeR的差异分析案例实操)与大家分享了使用edgeR检测肿瘤组织和正常组织之间的差异表达基因的差异分析案例,本期继续介绍第二个案例↓

案例二

本案例研究拟南芥接种丁香假单胞菌ΔhrcC突变体后产生的反应。在该实验中,用丁香假单胞菌ΔhrcC突变体接种6周龄的拟南芥,然后从叶片中提取总RNA。处理组和对照组各3个重复,共6个样本。

为了方便大家练习使用,本文用到的范例数据已上传到OmicShare,可直接下载使用。

范例数据下载:
https://www.omicshare.com/forum/thread-6985-1-1.html

1.读入数据

  1. #载入edgeR包;
  2. library(edgeR)
  3. arab <- read.csv("case02_arab.csv",row.names = 1)
  4. #查看前6行;
  5. head(arab)
复制代码


  1. #提取样本名(1~4个字符);
  2. Treat <- factor(substring(colnames(arab),1,4))
  3. Treat
  4. #指定"mock"为对照组(处于水平的前面);
  5. Treat <- relevel(Treat, ref="mock")
  6. Treat
  7. #提取样本编号(第5个字符)
  8. Time <- factor(substring(colnames(arab),5,5))
  9. Time

  10. #创建edgeR储存数据的DGElist对象;
  11. y <- DGEList(counts=arab, group=Treat)
  12. y
复制代码


2.数据过滤与标准化

首先,根据表达量对基因进行过滤。counts数非常低的基因几乎对差异分析没有意义,从生物学的角度来看,基因在翻译为蛋白质或具有生物学意义之前,必须具有一定的表达水平。

此外,这些低表达counts的离散性较大,干扰后续的似然估计分析。

因此,在进行进一步分析之前,应过滤掉这些基因。过滤条件一般是每个样本中counts大于5-10个才认为有表达。除了直接粗暴对原始counts过滤,还可以使用filterByExpr()函数对标准化后数据进行自动化过滤。

  1. #推荐根据CPM值进行过滤,此法考虑了文库大小;
  2. keep <- filterByExpr(y)
  3. table(keep)
  4. #FALSE TRUE
  5. #12292 13930
复制代码

Tips
关于过滤标准较复杂,具体可查看filterByExpr()函数的帮助信息:

1. Keeps genes that have count-per-million (CPM) above k in n samples, where k is determined by min.count and by the sample library sizes and n is determined by the design matrix;

2. Each kept gene is required to have at least min.total.count reads across all the samples.


  1. #过滤掉近一半(12292个)的低表达基因;
  2. y <- y[keep, , keep.lib.sizes=FALSE]
  3. dim(y)
  4. #[1] 13930 6

  5. #对数据进行TMM normalization,计算标准化因子;
  6. y <- calcNormFactors(y)
  7. y$samples
复制代码


3.数据探索

  1. #首先,检查样本中的异常值,使用plotMDS()作图;
  2. plotMDS(y, col=rep(1:2, each=3))
复制代码


上图中,样本之间的距离对应于这些样品之间主要生物变异系数(biological coefficient of variation,BCV);从纵轴方向(dimension 2)来看,处理和对照样本可以清楚分开,而在横轴方向上,样本的编号(Time)之间差异较大,即存在明显的批次效应(组内差异大)。

4.生成试验设计矩阵

在进行negative binomial GLMs拟合之前,需要先指定实验设计方案。

  1. #指定blocking factor和比较组,生成design矩阵;这种加性模型适用于配对试验设计或具有批次效应的试验;
  2. design <- model.matrix(~Time+Treat)
  3. rownames(design) <- colnames(y)
  4. design
复制代码


5.估计离散系数(dispersion)

  1. #注意安装"statmod"包否则会报错;
  2. #install.packages("statmod")
  3. y <- estimateDisp(y, design, robust=TRUE)
  4. y$common.dispersion
  5. #[1] 0.06383161

  6. 对common dispersion求平方根即为生物变异的变异系数;这里common dispersion为0.06383161,因此,生物变异系数(coefficient of biological variation)为sqrt(0.06383),约0.25左右。

  7. #可以用BCV plot查看离散系数,图中的红线即为common dispersion;
  8. plotBCV(y)
复制代码


6.差异分析

  1. #先进行quasi-likelihood (QL) glm拟合;
  2. fit <- glmQLFit(y, design, robust=TRUE)
  3. #绘制QL dispersion散点图;
  4. plotQLDisp(fit)
复制代码


  1. #首先,我们检查是否有必要去掉批次效应,我们检测这三次重复之间的差异表达;
  2. qlf <- glmQLFTest(fit, coef=2:3)
  3. #有相当大的差异,证明了校正批次效应的合理性;
  4. topTags(qlf)
复制代码


  1. #计算FDR值;
  2. FDR <- p.adjust(qlf$table$PValue, method="BH")
  3. #统计差异显著个数(TRUE=1);
  4. sum(FDR < 0.05)

  5. #然后重新执行QL F-test,对处理与对照组进行差异分析;
  6. #默认情况下,实验设计矩阵的最后一个coefficient(因子)为处理组;
  7. qlf <- glmQLFTest(fit)
  8. #展示Top10显著表达的差异基因(较小的p值,较大的logFC);
  9. topTags(qlf)

  10. #查看TOP10基因的CPM(counts-per-million)值;
  11. top <- rownames(topTags(qlf))
  12. cpm(y)[top,]
  13. #以FDR=0.05为阈值,统计差异基因个数;
  14. summary(decideTests(qlf))
复制代码


  1. #绘制差异基因散点图(类似火山图);
  2. plotMD(qlf)
  3. abline(h=c(-1,1), col="blue")
复制代码



本文作者:基迪奥-莫北

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
42 小时

发表于 2021.5.14 10:31:40 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
40 小时

发表于 2021.5.15 14:51:01 | 显示全部楼层
加油,加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2018.11.22
在线时间
13 小时

发表于 2021.5.16 11:44:43 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
1
注册时间
2016.8.25
在线时间
65 小时

发表于 2021.5.16 12:20:50 | 显示全部楼层
忘记签到了
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
40 小时

发表于 2021.5.17 08:01:22 | 显示全部楼层
加油,加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2020.2.12
在线时间
51 小时

发表于 2021.5.17 08:40:37 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表