Home
JoVE Journal
Cancer Research
RNA 测序的三种微分表达分析方法:利马、边缘、DESeq2

A subscription to JoVE is required to view this content. Sign in or start your free trial.

Research Article

RNA 测序的三种微分表达分析方法:利马、边缘、DESeq2

DOI: 10.3791/62528

September 18, 2021

, , , , ,

1Department of Obstetrics and Gynecology,Renmin Hospital of Wuhan University, 2Department of Pathology, Shanghai Skin Disease Hospital,Tongji University School of Medicine

Cite Watch Download PDF Download Material list
AI Banner

Please note that some of the translations on this page are AI generated. Click here for the English version.

In This Article

Summary Abstract Introduction Protocol Representative Results Discussion Disclosures Acknowledgements Materials References Reprints and Permissions

Summary

提供了RNA测序的微分表达分析方法的详细协议:伽利马、EdgeR、DESeq2。

Abstract

RNA测序（RNA-seq）是转录学中最广泛使用的技术之一，因为它可以揭示基因改变与复杂的生物过程之间的关系，在肿瘤的诊断、预后和治疗方面具有重要价值。RNA-seq 数据的微分分析对于识别异常转录至关重要，而 limma、EdgeR 和 DESeq2 是微分分析的有效工具。然而，RNA-seq差异分析需要一定的R语言技能和选择适当方法的能力，这是医学教育课程所缺乏的。

在此，我们提供详细的协议，以确定胆管癌 （CHOL） 和正常组织之间通过伽马， DESeq2 和 EdgeR 的差异表达基因 （DEG）， 结果在火山地块和维恩图中显示。伽马、DESeq2 和 EdgeR 这三种方案相似，但在分析过程中有不同的步骤。例如，线性模型用于伽马体的统计，而负二元分布用于边缘R和DESeq2。此外，正常化的RNA-seq计数数据对于 EdgeR 和 limma 是必要的，但对于 DESeq2 来说不是必要的。

在这里，我们为三种微分分析方法提供了详细的协议:伽利马、EdgeR 和 DESeq2。这三种方法的结果是部分重叠的。这三种方法都有各自的优势，方法的选择只取决于数据。

Introduction

RNA测序（RNA-seq）是转录学中应用最广泛的技术之一，具有许多优点（例如，高数据可重复性），并极大地增进了我们对复杂生物过程1、2的功能和动力学的理解。在不同的生物背景下识别异常记录（也称为微分表达基因 ）是RNA-seq分析的关键步骤。RNA-seq 使深入了解发病机制相关的分子机制和生物功能成为可能。因此，差异分析被认为是有价值的诊断，预后和治疗肿瘤3，4，5。目前，更多的开源R/生物导体包已经开发为RNA-seq差分表达分析，特别是利马，DESeq2和EdgeR 1，6，7。然而，差异分析需要一定的R语言技能和选择适当方法的能力，这是医学教育课程所缺乏的。

在本协议中，根据从癌症基因组图集 （TCGA） 中提取的胆管癌 （CHOL） RNA-seq 计数数据，R 程序11分别执行了三种最已知的方法（limma8、EdgeR9和 DESeq210），以确定 CHOL 和正常组织之间的 DEG。伽马、EdgeR 和 DESeq2 的三种方案相似，但在分析过程中有不同的步骤。例如，EdgeR 和 limma8、9需要规范化的 RNA-seq 计数数据，而 DESeq2 则使用自己的库差异来更正数据，而不是校正10。此外，edgeR 特别适用于 RNA-seq 数据，而 limma 则用于微阵列和 RNA-seq。Limma 采用线性模型来评估 DEG12，而 edgeR 中的统计数据基于负二元分布，包括经验贝叶估计、精确测试、通用线性模型和准可能性测试9。

总之，我们分别使用 limma、DESeq2 和 EdgeR 提供 RNA-seq 差分表达分析的详细协议。通过引用本文，用户可以轻松地执行 RNA-seq 差分分析，并为其数据选择适当的差分分析方法。

Protocol

注:打开 R 工作室程序并加载 R 文件"DEGs.R"，该文件可以从补充文件/脚本中获取。

