RNA-seq 详细教程:搞定count归一化(5)

news2025/1/12 7:49:47

学习目标

  1. 了解如何在归一化过程中列出不同的 uninteresting factors(无关因素)
  2. 了解常用的归一化方法,已经如何使用
  3. 了解如何创建 DESeqDataSet 对象及其结构
  4. 了解如何使用 DESeq2 进行归一化

1. 归一化

差异表达分析工作流程的第一步是计数归一化,这是对样本之间的基因表达进行准确比较所必需的。

Normalization
Normalization

每个基因的映射读数计数是 RNA 表达以及许多其他因素的结果。归一化是调整原始计数值以解决“无关”因素的过程。以这种方式,表达水平在样本之间或样本内更具可比性。

在归一化过程中经常考虑的“无关”因素:

1.1. 测序深度

考虑测序深度对于比较样本之间的基因表达是必要的。在下面的示例中,每个基因在样本 A 中的表达似乎是样本 B 的两倍。然而,这是样本 A 的测序深度加倍的结果。

sequencing depth
sequencing depth

1.2. 基因长度

计算基因长度对于比较同一样本中不同基因之间的表达是必要的。在下面的示例中,基因 X 和基因 Y 具有相似的表达水平,但映射到基因 X 的读数数量将比映射到基因 Y 的读数多得多,因为基因 X 更长。

Gene length
Gene length

1.3. RNA组成

样本之间一些高度差异表达的基因、样本之间表达的基因数量的差异或污染的存在可能会扭曲某些类型的归一化方法。建议考虑 RNA 组成以准确比较样本之间的表达,这在进行差异表达分析时尤为重要。

在下面的示例中,假设样本 A 和样本 B 之间的测序深度相似,并且除了基因差异表达之外的每个基因在样本之间呈现相似的表达水平。样本 B 中的计数会受到 差异表达基因的极大影响,它占据了大部分计数。因此,样本 B 的其他基因的表达似乎低于样本 A 中的相同基因。

RNA composition
RNA composition

归一化不仅对于差异表达分析必不可少,对于探索数据分析、数据可视化以及探索或比较样本之间或样本内的计数也是必要的。

2. 归一化方法

几种常见的归一化方法:

方法描述考虑因素使用场景
CPM (counts per million)按照reads总数缩放计数测序深度同一样本组重复之间的基因计数比较;不适用于样本内比较或差异表达分析
TPM (transcripts per kilobase million)每百万读取reads比对的转录本长度 (kb) 计数测序深度与基因长度样本内或同一样本组样本之间的基因计数比较;不适用于差异表达分析
RPKM/FPKM (reads/fragments per kilobase of exon per million reads/fragments mapped)类似于TPM测序深度与基因长度样本中基因之间的基因计数比较;不适用于样本比较或差异表达分析
DESeq2’s median of ratios计数除以特定于样本的大小因子,该因子由基因计数相对于每个基因的几何平均值的中位数比率确定测序深度和RNA组成样品之间的基因计数比较和差异表达分析;不适用于样本内比较
EdgeR’s trimmed mean of M values (TMM)使用样本之间对数表达比率的加权修剪平均值测序深度和RNA组成样品之间的基因计数比较和差异表达分析;不适用于样本内比较
  • RPKM/FPKM:不推荐用于样本间比较

虽然 TPM 和 RPKM/FPKM 归一化方法都考虑了测序深度和基因长度,但不推荐使用 RPKM/FPKM。原因是RPKM/FPKM方法输出的归一化计数值在样本之间没有可比性。

使用 RPKM/FPKM 归一化,每个样本的 RPKM/FPKM 归一化计数总数会有所不同。因此,您不能在样本之间平均比较每个基因的归一化计数。

RPKM-归一化计数表:

genesampleAsampleB
XCR15.55.5
WASHC173.421.8
Total RPKM-normalized counts1,000,0001,500,000

例如,在上表中,样本 A 的 XCR1 (5.5/1,000,000) 计数比例高于样本 B (5.5/1,500,000),即使 RPKM 计数值相同。因此,我们不能直接比较样本 A 和样本 B 之间 XCR1(或任何其他基因)的计数,因为样本之间的归一化计数总数不同。

  • DESeq2-归一化计数:比率方法的中值(Median of ratios method)

由于用于差异表达分析的工具正在比较样本组之间相同基因的计数,因此该工具不需要考虑基因长度。然而,确实需要考虑测序深度和 RNA 组成。为了标准化测序深度和 RNA 组成,DESeq2 使用比率中值方法。在用户端只有一个步骤,但在后端涉及多个步骤,如下所述。

  1. 创建一个伪参考样本(逐行几何平均值)

