RNA-seq 保姆教程：差异表达分析（二）

news2025/1/10 3:29:52

介绍

RNA-seq 目前是测量细胞反应的最突出的方法之一。RNA-seq 不仅能够分析样本之间基因表达的差异，还可以发现新的亚型并分析 SNP 变异。本教程[1]将涵盖处理和分析差异基因表达数据的基本工作流程，旨在提供设置环境和运行比对工具的通用方法。由于完整版过长，因此分为两部分，需要获取完整版的，请跳转文末。

7. 差异分析

  • 将基因计数导入 R/RStudio

工作流程完成后，您现在可以使用基因计数表作为 DESeq2 的输入，使用 R 语言进行统计分析。

7.1. 安装R包

source("https://bioconductor.org/biocLite.R")
biocLite("DESeq2") ; library(DESeq2)
biocLite("ggplot2") ; library(ggplot2)
biocLite("clusterProfiler") ; library(clusterProfiler)
biocLite("biomaRt") ; library(biomaRt)
biocLite("ReactomePA") ; library(ReactomePA)
biocLite("DOSE") ; library(DOSE)
biocLite("KEGG.db") ; library(KEGG.db)
biocLite("org.Mm.eg.db") ; library(org.Mm.eg.db)
biocLite("org.Hs.eg.db") ; library(org.Hs.eg.db)
biocLite("pheatmap") ; library(pheatmap)
biocLite("genefilter") ; library(genefilter)
biocLite("RColorBrewer") ; library(RColorBrewer)
biocLite("GO.db") ; library(GO.db)
biocLite("topGO") ; library(topGO)
biocLite("dplyr") ; library(dplyr)
biocLite("gage") ; library(gage)
biocLite("ggsci") ; library(ggsci)

7.2. 导入表达矩阵

  • 开始导入文件夹中的 featureCounts 表。本教程将使用 DESeq2 对样本组之间进行归一化和执行统计分析。 
# 导入基因计数表
# 使行名成为基因标识符
countdata <- read.table("example/final_counts.txt", header = TRUE, skip = 1, row.names = 1)

# 从列标识符中删除 .bam 和 '..'
colnames(countdata) <- gsub(".bam", "", colnames(countdata), fixed = T)
colnames(countdata) <- gsub(".bam", "", colnames(countdata), fixed = T)
colnames(countdata) <- gsub("..", "", colnames(countdata), fixed = T)

# 删除长度字符列
countdata <- countdata[ ,c(-1:-5)]

# 查看 ID
head(countdata)  # 如下图
countdata
countdata

7.3. 导入metadata

  • 导入元数据文本文件。 SampleID 必须是第一列。 
# 导入元数据文件
# 使行名称与 countdata 中的 sampleID 相匹配
metadata <- read.delim("example/metadata.txt", row.names = 1)

# 将 sampleID 添加到映射文件
metadata$sampleid <- row.names(metadata)

# 重新排序 sampleID 以匹配 featureCounts 列顺序。
metadata <- metadata[match(colnames(countdata), metadata$sampleid), ]

# 查看 ID
head(metadata)  # 如下图
metadata
metadata

7.4. DESeq2对象

  • 根据计数和元数据创建 DESeq2 对象 
# - countData : 基于表达矩阵
# - colData : 见上图
# - design : 比较
ddsMat <- DESeqDataSetFromMatrix(countData = countdata,
                                 colData = metadata,
                                 design = ~Group)


# 查找差异表达基因
ddsMat <- DESeq(ddsMat)

7.5. 统计

  • 获取基因数量的基本统计数据 
# 使用 FDR 调整 p-values 从检测中获取结果
results <- results(ddsMat, pAdjustMethod = "fdr", alpha = 0.05)

# 结果查看
summary(results)  # 如下图
results
results 
# 检查 log2 fold change
## Log2 fold change is set as (LoGlu / HiGlu)
## Postive fold changes = Increased in LoGlu
## Negative fold changes = Decreased in LoGlu
mcols(results, use.names = T)  # 结果如下
mcols_result
mcols_result

8. 注释基因symbol

经过比对和总结，我们只有带注释的基因符号。要获得有关基因的更多信息，我们可以使用带注释的数据库将基因符号转换为完整的基因名称和 entrez ID 以进行进一步分析。

  • 收集基因注释信息 
# 小鼠基因组数据库
library(org.Mm.eg.db) 