1. 数据的下载和预处理

  1. 从癌症基因组图集 （TCGA） 下载胆管癌 （CHOL） 的高通量测序 （HTSeq） 计数数据。此步骤可以通过以下 R 代码轻松实现。
    1. 单击 "运行" 以安装 R 包。
    2. 单击 "运行" 以加载 R 包。
      如果 （！ 需要命名空间 （"生物经理"， 悄悄地] 真实））
      • 安装.包（"生物经理"）
      生物管理器::安装（c（"TCGAbiolinks"、"总结经验"））
    3. 设置工作目录。
      图书馆（TCGAbiolinks）
      图书馆（总结经验）
      设置（"C:/用户/LIUSHIYI/桌面"）
    4. 选择癌症类型。
      癌症< - "Tcga - chol"
    5. 运行来自"GDCquery.R"文件的 R 代码以下载数据。文件"GDCquery.R"可以从补充文件/脚本中获取:
      来源（"补充文件/脚本/GDCquery.R"）
      头 （cnt）
      ##TCGA-3X-AAVA-01A-11r-A41I-07
      ##ENSG00000000003 4262
      ##ENSG00000000005 1
      ##ENSG00000000419 1254
      ##ENSG00000000457 699
      ##ENSG00000000460 239
      ##ENSG00000000938 334
      注:执行后，将下载 CHOLHTSeq 计数数据并命名为"cnt"，其中行表示合奏基因 ID，列代表样本 ID。请注意样本 ID 中 14-15 位置的数字;数字从01到09表示肿瘤，从10到19表示正常组织。
  2. 将合奏基因 ID 转换为基因符号。
    1. 根据存储路径将注释文件导入 R 中。注释文件（基因代码.v22.注释.gtf）可以从补充文件中获取。
      gtf_v22 <-跟踪层::进口（'补充文件/基因代码.v22.注释.gtf'）
    2. 从"gtf_v22运行 R 代码。R" 文件，可从补充文件/脚本中获取:
      来源（"补充文件/脚本/gtf_v22。R")
    3. 应用功能"ann"将合奏基因 ID 转换为基因符号。
      cnt=ann （gtf_v22）
  3. 过滤低表达基因
    1. 单击 "运行 "以安装 R 包"边缘"。
      生物管理器::安装（"边缘"）
    2. 单击 "运行 "以加载 R 包"边缘"。
      库（边缘）
    3. 运行以下 R 代码，使基因的计数值达到百万分之一（CPM），至少两个样本中的一个以上。
      保持<排（cpm（cnt）>1）>+2
      cnt < - 矩阵 （cnt [保持]）
      注:使用每百万（CPM）值的计数，而不是读数，以消除不同测序深度造成的偏差。

