专栏首页生物信息云基因芯片数据分析(七):edgeR差异分析实战案例

基因芯片数据分析(七):edgeR差异分析实战案例

我们在前2篇文章分别介绍了edgeR和DESeq2包的基本原理:

基因芯片数据分析(五):edgeR包的基本原理

基因芯片数据分析(六):DESeq2包的基本原理

我们接下来通过一个案例介绍利用edgeR和DESeq2包进行差异分析,本文先介绍edgeR。

包的安装和加载

BiocManager::install( "edgeR" )
library( "edgeR" )

读入数据

这里我们用的数据是一个Excel文件,是原始的counts文件(想运行案例,文末获取文件)

# 读取数据
counts <- read.table("gene_counts.xls", sep = "\t", header = T, row.names = 1)

导入的数据我们通过head()函数查看前6行。

行名A1,A2,A3,B1,B2,B3为样本名,列名是基因名。

创建分组

设置实验组别,在基因芯片数据分析(五):edgeR包的基本原理这篇文章中我们介绍基本原理时,有一步需要选择参考样本,在实际分析中,我们可以自己选择参考样本,一般都是对照组作为参考样本,在edgeR中,1代表control样本,2代表case样本。原始数据中有6个样本,对照和实验组各三个。所以我们可以创建一个分组向量。

# 创建分组
group <- c(rep(1,3),rep(2,3))

数据处理

# 创建DGEList类型变量
y <- DGEList(counts=counts, group=group)

这一步相当于创建一列表。注意group中的顺序和counts中行名要对应,也就是对照组和实验组要指定正确。这里A1,A2,A3为control,B1,B2,B3为case。

# 数据过滤
keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes=FALSE]

filterByExpr函数用于确定哪些基因的计数足够大,可以保留在统计分析中。也就是我们在介绍原理(基因芯片数据分析(五):edgeR包的基本原理)中提到的去除表达值为0的基因,实际分析中不是0,我个人理解,counts为只为个位数的也认为是不表达的,所以这里默认的min.count = 10, min.total.count = 15,当然也可以通过参数自己设定。

filterByExpr函数返回的是一个逻辑值类型数据,所以我们需要通过索引获取过滤后的数据。即去除FALSE的基因。

# 计算标准化因子
y <- calcNormFactors(y)
# 查看标准化因子
y$samples
# 计算离散度
y <- estimateDisp(y)
# 显著性检验
et <- exactTest(y)
# 获取排名靠前的基因,这里设置n=100000是为了输出所以基因
et <- topTags(et, n=100000)
# 转换为数据框类型
et <- as.data.frame(et)
# 将行名粘贴为数据框的第一列
et <- cbind(rownames(et),et)
# 指定列名
colnames(et) <- c("gene_id", "log2FoldChange", "log2CPM", "PValue", "FDR")
# 保存数据到本地
write.table(et, "case-vs-control-all.gene.xls", sep = "\t", col.names = TRUE,
            row.names = FALSE, quote = FALSE, na = "")

差异基因筛选

# 差异基因筛选
etSig <- et[which(et$PValue < 0.05 & abs(et$log2FoldChange) > 1),]
# 加入一列,up_down 体现上下调信息
etSig[which(etSig$log2FoldChange > 0), "up_down"] <- "Up"
etSig[which(etSig$log2FoldChange < 0), "up_down"] <- "Down"
# 保存文件
write.table(etSig, "case-vs-control-diff-pval-0.05-FC-2-edgeR.gene.xls", sep = "\t",
            col.names = TRUE, row.names = FALSE, quote = FALSE, na = "")

这样我们就可以获取差异表达的基因了,可以用于后面的作图,比如火山图!

案例文件和代码,后台回复:chip-data-edgeR,获取!

本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

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

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

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

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 蛋白质互作数据库汇总

    Licata L, Briganti L, Peluso D, et al. MINT, the molecular interaction database:...

    DoubleHelix
  • R语言基础教程——第五章:循环结构

    在编程的时候,当我们需要执行一段代码多次时就不可能重复输入该代码多次,这就有了循环编程结构。一般来说,语句按顺序执行。首先执行函数中的第一个语句,然后执行第二个...

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

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

    DoubleHelix
  • 【资源】深度学习 Top100:近 5 年被引用次数最高论文(下载)

    【新智元导读】这里是近5年100篇被引用次数最多的深度学习论文,覆盖了优化/训练方法、无监督/生成模型、卷积网络模型和图像分割/目标检测等十大子领域。重要的论文...

    新智元
  • IJCAI 2018 | 阿里提出IncepText:全新多向场景文本检测模块

    机器之心
  • ASP.NET AJAX(12)__浏览器兼容功能判断浏览器的类型和版本Sys.Browser针对DOM元素的兼容操作针对DOM事件的兼容操作

    目前,常见的浏览器IE(6,8,9),chrome,firefox,safari等,还有国内的一些曾经靠恐吓用户来提高使用率的某浏览器(河蟹社会),这些浏览器对...

    小白哥哥
  • 如何定制 Spring Boot 的 Banner?

    相信用过 Spring Boot 的朋友们一定在启动日志中见过类似如下的内容,比如在启动 Spring Boot 时,控制台默认会打印 Spring Boot ...

    武培轩
  • 一行代码引发的 CI 悲剧

    周五时候,升级通信框架的剥离后,CI主机运行缓慢。增量编译情况下,整个整个流程运行26分钟,以前正常的情况为7-10分钟左右。整个机子卡顿严重。特别是编译环境,...

    DevOps时代
  • 基于 Kubernetes 的微服务项目设计与实现

    随着互联网的发展,后端服务和容器编排技术的日益成熟,微服务成为了后端服务的首选,Kubernetes 也已经成为目前容器编排的事实标准, 微服务拥抱容器时代已经...

    DevOps持续交付
  • Spring boot 自定义banner

    Spring Boot启动的时候会在命令行生成一个banner,其实这个banner是可以自己修改的,本文将会将会讲解如何修改这个banner。

    程序那些事

扫码关注云+社区

领取腾讯云代金券