# 添加基因全名
results$description <- mapIds(x = org.Mm.eg.db,
                              keys = row.names(results),
                              column = "GENENAME",
                              keytype = "SYMBOL",
                              multiVals = "first")

# 添加基因 symbol
results$symbol <- row.names(results)

# 添加 ENTREZ ID
results$entrez <- mapIds(x = org.Mm.eg.db,
                         keys = row.names(results),
                         column = "ENTREZID",
                         keytype = "SYMBOL",
                         multiVals = "first")

# 添加 ENSEMBL
results$ensembl <- mapIds(x = org.Mm.eg.db,
                          keys = row.names(results),
                          column = "ENSEMBL",
                          keytype = "SYMBOL",
                          multiVals = "first")

# 取 (q < 0.05) 的基因
results_sig <- subset(results, padj < 0.05)

# 查看结果
head(results_sig)  # 如下图
alt
  • 将所有重要结果写入 .txt 文件 
# 将归一化基因计数写入 .txt 文件
write.table(x = as.data.frame(counts(ddsMat), normalized = T), 
            file = 'normalized_counts.txt', 
            sep = '\t', 
            quote = F,
            col.names = NA)

# 将标准化基因计数写入 .txt 文件
write.table(x = counts(ddsMat[row.names(results_sig)], normalized = T), 
            file = 'normalized_counts_significant.txt', 
            sep = '\t', 
            quote = F, 
            col.names = NA)

# 将带注释的结果表写入 .txt 文件
write.table(x = as.data.frame(results), 
            file = "results_gene_annotated.txt", 
            sep = '\t', 
            quote = F,
            col.names = NA)

# 将重要的注释结果表写入 .txt 文件
write.table(x = as.data.frame(results_sig), 
            file = "results_gene_annotated_significant.txt", 
            sep = '\t', 
            quote = F,
            col.names = NA)

9. 绘图

有多种方法可以绘制基因表达数据。下面只列出了一些流行的方法。

9.1. PCA

# 将所有样本转换为 rlog
ddsMat_rlog <- rlog(ddsMat, blind = FALSE)

# 按列变量绘制 PCA
plotPCA(ddsMat_rlog, intgroup = "Group", ntop = 500) +
  theme_bw() +
  geom_point(size = 5) + 
  scale_y_continuous(limits = c(-5, 5)) +
  ggtitle(label = "Principal Component Analysis (PCA)", 
          subtitle = "Top 500 most variable genes") 
plotPCA
plotPCA

9.2. Heatmap

# 将所有样本转换为 rlog
ddsMat_rlog <- rlog(ddsMat, blind = FALSE)

# 收集30个显著基因，制作矩阵
mat <- assay(ddsMat_rlog[row.names(results_sig)])[1:40, ]

# 选择您要用来注释列的列变量。
annotation_col = data.frame(
  Group = factor(colData(ddsMat_rlog)$Group), 
  Replicate = factor(colData(ddsMat_rlog)$Replicate),
  row.names = colData(ddsMat_rlog)$sampleid
)

# 指定要用来注释列的颜色。
ann_colors = list(
  Group = c(LoGlu = "lightblue", HiGlu = "darkorange"),
  Replicate = c(Rep1 = "darkred", Rep2 = "forestgreen")
)

# 使用 pheatmap 功能制作热图。
pheatmap(mat = mat, 
         color = colorRampPalette(brewer.pal(9, "YlOrBr"))(255), 
         scale = "row",
         annotation_col = annotation_col, 
         annotation_colors = ann_colors,
         fontsize = 6.5, 
         cellwidth = 55,
         show_colnames = F)
pheatmap
pheatmap

9.3. Volcano

# 从 DESeq2 结果中收集倍数变化和 FDR 校正的 pvalue
## - 将 pvalues 更改为 -log10 (1.3 = 0.05)
data <- data.frame(gene = row.names(results),
                   pval = -log10(results$padj), 
                   lfc = results$log2FoldChange)

# 删除任何以 NA 的行
data <- na.omit(data)

## If fold-change > 0 and pvalue > 1.3 (Increased significant)
## If fold-change < 0 and pvalue > 1.3 (Decreased significant)
data <- mutate(data, color = case_when(data$lfc > 0 & data$pval > 1.3 ~ "Increased",
                                       data$lfc < 0 & data$pval > 1.3 ~ "Decreased",
                                       data$pval < 1.3 ~ "nonsignificant"))

# 用 x-y 值制作一个基本的 ggplot2 对象
vol <- ggplot(data, aes(x = lfc, y = pval, color = color))