2. 通过"利马"进行差异表达分析

  1. 单击 "运行 "以安装 R 包"limma"。
    生物管理器::安装（"利玛"）
  2. 单击 "运行 "以加载 R 包"limma"、"边缘"。
    图书馆（利玛）
    库（边缘）
  3. 运行以下 R 代码以创建设计矩阵。
    群 <- substring(colnames(cnt),14,15) # Extract group information
    组 [组百分比在% "01"] <- "Cancer" # set '01' as tumor tissue
    组 [组百分比在% "11"] <- "Normal" # set '11' as normal tissue
    群 <- factor (group, levels = c("Normal","Cancer"))
    1. 创建设计矩阵。
      设计<-模型.矩阵（+组）
      行名（设计）<-结名（cnt）
    2. 创建 DGE 列表对象。
      dge < - Dgelist （计数 = cnt， 组 = 组）
    3. 使数据正常化。
      dge < - 钙计算器 （dge， 方法 = "Tmm"）
    4. 运行以下 R 代码以执行基于利马趋势方法的差分表达分析。
      dge
      ##An类"DGEList"的对象
      #$counts
      ##TCGA-3X-AAVA-01A-11r-A41I-07
      ##TSPAN6 4262
      ##DPM1 1254
      ##SCYL3 699
      ##C1orf112 239
      ##FGR 334
    5. 计算 CPM 值。
      日志< - cpm （dge， 日志 = 真实， 前。
    6. 单击 "运行" 以适合线性模型以预测数据或推断变量之间的关系。
      适合<-lmFit（日志，设计）
    7. 根据贝叶斯计算 T 值、F 值和日志赔率。
      适合<- 易趣 （适合， 趋势+ 真实）
    8. 提取结果表。
      res_limma<- 数据. 框架 （顶部表 （适合， n = Inf））

      头（res_limma）
      ## 日志开发 t P. 价值 adj.P.瓦尔 B
      ##RP11-252E2.2 -4.899493 -2.488589 -20.88052 2.386656e-25 4.931786e-21 47.28823
      ##BX842568.1 -4.347930 -2.595205 -20.14532 1.082759e-24 1.118706e-20 45.83656
      ##CTC-537E7.3 -5.154894 -2.143292 -19.59571 3.452354e-24 2.216114e-20 44.72001
      ##RP11-468N14.3 -6.532259 -2.029714 -19.49409 4.289807e-24 2.216114e-20 44.51056
      ##AP006216.5 -4.507051 -2.670915 -19.25649 7.153356e-24 2.956339e-20 44.01704
      ##RP11-669E14.4 -4.107204 -2.828311 -18.93246 1.448209e-23 4.987633e-20 43.33543
      差异表达分析的#The结果保存在"res_limma"中，其中包括基因 ID、日志 2 倍变化值 （logFC）、实验中基因的平均日志 2 表达水平 （AveExpr）、修改的 t 统计 （t）、重新变现 p 值 （P.Value）、错误发现率 （FDR） 更正 p 值 （adj） 。P.Val）和微分表达基因的日志奇数（B）
      注:使用"边缘"的函数"钙定量（）"使数据正常化，以消除样品制备或库构建和测序的影响。在设计矩阵的构造中，有必要将实验设计（例如组织类型:正常组织或肿瘤组织）与矩阵的 ID 样本进行匹配。limma 趋势适用于测序深度相同的数据，而 limma-voom 则适合:（i） 样本库大小不同时:（二） TMM 未正常化的数据;（三） 数据中有很多"噪音"。正日志FC意味着基因在实验中受到调节，而负数意味着基因被调节。
    9. 识别 DEG。
      res_limma美元西格< - 作为因素 （
      伊菲尔斯（res_limma美元）P. Val < 0.05 = 腹肌 （res_limma $logfc） > 2，
      ifelse （res_limma $logFC > 2 ， '向上'， '向下'）， '不是'） * 季< p 值 0.05 和 |log2FC|>= 2 是识别 DEG 的阈值
      摘要（res_limma美元）
      ##down不起来
      ##1880 ​17341 1443
    10. 将结果表输出到文件。
      写.csv （res_limma， 文件 = "result_limma.csv"）
    11. 单击 "运行 "以安装 R 包"ggplot2"。
      安装.包（"ggplot2"）
    12. 单击 "运行 "以加载 R 包"ggplot2"。
      图书馆（ggplot2）
    13. 运行"火山"中的 R 代码。R"创建火山图。文件"火山。R" 可从补充文件中获取。
      来源（"补充文件/脚本/火山。R")
      火山（res_limma，"日志"，"阿杰。P.瓦尔"，2，0.05）
      注:基因可以根据其 log2FC 和 adj-p 值映射到不同位置，向上调节的 DEG 以红色着色，向下调节的 DEG 以绿色着色。
    14. 单击 "导出 "以保存火山图。
      注:火山地块可以以不同格式生成和下载（例如，pdf、TIFF、PNG、JPEG 格式）。基因可以根据其 log2FC 和 adj p 值映射到不同位置，向上调节的 DEG（日志 2FC > 2，adj p < 0.05）为红色，向下调节的 DEG（log2FC < -2，adj p < 0.05）为绿色，非 DEG 为灰色。

