R语言可视化5: 火山图 - EnhancedVolcano

作者: cc的生信随记 | 来源:发表于2023-03-05 00:00 被阅读0次

R可视化：使用EnhancedVolcano包可视化火山图
R语言之可视化①④一页多图（1）
R语言之可视化⑥R图形系统续
R语言之可视化⑩坐标系统
R语言之可视化⑨火山图
R语言之可视化⑦easyGgplot2散点图
R语言之可视化⑧easyGgplot2散点图续
R语言之可视化①③散点图+拟合曲线
R语言之可视化①⑥一页多图（2）
R语言之可视化④韦恩图upsetR

1. 使用 $\color{green}{EnhancedVolcano}$ 包绘制火山图

1.1 基本用法

# 安装并加载所需的R包
# if (!requireNamespace('BiocManager', quietly = TRUE))
#    install.packages('BiocManager')
#  BiocManager::install('EnhancedVolcano')

# devtools::install_github('kevinblighe/EnhancedVolcano')
library(EnhancedVolcano)

# 加载数据
library(airway)
library(magrittr)

data('airway') # 加载airway数据，其中不同的气道平滑肌细胞用地塞米松治疗
airway$dex %<>% relevel('untrt') # %<>%复合赋值操作符， 功能与 %>% 基本是一样的，但多了一项额外的操作，就是把结果写到左侧对象

# 将 Ensembl 基因 ID 注释为基因symbol
library(org.Hs.eg.db)
symbols <- mapIds(org.Hs.eg.db, keys = ens,
                  column = c('SYMBOL'), keytype = 'ENSEMBL')
symbols <- symbols[!is.na(symbols)]
symbols <- symbols[match(rownames(airway), names(symbols))]
rownames(airway) <- symbols
keep <- !is.na(rownames(airway))
airway <- airway[keep,]

# 使用DESeq2进行差异表达，以创建两组结果
library('DESeq2')

dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- DESeq(dds, betaPrior = FALSE)
res <- results(dds,
               contrast = c('dex', 'trt', 'untrt'))
res <- lfcShrink(dds,
                 contrast = c('dex', 'trt', 'untrt'), res = res, type  =  'normal')

# 绘制基本的火山图，log2FC 的默认截止值是 >|2|；P 值的默认截止值为 10e-6
EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue')

EnhancedVolcano-1

1.2 其他参数

EnhancedVolcano(res,
                lab = rownames(res),
                x = 'log2FoldChange',
                y = 'pvalue',
                xlim = c(-6, 6), # 修改横坐标范围
                title = 'N061011 versus N61311', # 更改标题
                pCutoff = 10e-16, # 更改P值
                FCcutoff = 1.5, # 更改log2FC值
                pointSize = 3.0, # 点的大小
                labSize = 6.0, # 标题字体的大小
                shape = c(1, 4, 23, 25), # 不同的点绘制成不同的形状
                col=c('black', 'green', 'blue', 'red3'), # 修改配色
                colAlpha = 1, # 控制绘制点的透明度：1 = 100% 不透明；0 = 100% 透明
                cutoffLineType = 'solid', # FCcutoff 和 pCutoff 的线型，有：'blank'、'solid'、'dashed'、'dotted'、'dotdash'、'longdash'、'twodash'
                cutoffLineCol = 'black', # FCcutoff 和 pCutoff 的线条颜色
                cutoffLineWidth = 0.8, # FCcutoff 和 pCutoff 的线宽
                hline = c(10e-20,
                          10e-20 * 10e-30,
                          10e-20 * 10e-60,
                          10e-20 * 10e-90), # 在y轴上绘制一条或多条穿过此/这些值的水平线
                hlineCol = c('pink', 'hotpink', 'purple', 'black'), # hline 的颜色
                hlineType = c('solid', 'longdash', 'dotdash', 'dotted'), # hline 的线型，有：'blank'、'solid'、'dashed'、'dotted'、'dotdash'、'longdash'、'twodash'
                hlineWidth = c(1.0, 1.5, 2.0, 2.5), # hline 的线宽
                gridlines.major = FALSE, # 删除绘制主要网格线
                gridlines.minor = TRUE, # 删除次网格线
                legendLabels=c('Not sig.','Log (base 2) FC','p-value',
                               'p-value & Log (base 2) FC'), # 修改图例的文字
                legendPosition = 'right', # 修改图例的位置
                legendLabSize = 16, # 修改图例的文本大小
                legendIconSize = 5.0, # 图例图标/符号的大小
                drawConnectors = TRUE,
                selectLab = c('VCAM1','KCTD12','ADAM12',
                              'CXCL12','CACNB2','SPARCL1','DUSP1','SAMHD1','MAOA'), # 标记感兴趣的关键变量/变量
                widthConnectors = 0.75, # 连接器（箭头）的线宽
                boxedLabels = TRUE) # 在标签周围绘制简单的框

