以解牛之法析生信观微雀之形览科研。今天我们一起来学习使用DESeq2进行转录组的差异分析。三、DEseq2差异分析及对比使用数据airway工具DESeq2相关R包DESeq2分析、ggplot2绘图1. 工具介绍DESeq2DESeq2是 RNA-seq 差异表达分析 最主流的 R 包之一发表于 2014 年已被引用数万次专门用于基于负二项分布模型的计数数据差异分析。它由 Bioconductor 项目维护是转录组分析的标准工具。输入数据原始整数计数矩阵不能是 TPM/FPKM行基因列样本标准化方法 median-of-ratios 方法通过计算每个基因的相对计数因子来校正文库大小差异对异常值更稳健差异检验 基于 Wald 检验 或 似然比检验LRT自动估计基因离散度处理小样本 使用 经验贝叶斯方法 将基因间离散度信息“收缩”shrinkage使小样本如 n3也能获得可靠结果log2FC 收缩 提供 apeglm 和 ashr 两种收缩方法对低计数基因的 log2FC 进行适度压缩提高结果稳定性DESeq2说明文档Bioconductor 手册https://bioconductor.org/packages/devel/bioc/manuals/DESeq2/man/DESeq2.pdfgithub仓库https://github.com/thelovelab/DESeq2使用说明# 如何使用DESeq(object,testc(Wald,LRT),fitTypec(parametric,local,mean,glmGamPoi),sfTypec(ratio,poscounts,iterate),betaPrior,fulldesign(object),reduced,quietFALSE,minReplicatesForReplace7,modelMatrixType,useTFALSE,minmuif(fitTypeglmGamPoi)1e-06else0.5,parallelFALSE,BPPARAMbpparam())2. 下载及加载包# # 第一次安装可以选中ctrlshiftc取消注释运行安装# 安装biocManager# if (!require(BiocManager, quietly TRUE))# install.packages(BiocManager)# # 差异分析核心DESeq2包# BiocManager::install(c(DESeq2, apeglm), update FALSE)# # 数据处理# BiocManager::install(c(tidyverse, tibble), update FALSE)# # 绘图包# BiocManager::install(c(ggplot2, EnhancedVolcano, ggrepel, pheatmap, ggpubr, RColorBrewer), update FALSE)# 1.2.3 加载包确定已经安装了包后再加载# 如果报错或者提示没有安装返回上一步library(DESeq2)library(apeglm)library(tidyverse)library(ggplot2)library(EnhancedVolcano)library(ggrepel)library(pheatmap)library(ggpubr)3. 示例数据加载airway数据if(!require(BiocManager,quietlyTRUE))install.packages(BiocManager)options(BioC_mirrorhttps://mirrors.ustc.edu.cn/bioc/)# 切换中科大镜像下载BiocManager::install(airway)#下载airway数据2014年发表library(airway)#加载包data(airway)#加载airway数据集airway#查看数据集是否加载成功Bioconductor 页面https://bioconductor.org/packages/release/data/experiment/html/airway.html原始 GEO 项目https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?accGSE52778airway 数据集来源于 GEO 项目 GSE52778是该项目的子集包含 4 个细胞系各 1 对处理/对照样本共 8 个样本数据处理方式与原始研究一致。注如果加载不成功也可以手动下载但是此步需要进一步处理才可使用。地址https://bioconductor.org/packages/release/data/experiment/html/airway.html或者使用下载的GSE52778但是也需要适当的处理才可使用。4. 标准的 DESeq2 差异分析流程可从airway中获取原始整数技术矩阵并构建分组信息——过滤——构建 DESeq2 对象——运行 DESeq2进行差异基因的比对。# 标准的分析流程# 2.1 提取数据library(airway)data(airway)count-assay(airway)# 原始整数计数矩阵col_data-colData(airway)# 样本注释含 cell 和 dex# 构建分组信息也可以直接从 col_data 提取group-data.frame(samplecolnames(count),groupcol_data$dex,# trt / untrtcellcol_data$cell# 细胞系信息用于校正批次)# 确认样本顺序一致必须为 TRUEstopifnot(all(colnames(count)group$sample))# 2.2 过滤低表达基因RNA-seq 标准操作提高统计功效# 保留至少在 3 个样本中 count 10 的基因可根据实际情况调整阈值keep-rowSums(count10)3count_filter-count[keep,]cat(过滤前基因数,nrow(count),\n)cat(过滤后基因数,nrow(count_filter),\n)# 2.3 构建 DESeq2 对象配对设计校正细胞系效应dds-DESeqDataSetFromMatrix(countDatacount_filter,# 使用过滤后的矩阵colDatagroup,design~cellgroup# 加入 cell 作为协变量)# 明确指定对照组确保 log2FC 方向正确dds$group-relevel(dds$group,refuntrt)# 2.4 运行 DESeq2dds-DESeq(dds)# 查看可用的系数名用于后续对比resultsNames(dds)可以看到过滤前后基因数为过滤前基因数 63677 过滤后基因数 16596 。此步完成了差异分析核心包为DESeq2下步进行差异分析的提取。5. 提取差异基因提取差异基因并且分别保存全部的分析结果和显著的差异基因。# 2.5 提取两组比较的结果# contrast 参数c(分组列名, 实验组, 对照组)res-results(dds,contrastc(group,trt,untrt))# 按 p 值排序最显著的在前res-res[order(res$pvalue),]# 转为数据框并添加基因名列res_df-as.data.frame(res)res_df$gene-rownames(res_df)# 保存全部基因的差异结果write.csv(res_df,fileDESeq2_all_gene_result.csv,row.namesFALSE)# 2.6 筛选显著差异基因阈值|log2FC| 1 且 padj 0.05res_df$change-Stableres_df$change[res_df$log2FoldChange1res_df$padj0.05]-Upres_df$change[res_df$log2FoldChange-1res_df$padj0.05]-Down# 统计上下调基因数量table(res_df$change)cat(显著上调基因数,sum(res_df$changeUp),\n)cat(显著下调基因数,sum(res_df$changeDown),\n)# 2.7 可选保存显著差异基因子集sig_genes-res_df[res_df$change!Stable,]write.csv(sig_genes,fileDESeq2_sig_genes.csv,row.namesFALSE)# 打印前几个显著基因head(res_df[res_df$change!Stable,])可以看到 有显著上调基因511下调基因489。运行到这里我们完成了提取差异基因并且保存全部的分析结果和显著的差异基因结果。后续绘制火山图、富集分析就可直接使用这数据啦今天的生信知识分享就到这里啦咱们下期见以解牛之法析生信观微雀之形览科研。