3. 通过"边缘"进行差异表达分析

  1. 单击 "运行 "以加载 R 包"边缘"。
    库（边缘）
  2. 运行以下 R 代码以创建设计矩阵。
    组<子弦（同名（cnt），14，15）
    组 [组百分比在% "01"] < - "癌症"
    组 [组百分比在% "11"] < - "正常"
    组=因子（组，级别 = c（"正常"，"癌症"））
    设计<模型.矩阵（+组）
    行名（设计）= 共名（cnt）
  3. 单击 "运行 "以创建 DGEList 对象。
    dge < - Dgelist （计数= cnt）
  4. 使数据正常化。
    dge < - 钙计算器 （dge， 方法 = "Tmm"）
  5. 单击 "运行" 以估计基因表达值的分散性。
    dge <估计放款 （dge， 设计， 坚固 = T）
  6. 单击 "运行 "以适合模型以计数数据。
    适合< - glmqlfit （dge， 设计）
  7. 进行统计测试。
    适合< - 格姆克洛夫特 （适合）
  8. 提取结果表。结果保存在"res_edgeR"中，其中包括日志折叠更改值、日志 CPM、F、p 值和 FDR 更正 p 值。
    res_edgeR=数据.帧（顶部塔格（适合，n=Inf））
    头（res_edgeR）
    ## 日志日志下午 F 价值 Fdr
    ##GCDH -3.299633 5.802700 458.5991 1.441773e-25 2.979280e-21
    ##MSMO1 -3.761400 7.52111 407.0416 1.730539e-24 1.787993e-20R
    ##CL1 -3.829504 5.319641 376.5043 8.652474e-24 5.516791e-20
    ##ADI1 -3.533664 8.211281 372.6671 1.067904e-23 5.516791e-20
    ##KCNN2 -5.583794 3.504017 358.6525 2.342106e-23 9.679455e-20
    ##GLUD1 -3.287447 8.738080 350.0344 3.848408e-23 1.194406e-19
    #The结果保存在"res_edgeR"中，其中包括日志折叠更改值（日志FC）、日志 CPM、F、p 值和 FDR 更正 p 值
  9. 识别 DEG。
    res_edgeR美元 = 作为因素
    > 2， res_edgeR $Fdr < 0.05 + 腹肌 （res_edgeR $logfc）
    如果（res_edgeR$logFC>2，'向上'，'向下'），'不'））
    摘要（res_edgeR美元）
    ##down不起来
    ##1578 15965 3121
  10. 将结果表输出到文件。
    写.csv （res_edgeR， 文件 = "res_edgeR.csv"）
  11. 创建火山图。
    火山 （res_edgeR，"日志"，"FDR"，2，0.05）
  12. 单击 "导出 "以保存火山图。

