使用DESeq2进行两组间的差异分析
欢迎关注”生信修炼手册”!
DESeq2 接受raw count的定量表格,然后根据样本分组进行差异分析,具体步骤如下
1. 读取数据
读取基因的表达量表格和样本的分组信息两个文件,其中表达量的文件示例如下
gene_id ctrl-1 ctrl-2 ctrl-3 case-1 case-2 case-3
geneA 14 0 11 4 0 12
geneB 125 401 442 175 59 200
每一行为一个基因,每一列代表一个样本。
分组信息的文件示例如下
sample condition
ctrl-1 control
ctrl-2 control
ctrl-3 control
case-1 case
case-2 case
case-3 case
第一列为样本名,第二列为样本的分组信息。
读取文件的代码如下
# 读取表达量的表格
count <- read.table("gene.counts.tsv",header=T,sep="\t",row.names=1,comment.char="",check.names=F)
# 预处理,过滤低丰度的数据
countData <- count[apply(count, 1, sum) > 0 , ]
# 读取样本分组信息
colData <- read.table("sample.group.tsv",header=T,sep="\t",row.names=1,comment.char="",check.names=F)
# 构建DESeq2中的对象
dds <- DESeqDataSetFromMatrix(countData = countData,colData = colData,design = ~ condition)
# 指定哪一组作为control
dds$condition <- relevel(dds$condition, ref = "control")
在读取数据的过程中,有两点需要注意,第一个就是根据表达量对基因进行过滤,通常是过滤低表达量的基因,这一步是可选的,阈值可以自己定义;另外一个就是指定哪一组作为control组,在计算log2FD时 ,需要明确control组,默认会字符串顺序对分组的名字进行排序,排在前面的作为control组,这种默认行为选出的control可能与我们的实验设计不同,所以必须明确指定control组。
2. 归一化
计算每个样本的归一化系数,代码如下
dds <- estimateSizeFactors(dds)
将原始的表达量除以每个样本的归一化系数,就得到了归一化之后的表达量。
3. 估计基因的离散程度
DESeq2假定基因的表达量符合负二项分布,有两个关键参数,总体均值和离散程度α值, 如下图所示
这个α
值衡量的是均值和方差之间的关系,表达式如下
通过下列代码估算每个基因的α值
dds <- estimateDispersions(dds)
4. 差异分析
代码如下
dds <- nbinomWaldTest(dds)
res <- results(dds)
为了简化调用,将第二部到第四部封装到了DESeq
这个函数中,代码如下
dds <- DESeq(dds)
res <- results(dds)
write.table(res,"DESeq2.diff.tsv",sep="\t",quote=F,col.names = NA)
在DESeq2差异分析的过程中,已经考虑到了样本之间已有的差异,所以可以发现,最终结果里的log2FD值和我们拿归一化之后的表达量计算出来的不同, 示意如下
> head(results(dds)[, 1:2])
log2 fold change (MLE): condition B vs A
DataFrame with 6 rows and 2 columnsbaseMean log2FoldChange<numeric> <numeric>
gene1 7.471250 -0.8961954
gene2 18.145279 0.4222174
gene3 2.329461 -2.3216915
gene4 165.634256 -0.1974001
gene5 38.300621 1.3573162
gene6 7.904819 1.8384322
提取归一化之后的表达量,自己计算baseMean和logFoldChange, 示例数据包含了12个样本,其中前6个样本为control, 后6个样本为case , 代码如下
> nor_count <- counts(dds, nor = T)
> res <- data.frame(baseMean = apply(nor_count, 1, mean),log2FoldChange = apply(nor_count, 1, function(t){ mean(t[7:12]) / mean(t[1:6])}))
> head(res)baseMean log2FoldChange
gene1 7.471250 0.5380191
gene2 18.145279 1.3404422
gene3 2.329461 0.1991979
gene4 165.634256 0.8719078
gene5 38.300621 2.5621035
gene6 7.904819 3.5365201
对比DESeq2和自己计算的结果,可以发现,baseMeans是一致的,而log2Foldchange 差异很大,甚至连数值的正负都发生了变化。
log2FD 反映的是不同分组间表达量的差异,这个差异由两部分构成,一种是样本间本身的差异,比如生物学重复样本间基因的表达量就有一定程度的差异,另外一部分就是我们真正感兴趣的,由于分组不同或者实验条件不同造成的差异。
用归一化之后的数值直接计算出的log2FD包含了以上两种差异,而我们真正感兴趣的只有分组不同造成的差异,DESeq2在差异分析的过程中已经考虑到了样本本身的差异,其最终提供的log2FD只包含了分组间的差异,所以会与手动计算的不同。
·end·
—如果喜欢,快分享给你的朋友们吧—
扫描关注微信号,更多精彩内容等着你!
使用DESeq2进行两组间的差异分析相关推荐
- edger多组差异性分析_使用edgeR进行两组间的差异分析
edgeR 接受raw count的定量表格,然后根据样本分组进行差异分析,具体步骤如下 1. 读取文件 需要读取基因在所有样本中的表达量文件,示例如下gene_id ctrl-1 ctrl-2 c ...
- graphpad两组t检验_SPSS如何比较样本两组样本的组内和组间差异(含GraphPad Prism绘图)...
在生物医药统计分析中,经常会遇到这样一类问题.样本分了实验组和对照组,而又同时进行了一种干预手段,在干预前和干预后进行了分别测量数据.这时候既要考虑实验组和对照组之间的比较,又要考虑干预前和干预后的对 ...
- 单细胞测序两组差异分析—seurat包
20210829修改 之前是根据官网+别人帖子写的总结,自己做了一段时间,把之前的再完善一下 尝试使用seurat包进行两组间差异分析 使用的是seurat包自带的数据 创建seurat对象 #首先载 ...
- graphpad如何检测方差齐_SPSS和GraphPad如何比较组内和组间差异「杏花开生物医药统计」...
在生物医药统计分析中,经常会遇到这样一类问题.样本分了实验组和对照组,而又同时进行了一种干预手段,在干预前和干预后进行了分别测量数据.这时候既要考虑实验组和对照组之间的比较,又要考虑干预前和干预后的对 ...
- 【统计学习】一篇文章理解什么是组间差异检验
理解什么是组间差异检验 参数检验与非参数检验 抽样分布 展示差异的常用图表 箱线图(boxplot) 散点图(Scatter plot) 热图(heatmap) 树状图 如何寻找差异? 基于类别标签的 ...
- R堆叠柱状图各成分连线画法:突出展示组间物种丰度变化
作者:朱微金 李陈浩 堆叠柱状图连线画法 提出问题 18年1月29日宏基因组转载了中科院生态中心邓晔组的文章<土壤细菌定量方法结合相对丰度分析揭示种群的真实变化 >.其中的图3基于堆叠柱状 ...
- 扩增子统计绘图6韦恩图:比较组间共有和特有OTU或分类单元
本网对Markdown排版支持较差,对格式不满意的用户请跳转至 或"宏基因组"公众号阅读: 写在前面 优秀的作品都有三部分曲,如骇客帝国.教父.指环王等. 扩增子系列课程也分为三部 ...
- Stata:多个变量组间均值\中位数差异检验
2019暑期Stata现场班,7.17-26日,北京,连玉君+刘瑞明 主讲 作者:韩少真(西北大学) || 刘婉青(西北大学) Stata 连享会: 知乎 | 简书 | 码云 | CSDN ...
- 比较两组数据的差异用什么图更直观_扩增子图表解读7三元图:三组差异数量和关系...
点击上方蓝色「宏基因组」关注我们!专业干货每日推送! 背景介绍(Introduction) 宏基因组学 宏基因组学目前的主要研究方法包括:16S/ITS/18S扩增子.宏基因组.宏转录组和代谢组,其中 ...
最新文章
- 注解 java 反射_java注解和反射
- asp和php数据库怎么区分,asp与php的数据库有哪些区别
- 【实战篇】| 小鹿教你用动态规划撩妹的正确方式
- 2. TypeScript笔记
- 常用HTML标签元素介绍,常用的HTML标签元素总结简介
- R语言观察日志(part3)--repeat循环
- 李春雷 | 夜宿棚花村
- 阿姆斯特朗数_阿姆斯特朗的功能依赖公理 数据库管理系统
- solr 英文模拟mysql like查询xml_Solr实现类似MySQL的LIKE查询功能
- 安卓系统的文件管理神器Solid Explorer(v2.2)
- sprintf 函数
- 《java入门第一季》之面向对象(static关键字内存图解)
- 重t2加权是什么意思_魔兽世界怀旧服:详解盗贼T2.5套装,别犹豫真香
- 在线音乐网站网站开发项目 ,第一篇
- 计算机1级b知识点,初中信息技术等级考试知识点
- fiddler手机模拟器抓包_fiddler抓取手机模拟器数据
- 大众点评评论抓取-加密评论信息完整抓取
- html表单中怎么写年份,HTML表单
- IMAP协议定时监听接收邮件(QQ邮箱、网易邮箱都可)
- nginx中配置二级域名和ssl
热门文章
- 基于微信小程序的图书馆座位预约系统的设计与实现
- 中专计算机寒假作业,中职计算机寒假作业[优质课资]
- Unity 编辑器批量修改Prefab
- 关于HTML colgroup 标签介绍
- vim使用gf(go file)跳转文件
- android+播放器+螺旋效果,螺旋丸特效相机(抖音螺旋丸特效)V2.3.0.3 安卓版
- 北京地铁客流数据特征分析
- Converted Tween Animation Class in AS3
- Gartner 2018新技术成熟度曲线(五大技术趋势)
- 源发行版本 17 需要目标发行版 17