# 添加 ggplot2 图层
vol +   
  ggtitle(label = "Volcano Plot", subtitle = "Colored by fold-change direction") +
  geom_point(size = 2.5, alpha = 0.8, na.rm = T) +
  scale_color_manual(name = "Directionality",
                     values = c(Increased = "#008B00", Decreased = "#CD4F39", nonsignificant = "darkgray")) +
  theme_bw(base_size = 14) + 
  theme(legend.position = "right") + 
  xlab(expression(log[2]("LoGlu" / "HiGlu"))) + 
  ylab(expression(-log[10]("adjusted p-value"))) + 
  geom_hline(yintercept = 1.3, colour = "darkgrey") + 
  scale_y_continuous(trans = "log1p") 
Volcano
Volcano

9.4. MA

plotMA(results, ylim = c(-5, 5))
MA
MA

9.5. Dispersions

plotDispEsts(ddsMat)
plotDispEsts
plotDispEsts

9.6. 单基因图

# 将所有样本转换为 rlog
ddsMat_rlog <- rlog(ddsMat, blind = FALSE)

# 获得最高表达的基因
top_gene <- rownames(results)[which.min(results$log2FoldChange)]

# 画单基因图
plotCounts(dds = ddsMat, 
           gene = top_gene, 
           intgroup = "Group", 
           normalized = T, 
           transform = T)
单基因图
单基因图

10. 通路富集

  • 从差异表达基因中寻找通路

通路富集分析是基于单个基因变化生成结论的好方法。有时个体基因的变化是难以解释。但是通过分析基因的通路，我们可以收集基因反应的视图。

设置矩阵以考虑每个基因的 EntrezID 和倍数变化

# 删除没有任何 entrez 标识符的基因
results_sig_entrez <- subset(results_sig, is.na(entrez) == FALSE)

# 创建一个log2倍数变化的基因矩阵
gene_matrix <- results_sig_entrez$log2FoldChange

# 添加 entrezID 作为每个 logFC 条目的名称
names(gene_matrix) <- results_sig_entrez$entrez

# 查看基因矩阵的格式
##- Names = ENTREZ ID
##- Values = Log2 Fold changes
head(gene_matrix)  # 如下图
gene_matrix
gene_matrix

10.1. KEGG

  • 使用 KEGG 数据库丰富基因 
kegg_enrich <- enrichKEGG(gene = names(gene_matrix),
                          organism = 'mouse',
                          pvalueCutoff = 0.05, 
                          qvalueCutoff = 0.10)

# 结果可视化
barplot(kegg_enrich, 
        drop = TRUE, 
        showCategory = 10, 
        title = "KEGG Enrichment Pathways",
        font.size = 8)
KEGG
KEGG

10.2. GO

  • 使用 Gene Onotology 丰富基因 
go_enrich <- enrichGO(gene = names(gene_matrix),
                      OrgDb = 'org.Mm.eg.db', 
                      readable = T,
                      ont = "BP",
                      pvalueCutoff = 0.05, 
                      qvalueCutoff = 0.10)

# 结果可视化
barplot(go_enrich, 
        drop = TRUE, 
        showCategory = 10, 
        title = "GO Biological Pathways",
        font.size = 8)
GO
GO

11. 通路可视化

Pathview 是一个包，它可以获取显著差异表达基因的 KEGG 标识符，还可以与 KEGG 数据库中发现的其他生物一起使用，并且可以绘制特定生物的任何 KEGG 途径。

# 加载包
biocLite("pathview") ; library(pathview)

# 可视化通路 (用 fold change) 
## pathway.id : KEGG pathway identifier
pathview(gene.data = gene_matrix, 
         pathway.id = "04070", 
         species = "mouse")
pathview
pathview

欢迎Star -> 学习目录

国内链接 -> 学习目录

参考资料

[1]

Source: https://github.com/twbattaglia/RNAseq-workflow

本文由 mdnice 多平台发布

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

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

相关文章

计算机网络——数据链路层

计算机网络——数据链路层

数据链路层概述 数据链路层在网络体系结构中所处的地位 主机H1给主机H2发送数据时&#xff0c;需要通过路由器R1通过广域网链路转发到路由器R2&#xff0c;R2转发到主机H2。路由器转发只用到网络层及以下各层。【以上涉及数据包按网络体系结果逐层封装解封】 为了简单起见&am…
阅读更多...
DevC++的调试方法

DevC++的调试方法