4. 通过"DESeq2"进行差异表达分析

  1. 单击 "运行 "以安装 R 包"DESEq2"。
    生物管理器::安装（"DESEq2"）
  2. 单击 "运行 "以加载 R 包"DESEq2"。
    图书馆（德塞Q2）
  3. 运行以下 R 代码以确定分组因子。
    组<子弦（同名（cnt），14，15）
    组 [组百分比在% "01"] < - "癌症"
    组 [组百分比在% "11"] < - "正常"
    组=因子（组，级别 = c（"正常"，"癌症"））
  4. 创建解答数据集对象。
    dds < - 德塞克数据集从矩阵 （cnt， 数据框架 （组）， 设计 + 组）
    dds
    ##class: 德塞克数据集
    ##dim: 20664 45
    ##metadata （1）: 版本
    ##assays （1）: 计数
    ##rownames （20664）: TSPAN6 DPM1...RP11-274B21.13 LINC01144
    ##rowData名称 （0）:
    ##colnames （45）: TCGA-3X-AAVA-01A-11r-A41I-07...
    ##colData名称 （1）: 组
  5. 执行分析。
    dds < - 德塞克 （dds）
  6. 生成结果表。
    res_DESeq2 <-数据.框架（结果（dds））

    头（res_DESeq2）
    ## 基础米恩日志 2 折叠改变 lfcse 统计价值帕杰
    ##TSPAN6 4704.9243 -0.8204515 0.3371667 -2.433370 1.495899e-02 2.760180e-02
    ##DPM1 1205.9087 -0.3692497 0.1202418 -3.070894 2.134191e-03 4.838281e-03
    ##SCYL3 954.9772 0.2652530 0.2476441 1.071106 2.841218e-01 3.629059e-01
    ##C1orf112 277.7756 0.7536911 0.2518929 2.992109 2.770575e-03 6.101584e-03
    ##FGR 345.8789 -0.6423198 0.3712729 -1.730047 8.362180e-02 1.266833e-01
    ##CFH 27982.3546 -3.8761382 0.5473363 -7.081823 1.422708e-12 1.673241e-11
    注:结果保存在"res_DESeq2"中，其中包括规范化读数（基数）、日志折叠更改值（日志 2 折叠更改）、日志折叠更改标准错误 （lfcSE）、Wald 统计（统计）、原始 p 值（pvalue）和更正 p 值 （padj） 的均值
  7. 识别 DEG。
    res_DESeq2美元 = 作为因素 （
    res_DESeq2 $padj < 0.05 + 腹肌 （res_DESeq2 $log2 折叠改变） > 2，
    如果（res_DESeq2$log2foldchange>2，'向上'，'向下'），'不'））
    摘要（res_DESeq2美元）
    ##down不起来
    ##1616 16110 2938
  8. 将结果表输出到文件。
    写.csv （res_DESeq2， 文件 = "res_DESeq2.csv"）
  9. 创建火山图。
    火山 （res_DESeq2，" 日志 2 折叠变化"， 帕杰，2，0.05）
  10. 单击 "导出 "以保存火山图。

5. 维恩图

  1. 单击 "运行 "以安装 R 包"文图"。
    安装.包（"文图"）
  2. 单击 "运行 "以加载 R 包"文图"。
    图书馆（文图）
  3. 制作一个向上调节的 DEG 的 Venn 图。
    网格.新页面（）
    网格.绘制（文.图（列表（利马+行名（res_
    利玛 [res_limma美元] "向上"，
    边缘+行名（res_edgeR[res_edgeR美元="向上"，]），
    德塞克2=行名（res_DESeq2[res_DESeq2美元]]
    "向上"，[））
    NULL， 高度 = 3， 宽度 = 3， 单位 = "在"，
    col="黑色"，lwd=0.3，填充=c（"#FF6666"，"#FFFF00"，
    "#993366"),
    阿尔法+c （0.5， 0.5， 0.5），主 = "向上调节的 DEG"）
  4. 单击 "导出 "以保存 Venn 图。
  5. 制作向下调节的 DEG 的 Venn 图。
    网格.新页面（）
    网格.绘制（文.图（列表（利马+行名（res_
    利玛 [res_limma美元] "向下"，
    边缘+ 行名 （res_edgeR [res_edgeR美元]
    "向下"，]）
    DESeq2=行名（res_DESeq2[res_DESeq2美元="向下"，]）
    NULL， 高度 = 3， 宽度 = 3， 单位 = "在"，
    col="黑色"，lwd=0.3，填充=c（"#FF6666"，"#FFFF00"，
    "#993366"),
    alpha=c （0.5， 0.5， 0.5），主 = "向下调节的 DEG"）
  6. 单击 "导出 "以保存 Venn 图。

Representative Results

有各种方法可视化差异表达分析的结果，其中火山图和维恩图特别使用。利马用|logFC|≥2和adj识别了 CHOL 和正常组织之间的 3323 个 DEG。P.Val<0.05作为阈值，其中1880个在CHOL组织中被降低调节，1443个被调高（图1a）。同时，EdgeR 确定了 1578 个下监管 DEG 和 3121 个向上监管的 DEG（图 1b）:DESeq2 确定了 1616 个下行监管的 DEG 和 2938 个向上监管的 DEG（图 1c）。比较这三种方法的结果，1431个上调的DEG和1531个下调节的DEG是重叠的（图2）。

