前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >专栏 >单基因差异分析并绘制火山图和热图

单基因差异分析并绘制火山图和热图

作者头像
医学和生信笔记
发布2023-08-30 12:57:03
发布2023-08-30 12:57:03
77400
代码可运行
举报
运行总次数:0
代码可运行

介绍下绘制火山图和热图的方法,如何在火山图或者热图中标记特定的基因,顺便学习下EnhancedVolcano包绘制火山图。

前面已经介绍了单基因富集分析:单基因富集分析

数据准备

使用TCGA黑色素瘤的转录组数据,使用easyTCGA,1行代码下载,即可得到6种表达矩阵和临床信息,而且是官网最新的数据:

代码语言:javascript
代码运行次数:0
运行
复制
library(easyTCGA)
getmrnaexpr("TCGA-SKCM")

加载数据:

代码语言:javascript
代码运行次数:0
运行
复制
load(file = "G:/easyTCGA_test/output_mRNA_lncRNA_expr/TCGA-SKCM_mrna_expr_tpm.rdata")

这个数据是直接从GDC的官网数据中提取出来的,没有经过任何转化,所以我们先进行log2转换。

代码语言:javascript
代码运行次数:0
运行
复制
expr <- log2(mrna_expr_tpm+1)
dim(expr)
## [1] 19938   473

一共有19938个mRNA和473个样本。

我们这里根据HOPX表达量中位数进行分组,把所有样本分为高表达组和低表达组。

然后进行差异分析,这里也是用easyTCGA1行代码解决:

代码语言:javascript
代码运行次数:0
运行
复制
sample_group <- ifelse(expr["HOPX",] > median(t(expr["HOPX",])), "high","low")
sample_group <- factor(sample_group, levels = c("low","high"))

library(easyTCGA)
deg_res <- diff_analysis(exprset = expr, 
                         group = sample_group,
                         is_count = F,
                         save = F
                         )
## => log2 transform not needed
## => Running limma
## => Running wilcoxon test
## => Analysis done.

# 提取limma的结果
deg_limma <- deg_res$deg_limma

head(deg_limma)
##            logFC  AveExpr        t      P.Value    adj.P.Val         B
## HOPX    1.411491 1.330710 17.92159 3.648490e-55 7.274359e-51 114.22524
## IL18    1.579073 2.882566 15.19126 1.021080e-42 1.017914e-38  86.06859
## TNFSF10 1.627203 3.868731 14.82989 4.077339e-41 2.709799e-37  82.44504
## DAPP1   1.120535 1.551681 13.49540 2.489357e-35 1.240820e-31  69.35296
## ANKRD22 1.459943 1.975013 13.06324 1.661069e-33 6.623677e-30  65.22542
## ZBED2   1.241665 1.613864 12.12598 1.202387e-29 3.995531e-26  56.49488
##         genesymbol
## HOPX          HOPX
## IL18          IL18
## TNFSF10    TNFSF10
## DAPP1        DAPP1
## ANKRD22    ANKRD22
## ZBED2        ZBED2

ggplot2绘制火山图

绘制火山图需要差异分析的结果,我们再增加一列信息展示这个基因是上调、下调还是没意义。

代码语言:javascript
代码运行次数:0
运行
复制
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
library(ggrepel)

deg_limma <- deg_limma %>% 
  mutate(type = case_when(logFC > 1 & adj.P.Val < 0.05 ~ "Up-regulation",
                          logFC < -1 & adj.P.Val < 0.05 ~ "Down-regulation",
                          .default = "None"
                          ))
head(deg_limma)
##            logFC  AveExpr        t      P.Value    adj.P.Val         B
## HOPX    1.411491 1.330710 17.92159 3.648490e-55 7.274359e-51 114.22524
## IL18    1.579073 2.882566 15.19126 1.021080e-42 1.017914e-38  86.06859
## TNFSF10 1.627203 3.868731 14.82989 4.077339e-41 2.709799e-37  82.44504
## DAPP1   1.120535 1.551681 13.49540 2.489357e-35 1.240820e-31  69.35296
## ANKRD22 1.459943 1.975013 13.06324 1.661069e-33 6.623677e-30  65.22542
## ZBED2   1.241665 1.613864 12.12598 1.202387e-29 3.995531e-26  56.49488
##         genesymbol          type
## HOPX          HOPX Up-regulation
## IL18          IL18 Up-regulation
## TNFSF10    TNFSF10 Up-regulation
## DAPP1        DAPP1 Up-regulation
## ANKRD22    ANKRD22 Up-regulation
## ZBED2        ZBED2 Up-regulation

随便选择几个基因展示一下:

代码语言:javascript
代码运行次数:0
运行
复制
# 要标记的基因
marker <- c("CSTA","AQP3","CD3D","CXCL9","CXCL13","CXCL10","S100A9","S100A8","JCHAIN","KRT5")

下面画图即可:

代码语言:javascript
代码运行次数:0
运行
复制
ggplot(deg_limma, aes(logFC, -log10(adj.P.Val)))+
  geom_point(aes(color=type))+
  scale_color_manual(values = c("blue","gray70","red"),name = NULL)+
  geom_hline(yintercept = -log10(0.05),linetype=2)+
  geom_vline(xintercept = c(-1,1), linetype=2)+
  geom_label_repel(data = subset(deg_limma, genesymbol %in% marker), 
                  aes(label=genesymbol),col="black",alpha = 0.8)+
  theme_bw()+
  theme(legend.position = c(0.15,0.8),
        legend.background = element_blank()
        )

plot of chunk unnamed-chunk-7

EnhancedVolcano

一个专门用于绘制火山图的R包,增加了很多好用的功能,比如添加要展示的基因更加方便:

代码语言:javascript
代码运行次数:0
运行
复制
library(EnhancedVolcano)

EnhancedVolcano(deg_limma,
                x = "logFC",
                y = "adj.P.Val",
                lab = rownames(deg_limma),
                selectLab = marker,
                boxedLabels = TRUE,
                drawConnectors = TRUE,
                pCutoff = -log10(0.05),
                FCcutoff = 1,
                title = NULL,
                subtitle = NULL
                )

plot of chunk unnamed-chunk-8

pheatmap绘制热图

首先是准备热图需要的数据,其实就是表达矩阵的可视化而已。

我们选择311个差异基因的表达矩阵进行展示。

代码语言:javascript
代码运行次数:0
运行
复制
deg_genes <- deg_limma %>% 
  filter(abs(logFC)>1, adj.P.Val < 0.05) %>% 
  pull(genesymbol)

heat_df <- expr[deg_genes,]

# 自己进行标准化的好处就是可以自定义范围
heat_df <- t(scale(t(heat_df)))
heat_df[heat_df >= 3] <- 3
heat_df[heat_df < -3] <- -3

# 列注释
anno_col <- data.frame(sample_group = sample_group)
rownames(anno_col) <- colnames(heat_df)

有了数据,画图就很简单:

代码语言:javascript
代码运行次数:0
运行
复制
library(pheatmap)

pheatmap::pheatmap(heat_df,
         annotation_col = anno_col,
         show_colnames = F,
         show_rownames = F
         )

plot of chunk unnamed-chunk-10

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2023-07-10,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 医学和生信笔记 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 数据准备
  • ggplot2绘制火山图
  • EnhancedVolcano
  • pheatmap绘制热图
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档