目录 Dev C调试程序 Dev C调试注意事项对于修改后的程序&#xff0c;调试程序之前一定要先编译程序。 要想学会编程,第一步就是要学会调试(想我这种码龄一年的人还不会调试,丢死人). 今天,为了让你们的脸丢少点,特意写了这篇博文,给予需要帮助的人. 所谓调试程序&#xff0…
阅读更多...
[附源码]计算机毕业设计JAVA基于JSP健身房管理系统

[附源码]计算机毕业设计JAVA基于JSP健身房管理系统

[附源码]计算机毕业设计JAVA基于JSP健身房管理系统 项目运行 环境配置&#xff1a; Jdk1.8 Tomcat7.0 Mysql HBuilderX&#xff08;Webstorm也行&#xff09; Eclispe&#xff08;IntelliJ IDEA,Eclispe,MyEclispe,Sts都支持&#xff09;。 项目技术&#xff1a; SSM m…
阅读更多...
Python实现KNN(K近邻)分类模型(KNeighborsClassifier算法)并应用网格搜索算法寻找最优参数值以及数据标准化均衡化项目实战

Python实现KNN(K近邻)分类模型(KNeighborsClassifier算法)并应用网格搜索算法寻找最优参数值以及数据标准化均衡化项目实战

说明&#xff1a;这是一个机器学习实战项目&#xff08;附带数据代码文档视频讲解&#xff09;&#xff0c;如需数据代码文档视频讲解可以直接到文章最后获取。 1.项目背景 股票市场是已经发行的股票转让、买卖和流通的场所&#xff0c;包括交易所市场和场外交易市场两大类别。…
阅读更多...
【C++】C++11部分特性

【C++】C++11部分特性

目录 一、初始化列表 二、变量类型的推导 1、auto 2、decltype 三、右值引用 1、左值与右值 2、关于左值引用、右值引用的问题 1、左值引用可以引用右值吗&#xff1f; 2、右值引用可以引用左值吗&#xff1f; 3、右值引用之后的问题 3、移动构造、移动拷贝 1、引用…
阅读更多...
软件测试质量保证与测试

软件测试质量保证与测试

软件测试质量保证与测试 第一章 软件测试概述 1.1 软件测试背景 随着计算机技术的迅速发展和越来越广泛深入地应用于国民经济与社会生活的各个方面&#xff0c;软件系统的规模和复杂性与日俱增&#xff0c;软件的生产成本和软件中存在的缺陷与故障造成的各类损失也大大增加&…
阅读更多...
【应用】PostgreSQL 流复制配置

【应用】PostgreSQL 流复制配置

PostgreSQL 流复制配置centos7 安装 postgresql时序库 timescaleDB 的安装postgresql-14 主从流复制主库配置从库配置同步流复制与异步流复制异步流复制转换为同步流复制流复制的相关参数主从流复制原理PostgreSQL WAL 日志主从流复制架构主从流复制的过程基于 docker swarm 的…
阅读更多...
CMSC5713-IT项目管理之七、质量管理Quality Management

CMSC5713-IT项目管理之七、质量管理Quality Management

文章目录7.1 Quality7.2. Software Quality7.2.1. ISO/IEC 25010 Software Qualities7.2.2. Internal versus External qualities7.2.3. Software Metrics7.3. Quality Specification7.4. Project Quality Management7.4.1. Quality Planning7.4.2. Quality Assurance7.4.2.1. …
阅读更多...
vscode 阅读 linux kernel 源码

vscode 阅读 linux kernel 源码

前言 虽然身边的朋友大都在使用 source insight&#xff0c;但我却更喜欢 vscode。 不过 vscode 在代码搜索上确实不如 source insight&#xff0c;这点上我也是吃过亏的。阅读大型代码时&#xff0c;常常搜索不到关键代码&#xff0c;导致对代码的理解不充分。 当使用 vscode…
阅读更多...
Java-反射

Java-反射

前言 动态语言与静态语言 动态语言 是一类在运行时可以改变其结构的语言&#xff1a;例如新的函数、对象、甚至代码可以被引进&#xff0c;已有的函数可以被删除或是其他结构上的变化。通俗点说就是在运行时代码可以根据某些条件改变自身结构主要动态语言有&#xff1a;Object…
阅读更多...
【开源电路】STM32F401RCT6开发板

【开源电路】STM32F401RCT6开发板