Figure 1
图1。在 CHOL 和正常组织之间识别不同表达的基因 （DEG）。 （a-c）由伽马、边缘和DESeq2获得的所有基因的火山图，分别根据折叠变化（日志2）绘制 ， 红点表示受调节的DEG（调整后 的p 值<0.05和日志|FC|> 2） 和绿点表示向下调节的 DEG（调整后 的 p 值< 0.05 和日志|FC|< 2）。 请单击此处查看此图的较大版本。

Figure 2
图2。Venn 图显示从 limma、边缘和 DESeq2 中得出的结果之间重叠。请单击此处查看此图的较大版本。

补充文件。请点击这里下载此文件。

Discussion

该手稿以前没有出版过，也没有考虑在其他地方出版。所有作者都为撰写这份重要知识内容的手稿作出了贡献，并阅读并批准了最终手稿。我们宣布没有利益冲突。

Disclosures

提供了RNA测序的微分表达分析方法的详细协议:伽利马、EdgeR、DESeq2。

Acknowledgements

这项工作得到了中国国家自然科学基金（81860276号赠款）和国家重点研发计划重点专项资金项目（2018YFC1003200号赠款）的支持。

Materials

R版本 3.6.2免费软件
Rstudio免费软件

References

  1. Tambonis, T., Boareto, M., Leite, V. B. P. Differential Expression Analysis in RNA-seq Data Using a Geometric Approach. Journal of Computational Biology. 25, 1257-1265 (2018).
  2. Wang, Z., Gerstein, M., Snyder, M. RNA-Seq: a revolutionary tool for transcriptomics. Nature Reviews. Genetics. 10, 57-63 (2009).
  3. Anders, S., et al. Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. Nature Protocols. 8, 1765-1786 (2013).
  4. McDermaid, A., Monier, B., Zhao, J., Liu, B., Ma, Q. Interpretation of differential gene expression results of RNA-seq data: review and integration. Briefings in Bioinformatics. 20, 2044-2054 (2019).
  5. Costa-Silva, J., Domingues, D., Lopes, F. M. RNA-Seq differential expression analysis: An extended review and a software tool. PloS One. 12, 0190152 (2017).
  6. Law, C. W., et al. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Research. 5, (2016).
  7. Varet, H., Brillet-Guéguen, L., Coppée, J. Y., Dillies, M. A. SARTools: A DESeq2- and EdgeR-Based R Pipeline for Comprehensive Differential Analysis of RNA-Seq Data. PloS One. 11, 0157022 (2016).
  8. Ritchie, M. E., et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research. 43, 47 (2015).
  9. Robinson, M. D., McCarthy, D. J., Smyth, G. K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 26, 139-140 (2010).
  10. Love, M. I., Huber, W., Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 15, 550 (2014).
  11. Gentleman, R. C., et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biology. 5, 80 (2004).
  12. Law, C. W., Chen, Y., Shi, W., Smyth, G. K. voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology. 15, 29 (2014).
  13. Smyth, G. K. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology. 3, (2004).
  14. Lund, S. P., Nettleton, D., McCarthy, D. J., Smyth, G. K. Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates. Statistical Applications in Genetics and Molecular Biology. 11, (2012).
  15. Reeb, P. D., Steibel, J. P. Evaluating statistical analysis models for RNA sequencing experiments. Frontiers in Genetics. 4, 178 (2013).
  16. Rocke, D. M., et al. Excess False Positive Rates in Methods for Differential Gene Expression Analysis using RNA-Seq Data. bioRxiv. , (2015).
  17. Agarwal, A., et al. Comparison and calibration of transcriptome data from RNA-Seq and tiling arrays. BMC genomics. 11, 383 (2010).
  18. Leng, N., et al. EBSeq: an empirical Bayes hierarchical model for inference in RNA-seq experiments. Bioinformatics. 29, 1035-1043 (2013).

Reprints and Permissions

Request permission to reuse the text or figures of this JoVE article
Request Permission

Play Video

RNA 测序的三种微分表达分析方法:利马、边缘、DESeq2
Contact Us Recommend to Library
Research
Business
Education
Solutions
About JoVE
Others
Contact Us Recommend to Library
JoVE logo

Copyright © 2026 MyJoVE Corporation. All rights reserved

Privacy Terms of Use Policies