EnhancedVolcano-2

1.3 斜体标签并在其一侧翻转火山

# 创建一个新的向量，并对标签进行设置
lab_italics <- paste0("italic('", rownames(res), "')")
  selectLab_italics = paste0(
    "italic('",
    c('VCAM1','KCTD12','ADAM12', 'CXCL12','CACNB2','SPARCL1','DUSP1','SAMHD1','MAOA'),
    "')")

EnhancedVolcano(res,
    lab = lab_italics,
    x = 'log2FoldChange',
    y = 'pvalue',
    selectLab = selectLab_italics,
    xlab = bquote(~Log[2]~ 'fold change'),
    pCutoff = 10e-14,
    FCcutoff = 1.0,
    pointSize = 3.0,
    labSize = 6.0,
    labCol = 'black',
    labFace = 'bold',
    boxedLabels = TRUE,
    parseLabels = TRUE,
    col = c('black', 'pink', 'purple', 'red3'),
    colAlpha = 4/5,
    legendPosition = 'bottom',
    legendLabSize = 14,
    legendIconSize = 4.0,
    drawConnectors = TRUE,
    widthConnectors = 1.0,
    colConnectors = 'black') + coord_flip() # 要翻转火山，需使用EnhancedVolcano(…) + coord_flip()

EnhancedVolcano-3

1.4 自定义方案

# 定义将被着色的不同细胞类型
celltype1 <- c('VCAM1','KCTD12','ADAM12','CXCL12')
celltype2 <- c('CACNB2','SPARCL1','DUSP1','SAMHD1','MAOA')

# 为不同的单元格类型创建自定义键值对
keyvals.shape <- ifelse(
  rownames(res) %in% celltype1, 17,
  ifelse(rownames(res) %in% celltype2, 64,
         3))
keyvals.shape[is.na(keyvals.shape)] <- 3
names(keyvals.shape)[keyvals.shape == 3] <- 'PBMC'
names(keyvals.shape)[keyvals.shape == 17] <- 'Cell-type 1'
names(keyvals.shape)[keyvals.shape == 64] <- 'Cell-type 2'

p1 <- EnhancedVolcano(res,
                      lab = rownames(res),
                      x = 'log2FoldChange',
                      y = 'pvalue',
                      selectLab = rownames(res)[which(names(keyvals) %in% c('high', 'low'))],
                      xlab = bquote(~Log[2]~ 'fold change'),
                      title = 'Custom shape over-ride',
                      pCutoff = 10e-14,
                      FCcutoff = 1.0,
                      pointSize = 4.5,
                      labSize = 4.5,
                      shapeCustom = keyvals.shape, # 用自己的形状方案覆盖默认形状
                      colCustom = NULL,
                      colAlpha = 1,
                      legendLabSize = 15,
                      legendPosition = 'left',
                      legendIconSize = 5.0,
                      drawConnectors = TRUE,
                      widthConnectors = 0.5,
                      colConnectors = 'grey50',
                      gridlines.major = TRUE,
                      gridlines.minor = FALSE,
                      border = 'partial',
                      borderWidth = 1.5,
                      borderColour = 'black')

# 创建自定义键值对，通过fold-change分为'high', 'low', 'mid'
keyvals.colour <- ifelse(
  res$log2FoldChange < -2.5, 'royalblue',
  ifelse(res$log2FoldChange > 2.5, 'gold',
         'black'))