对于每个基因,都会创建一个伪参考样本,该样本等于所有样本的几何平均值。

genesampleAsampleBpseudo-reference sample
EF2A1489906sqrt(1489 * 906) = 1161.5
ABCD12213sqrt(22 * 13) = 17.7
  1. 计算每个样本与参考的比率

对于每个样本中的每个基因,计算比率(样本/参考)(如下所示)。由于大多数基因没有差异表达,因此每个样本中的大多数基因在样本中的比例应该相似。

genesampleAsampleBpseudo-reference sampleratio of sampleA/refratio of sampleB/ref
EF2A14899061161.51489/1161.5 = 1.28906/1161.5 = 0.78
ABCD1221316.922/16.9 = 1.3013/16.9 = 0.77
MEFV793410570.2793/570.2 = 1.39410/570.2 = 0.72
BAG1764256.576/56.5 = 1.3542/56.5 = 0.74
MOV105211196883.7521/883.7 = 0.5901196/883.7 = 1.35
  1. 计算每个样本的归一化因子(大小因子)

给定样本的所有比率的中值(上表中的列)被视为该样本的归一化因子(大小因子),计算如下。请注意,差异表达的基因不应影响中值:

normalization_factor_sampleA <- median(c(1.28, 1.3, 1.39, 1.35, 0.59))

normalization_factor_sampleB <- median(c(0.78, 0.77, 0.72, 0.74, 1.35))

下图说明了单个样本的所有基因比率分布的中值(y 轴是频率)。

figure
figure

比率中位数法假设并非所有基因都差异表达;因此,归一化因子应考虑样本的测序深度和 RNA 组成(大的离群基因不会影响中值比率值)。该方法对上调/下调和大量差异表达基因的不平衡具有鲁棒性。

  1. 使用归一化因子计算归一化计数值

这是通过将给定样本中的每个原始计数值除以该样本的归一化因子来执行的,生成归一化计数值。这是针对所有计数值(每个样本中的每个基因)执行的。例如,如果样本 A 的中值比率为 1.3,样本 B 的中值比率为 0.77,则可以按如下方式计算归一化计数:

  • Raw Counts
genesampleAsampleB
EF2A1489906
ABCD12213
  • Normalized Counts
genesampleAsampleB
EF2A1489 / 1.3 = 1145.39906 / 0.77 = 1176.62
ABCD122 / 1.3 = 16.9213 / 0.77 = 16.88

请注意,归一化计数值不是整数。


以上步骤仅作为演示,在实际使用DESeq2过程中,只需要一步命令,即可完成计算。

3. Mov10 归一化

现在我们知道了计数归一化的理论,我们将使用 DESeq2Mov10 数据集的计数进行归一化。这需要几个步骤:

  1. 确保 metadata 数据框的行名存在,并且与 counts 数据框的列名顺序相同。
  2. 创建一个 DESeqDataSet 对象
  3. 生成归一化 counts

3.1. 数据匹配

我们应该始终确保样本名称在两个文件之间匹配,并且样本的顺序相同。如果不是这种情况,DESeq2 将输出错误。

# 检查两个文件中的样本名称是否匹配
all(colnames(txi$counts) %in% rownames(meta))
all(colnames(txi$counts) == rownames(meta))

如果数据不匹配,可以使用 match() 函数重新排列它们。

3.2. 创建对象

让我们从创建 DESeqDataSet 对象开始,然后可以更多地讨论其中存储的内容。要创建对象,我们需要将计数矩阵和元数据表作为输入。我们还需要指定一个设计公式。设计公式指定元数据表中的列以及它们在分析中的使用方式。对于我们的数据集,我们只有一列感兴趣,即 ~sampletype。此列具有三个因子水平,它告诉 DESeq2 对于每个基因,我们要评估相对于这些不同水平的基因表达变化。

我们的计数矩阵输入存储在 txi 列表对象中。所以我们需要指定使用 DESeqDataSetFromTximport() 函数,这将提取计数组件并将值四舍五入到最接近的整数。

# 对象创建
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ sampletype)

3.3. 生成归一化counts

下一步是对计数数据进行归一化,以便在样本之间进行正确的基因比较。

normalize
normalize

为了执行归一化比率方法的中位数,DESeq2 有一个 estimateSizeFactors() 函数可以生成大小因子。我们将在下面的示例中演示此功能,但在典型的 RNA-seq 分析中,此步骤由 DESeq() 函数自动执行,我们将在后面讨论。

dds <- estimateSizeFactors(dds)

