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

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

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

数据准备

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

library(easyTCGA)
getmrnaexpr("TCGA-SKCM")

加载数据：

load(file = "G:/easyTCGA_test/output_mRNA_lncRNA_expr/TCGA-SKCM_mrna_expr_tpm.rdata")

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

expr <- log2(mrna_expr_tpm+1)
dim(expr)
## [1] 19938   473

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

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

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

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绘制火山图

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

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

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

# 要标记的基因
marker <- c("CSTA","AQP3","CD3D","CXCL9","CXCL13","CXCL10","S100A9","S100A8","JCHAIN","KRT5")

下面画图即可：

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()
        )

EnhancedVolcano

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

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
                )

pheatmap绘制热图

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

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

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)

有了数据，画图就很简单：

library(pheatmap)

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