keyvals.colour[is.na(keyvals.colour)] <- 'black'
names(keyvals.colour)[keyvals.colour == 'gold'] <- 'high'
names(keyvals.colour)[keyvals.colour == 'black'] <- 'mid'
names(keyvals.colour)[keyvals.colour == 'royalblue'] <- 'low'

p2 <- EnhancedVolcano(res,
                      lab = rownames(res),
                      x = 'log2FoldChange',
                      y = 'pvalue',
                      selectLab = rownames(res)[which(names(keyvals) %in% c('High', 'Low'))],
                      xlab = bquote(~Log[2]~ 'fold change'),
                      title = 'Custom shape & colour over-ride',
                      pCutoff = 10e-14,
                      FCcutoff = 1.0,
                      pointSize = 5.5,
                      labSize = 0.0,
                      shapeCustom = keyvals.shape, # 用自己的形状方案覆盖默认形状
                      colCustom = keyvals.colour,
                      colAlpha = 1,
                      legendPosition = 'right',
                      legendLabSize = 15,
                      legendIconSize = 5.0,
                      drawConnectors = TRUE,
                      widthConnectors = 0.5,
                      colConnectors = 'grey50',
                      gridlines.major = TRUE,
                      gridlines.minor = FALSE,
                      border = 'full',
                      borderWidth = 1.0,
                      borderColour = 'black')

library(gridExtra)
library(grid)
grid.arrange(p1, p2,
             ncol=2,
             top = textGrob('EnhancedVolcano',
                            just = c('center'),
                            gp = gpar(fontsize = 32)))

EnhancedVolcano-4

1.5 圈出/突出某些变量（此功能最适合仅对 1 或 2 个关键变量进行着色）

# 定义将被着色的不同细胞类型
  celltype1 <- c('VCAM1','CXCL12')
  celltype2 <- c('SORT1', 'KLF15')

EnhancedVolcano(res,
    lab = rownames(res),
    x = 'log2FoldChange',
    y = 'pvalue',
    selectLab = c(celltype1, celltype2),
    xlab = bquote(~Log[2]~ 'fold change'),
    title = 'Shading cell-type 1|2',
    pCutoff = 10e-14,
    FCcutoff = 1.0,
    pointSize = 8.0,
    labSize = 6.0,
    labCol = 'black',
    labFace = 'bold',
    boxedLabels = TRUE,
    shape = 42,
    colCustom = keyvals,
    colAlpha = 1,
    legendPosition = 'right',
    legendLabSize = 20,
    legendIconSize = 20.0,
    # encircle
      encircle = celltype1,
      encircleCol = 'black',
      encircleSize = 2.5,
      encircleFill = 'pink',
      encircleAlpha = 1/2,
    # shade
      shade = celltype2,
      shadeAlpha = 1/2,
      shadeFill = 'skyblue',
      shadeSize = 1,
      shadeBins = 5,
    drawConnectors = TRUE,
    widthConnectors = 2.0,
    gridlines.major = TRUE,
    gridlines.minor = FALSE,
    border = 'full',
    borderWidth = 5,
    borderColour = 'black')

EnhancedVolcano-5

1.6 自定义关键变量

library("pasilla")
  pasCts <- system.file("extdata", "pasilla_gene_counts.tsv",
    package="pasilla", mustWork=TRUE)
  pasAnno <- system.file("extdata", "pasilla_sample_annotation.csv",
    package="pasilla", mustWork=TRUE)
  cts <- as.matrix(read.csv(pasCts,sep="\t",row.names="gene_id"))
  coldata <- read.csv(pasAnno, row.names=1)
  coldata <- coldata[,c("condition","type")]
  rownames(coldata) <- sub("fb", "", rownames(coldata))
  cts <- cts[, rownames(coldata)]
  library("DESeq2")
  dds <- DESeqDataSetFromMatrix(countData = cts,
    colData = coldata,
    design = ~ condition)

  featureData <- data.frame(gene=rownames(cts))
  mcols(dds) <- DataFrame(mcols(dds), featureData)
  dds <- DESeq(dds)
  res <- results(dds)