通过将结果分配回 dds 对象,我们用适当的信息填充了 DESeqDataSet 对象。我们可以使用以下方法查看每个样本的归一化因子:

sizeFactors(dds)

现在,要从 dds 中检索归一化计数矩阵,我们使用 counts() 函数并添加参数 normalized=TRUE

normalized_counts <- counts(dds, normalized=TRUE)

我们可以将这个归一化的数据矩阵保存到文件中以备后用:

write.table(normalized_counts, file="data/normalized_counts.txt", sep="\t", quote=F, col.names=NA)

注意:DESeq2 实际上并不使用归一化计数,而是使用原始计数并对广义线性模型 (GLM) 中的归一化进行建模。这些归一化计数对于结果的下游可视化很有用,但不能用作 DESeq2 或任何其他使用负二项式模型执行差异表达分析的工具的输入。


欢迎Star -> 学习目录

更多教程 -> 学习目录


本文由 mdnice 多平台发布

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/60340.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

exploit-db图文教程

一、ExploitDB 简介 ExploitDB 是一个面向全世界黑客的漏洞提交平台&#xff0c;该平台会公布最新漏洞的相关情况&#xff0c;这些可以帮助企业改善公司的安全状况&#xff0c;同时也以帮助安全研究者和渗透测试工程师更好的进行安全测试工作。Exploit-DB 提供一整套庞大的归档…

Spring框架(十):Spring注解开发配置MyBatis框架等第三方框架

Spring注解开发配置MyBatis框架等第三方框架引子注解配置MyBatis注解整合Mapper的原理Import注解整合第三方框架引子 痛定思痛&#xff0c;主要问题出现在自己雀氏不熟悉框架底层、一些面试题&#xff0c;以及sql的一些情况淡忘了。 本章节的开始是对于过去的重新回顾&#xf…

IDEA常用插件(持续更新中......)

再使用idea办公后&#xff0c;会发现安装一些插件&#xff0c;对我们来说会事半功倍。同时也会让同事也开始情不自禁的嘚瑟了。 1.Free Mybatis plugin mybatis和mybatis-plus基本上是目前最主流的ORM框架了&#xff0c;相比于Hibernate更加灵活&#xff0c;性能也更好。所以…

爬虫入门——1、基本概念

1、爬虫是什么 网络爬虫&#xff08;又称网络机器人&#xff09;&#xff0c;是一种按照一定的规则&#xff0c;自动地抓取网络信息的程序或者脚本。 通俗地讲&#xff0c;我们把互联网比作一张大蜘蛛网&#xff0c;每个站点资源比作蜘蛛网上的一个结点&#xff0c;爬虫就像一…

【spark】第二章——SparkCore之运行架构及核心编程

文章目录1. Spark 运行架构1.1 1 运行架构1.2 核心组件1.2.1 Driver1.2.2 Executor1.2.3 Master & Worker1.2.4 ApplicationMaster1.3 核心概念1.3.1 Executor 与 Core1.3.2 并行度&#xff08;Parallelism&#xff09;1.3.3 有向无环图&#xff08;DAG&#xff09;1.4 提交…

Compiere的应用字典介绍

模型驱动架构介绍 本章介绍了Compiere的模型驱动架构和Compiere的数据字典功能。 在大多数应用程序中&#xff0c;开发人员必须设计代码并测试每个屏幕。这可能是非常耗时的&#xff0c;并导致整个应用程序在外观和感觉以及功能方面的不一致。 这也会使用户难以学习像ERP这样复…

H5 Canvas 垂直箭头绘制

效果 ⚠ 因为使用的是斜率来处理的垂直逻辑 tan&#xff0c;当为被除数为0时做了特殊处理&#xff0c;两点自由变换时到达零界点会有卡顿。 推导 开始复习初中二年级数学知识 斜率k的公式&#xff1a;k(y1−y2)(x1−x2)k \dfrac{(y_1 -y_2)}{(x_1 - x_2)}k(x1​−x2​)(y1​…

【三维目标检测】VoteNet(一)

VoteNet是用于点云三维目标检测模型算法&#xff0c;发表在ICCV 2019《Deep Hough Voting for 3D Object Detection in Point Clouds》&#xff0c;论文地址为“https://arxiv.org/abs/1904.09664”。VoteNet核心思想在于通过霍夫投票的方法实现了端到端3D对象检测网络&#xf…

LeetCode 852. 山脉数组的峰顶索引

&#x1f308;&#x1f308;&#x1f604;&#x1f604; 欢迎来到茶色岛独家岛屿&#xff0c;本期将为大家揭晓LeetCode 852. 山脉数组的峰顶索引 &#xff0c;做好准备了么&#xff0c;那么开始吧。 &#x1f332;&#x1f332;&#x1f434;&#x1f434; 一、题目名称 Leet…