【开源电路】STM32F401RCT6开发板&#x1f337;实物PCBA&#xff1a; &#x1f33c;优化后的3D效果图 &#x1f4da;STM32F401RCT6开发板简介 &#x1f4d1;主控是LQFP-64封装的STM32F401RCT6芯片&#xff0c;Micro USB接口供电&#xff0c;功能引脚全部引出&#xff0c;一个…
阅读更多...
金融强化学习与finRL开发包

金融强化学习与finRL开发包

原创文章第110篇&#xff0c;专注“个人成长与财富自由、世界运作的逻辑&#xff0c; AI量化投资”。 01 一些感受 时代的一粒沙&#xff0c;落在每个人身上就是一座山。 这三年&#xff0c;对于这句话&#xff0c;相信很多人更能感同身受。 看历史风云变幻&#xff0c;轻轻…
阅读更多...
力扣(LeetCode)21. 合并两个有序链表(C++)

力扣(LeetCode)21. 合并两个有序链表(C++)

迭代 同时遍历两个链表 &#xff0c; 当前结点值较小的结点插入新的链表尾部。直到有一个链表为空 &#xff0c; 我们将另一个非空链表插入新的链表尾部。 提示 : 使用哑结点&#xff0c;避免特判头结点。二路归并思想应用于链表~ class Solution { public:ListNode* mergeT…
阅读更多...
gRPC gateway - Http Restful API gRPC 协议转换

gRPC gateway - Http Restful API gRPC 协议转换

gRPC gateway - http restful gRPC gateway 介绍 gRPC-Gateway 是protocalBufffer的编译插件,根据protobuf的服务定义自动生成反向代理服务器&#xff0c;并将Restful Http API 转化为 gRPC通信协议。反向代理服务器根据 google.api.http 注解生成。 gRPC底层是使用HTTP2 协…
阅读更多...
mybatis # $

mybatis # $

总结&#xff1a; # 你传入的变量类型会被保留 $ 本质就是拼接 不会考虑拼串的 $ 情况下 参数是整数 跟参数是字符串 字符串情况&#xff1a; 缺少’’ 相当于字符串拼接 ”select * from t_user where username “ “张三” ”select * from t_user where username 张三&…
阅读更多...
java 容器

java 容器

java 容器 数组 数组的扩容问题 ArrayList 的默认初始化容量为0&#xff0c;首次添加元素时&#xff0c;创建容量为&#xff08;10 || 添加集合大小) ,以后每次扩容的话&#xff0c;为当前容量的1.5倍 public ArrayList() {/*初始化容量大小为0private static final Object…
阅读更多...
CEAC之《计算机应用助理工程师》2

CEAC之《计算机应用助理工程师》2

&#x1f468;‍&#x1f4bb;个人主页&#xff1a;微微的猪食小窝 欢迎 点赞&#x1f44d; 收藏⭐ 留言&#x1f4dd; 加关注✅! 本文由 微微的猪食小窝 原创 收录于专栏 【CEAC证书】 1组合框的常用属性有 ____________ 。 A、Index B、Text C、Caption D、ListCountA,B,D2在…
阅读更多...
ES6 入门教程 16 Reflect 16.2 静态方法 16.3 实例：使用 Proxy 实现观察者模式

ES6 入门教程 16 Reflect 16.2 静态方法 16.3 实例：使用 Proxy 实现观察者模式

ES6 入门教程 ECMAScript 6 入门 作者&#xff1a;阮一峰 本文仅用于学习记录&#xff0c;不存在任何商业用途&#xff0c;如侵删 文章目录ES6 入门教程16 Reflect16.2 静态方法16.2.1 Reflect.get(target, name, receiver)16.2.2 Reflect.set(target, name, value, receiver)1…
阅读更多...
数据结构之：数组

数据结构之：数组

数组初体验之数组中重复的数字 数组 &#xff1a; 有限个 相同类型的变量 组成的有序集合 int[] arr; int arr[]; // 静态初始化 String[] strArr {"和平精英","王者荣耀","开心消消乐","欢乐斗地主"}; // 动态初始化 String[] strAr…
阅读更多...
自学 TypeScript 第三天 使用webpack打包 TS 代码

自学 TypeScript 第三天 使用webpack打包 TS 代码

安装&#xff1a; 首先第一步&#xff0c;我们要初始化我们项目&#xff0c;在目录下输入 npm init 接下来&#xff0c;我们的安装几个工具 npm i -D webpack webpack-cli typescript ts-loader -D 意思是 开发依赖&#xff0c;也就是我们现在所安装的依赖都是开发依赖&am…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023