p1 <- EnhancedVolcano(res,
    lab = rownames(res),
    x = "log2FoldChange",
    y = "pvalue",
    pCutoff = 10e-4,
    FCcutoff = 2,
    ylim = c(0, -log10(10e-12)),
    pointSize = c(ifelse(res$log2FoldChange>2, 8, 1)), # 只更改 log 2 FC>2 的那些变量的大小
    labSize = 6.0,
    shape = c(6, 6, 19, 16),
    title = "DESeq2 results",
    subtitle = "Differential expression",
    caption = bquote(~Log[2]~ "fold change cutoff, 2; p-value cutoff, 10e-4"),
    legendPosition = "right",
    legendLabSize = 14,
    col = c("grey30", "forestgreen", "royalblue", "red2"),
    colAlpha = 0.9,
    drawConnectors = TRUE,
    hline = c(10e-8),
    widthConnectors = 0.5)

p2 <- EnhancedVolcano(res,
    lab = rownames(res),
    x = "log2FoldChange",
    y = "pvalue",
    pCutoff = 10e-4,
    FCcutoff = 2,
    ylim = c(0, -log10(10e-12)),
    pointSize = c(ifelse(res$log2FoldChange>2, 8, 1)),
    labSize = 6.0,
    shape = c(6, 6, 19, 16),
    title = "DESeq2 results",
    subtitle = "Differential expression",
    caption = bquote(~Log[2]~ "fold change cutoff, 2; p-value cutoff, 10e-4"),
    legendPosition = "right",
    legendLabSize = 14,
    colAlpha = 0.9,
    colGradient = c('red3', 'royalblue'), # 更改为连续配色方案
    drawConnectors = TRUE,
    hline = c(10e-8),
    widthConnectors = 0.5)

grid.arrange(p1, p2,
    ncol=2,
    top = textGrob('EnhancedVolcano',
      just = c('center'),
      gp = gpar(fontsize = 32)))

EnhancedVolcano-6

参考：

https://github.com/kevinblighe/EnhancedVolcano

R可视化：使用EnhancedVolcano包可视化火山图
使用EnhancedVolcano包可视化火山图火山图是用来展示差异表达的基因或者物种，其常常出现RNA-seq...
R语言之可视化①④一页多图（1）
目录 R语言之可视化①误差棒 R语言之可视化②点图 R语言之可视化③点图续 R语言之可视化④点韦恩图upsetR ...
R语言之可视化⑥R图形系统续
目录 R语言之可视化①误差棒 R语言之可视化②点图 R语言之可视化③点图续 R语言之可视化④点韦恩图upsetR ...
R语言之可视化⑩坐标系统
目录 R语言之可视化①误差棒 R语言之可视化②点图 R语言之可视化③点图续 R语言之可视化④点韦恩图upsetR ...
R语言之可视化⑨火山图
目录 R语言之可视化①误差棒 R语言之可视化②点图 R语言之可视化③点图续 R语言之可视化④点韦恩图upsetR ...
R语言之可视化⑦easyGgplot2散点图
目录 R语言之可视化①误差棒 R语言之可视化②点图 R语言之可视化③点图续 R语言之可视化④点韦恩图upsetR ...
R语言之可视化⑧easyGgplot2散点图续
目录 R语言之可视化①误差棒 R语言之可视化②点图 R语言之可视化③点图续 R语言之可视化④点韦恩图upsetR ...
R语言之可视化①③散点图+拟合曲线
目录 R语言之可视化①误差棒 R语言之可视化②点图 R语言之可视化③点图续 R语言之可视化④点韦恩图upsetR ...
R语言之可视化①⑥一页多图（2）
目录 R语言之可视化①误差棒 R语言之可视化②点图 R语言之可视化③点图续 R语言之可视化④点韦恩图upsetR ...
R语言之可视化④韦恩图upsetR
目录 R语言之可视化①误差棒 R语言之可视化②点图 R语言之可视化③点图续 R语言之可视化④点韦恩图upsetR ...