[Error]适配iPad时调用UIAlertController和UIActivityViewController软件崩溃问题

问题&#xff1a; 适配iPad时&#xff0c;调用UIAlertController和UIActivityViewController软件崩溃。iPhone设备上软件运行正常。 错误打印如下&#xff1a; Thread 1: "UIPopoverPresentationController (<UIPopoverPresentationController: 0x12f33b020>) sho…

管道模式 流处理

&#xff08;一&#xff09;介绍 管道这个名字源于自来水厂的原水处理过程。原水要经过管道&#xff0c;一层层地过滤、沉淀、去杂质、消毒&#xff0c;到管道另一端形成纯净水。我们不应该把所有原水的过滤都放在一个管道中去提纯&#xff0c;而应该把处理过程进行划分&#…

Pytorch Bert 中文分类 运行代码时候遇到的问题

问题1 bert AutoModel.from_pretrained(bert-base-chinese) 报错信息如下&#xff1a; RuntimeError: Error(s) in loading state_dict for BertModel: size mismatch for bert.embeddings.word_embeddings.weight: copying a param with shape torch.Size([21128, 768…

m基于GA遗传算法的PMSM永磁同步电机参数最优计算matlab仿真

目录 1.算法描述 2.仿真效果预览 3.MATLAB核心程序 4.完整MATLAB 1.算法描述 永磁同步电机&#xff08;PMSM&#xff09;基本结构为定子、转子和端盖。其中转子磁路结构是永磁同步电机与其它电机最主要的区别&#xff0c;其在很大程度上决定了永磁同步电机的实际性能指标。…

AtCoder Beginner Contest 280 老年人复建赛

好久没更新了&#xff0c;因为最近p事实在是有点多&#xff0c;让人心烦意乱 还是安安心心打比赛舒服 A&#xff0c;B&#xff0c;C就不讲啦 D - Factorial and Multiple 大意&#xff1a; 给定一个数字k<1e12&#xff0c;求最小的数字n满足n!%k0; 思路1&#xff1a; 不…

hadoop完全分布式环境搭建详细版

1. hadoop集群规划 1.准备3台客户机(关闭防火墙&#xff0c;静态ip&#xff0c;主机名称) 2.安装jdk 3.配置环境变量 4.安装hadoop&#xff0c;hadoop版本是3.1.3,包名为hadoop-3.1.3.tar.gz 5.配置环境变量 6.配置集群 7.单点启动 8.配置ssh 9.群起集群并测试集群 注意: NameN…

Ubuntu20.04静态路由表连通局域网各网段主机 Vmware WorkStation

文章目录示例拓扑虚拟机的三种网络模式虚拟网络编辑器的设置虚拟主机与虚拟路由设置细节Ubuntu20.04设置静态IP给R1添加双网卡给R1、R2开启转发功能配置路由表References示例拓扑 宿主机是Windows11 PC与Router均为 Ubuntu20.04系统。 虚拟机的三种网络模式 虚拟机默认是只初…

如何利用InVest模型估算区域产水量

1.什么是InVEST模型 InVEST模型&#xff08;Integrated Valuation of Ecosystem Services and Tradeoffs &#xff09;是生态系统服务评估与权衡模型的简称&#xff0c;是美国自然资本项目组开发的、用于评估生态系统服务功能量及其经济价值、支持生态系统管理和决策的一套模型…

十四、使用 Vue Router 开发单页应用(1)

本章概要 感受前端路由 HTML 使用路由模块开发使用路由 传统的 Web 应用程序不同页面间的跳转都是向服务器发起请求&#xff0c;服务器处理请求后向浏览器推送页面。 在单页应用程序中&#xff0c;不同视图&#xff08;组件的模板&#xff09;的内容都是在同一个页面中渲染&…

golang 琐碎知识

golang 琐碎知识&#xff08;持续进行&#xff09; 时间格式 time.now.Format("2006-01-02T 15:04:05")make声明切片bug Golang&#xff1a;statusList : make([]*model.StatusList, 6) 会声明一个长为6的null切片&#xff0c;使用append添加时不会将null覆盖掉去掉切…

JMeter入门教程(10) --函数助手

文章目录1.CSVRead2.Random3.RandomString4.RandomDate5.time在JMeter的选项菜单中有一个“函数助手对话框”&#xff0c;点击打开“函数助手”对话框&#xff0c;使用函数助手&#xff0c;我们可以从“选择一个功能”下拉列表中选择一个函数&#xff0c;并为其参数设定值。表格…