绘制超漂亮的基因差异表达火山图

基因表达，漂亮火山图

rm(list = ls())
getwd()
#ctrl+shift+h切换工作目录
library(EnhancedVolcano)
library(DESeq2)
library(airway)
library(magrittr)
data('airway')
# %<>%复合赋值操作符， 功能与 %>% 基本是一样的，但多了一项额外的操作，就是把结果写到左侧对象。
# 对dex列进行relevel，再把revel后的结果赋值到airway$dex。
airway$dex %<>% relevel('untrt')
#使用DESeq2进行差异表达，以创建两组结果(DESeq2差异基因分析和批次效应移除)
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds <- DESeq(dds, betaPrior=FALSE)
# estimating size factors
# estimating dispersions
# gene-wise dispersion estimates
# mean-dispersion relationship
# final dispersion estimates
# fitting model and testing

# compare trt & untrt
res1 <- results(dds,contrast = c('dex','trt','untrt'))
head(res1)
# shrink(收缩) log2 fold change
res1 <- lfcShrink(dds,contrast = c('dex','trt','untrt'), res=res1)
head(res1)
# compare different cells
res2 <- results(dds,contrast = c('cell', 'N061011', 'N61311'))
res2 <- lfcShrink(dds,contrast = c('cell', 'N061011', 'N61311'), res=res2)
head(res2)
head(as.data.frame(res2))
#                    baseMean log2FoldChange      lfcSE        stat     pvalue      padj
# ENSG00000000003 708.6021697     0.29018323 0.13601829  2.13421942 0.03282482 0.2201905
# ENSG00000000005   0.0000000             NA         NA          NA         NA        NA
# ENSG00000000419 520.2979006    -0.05060086 0.14954177 -0.33839235 0.73506754 0.9390635
##saving result
result <- as.data.frame(res2)
result$Gene <- rownames(result)
write.table(result,file = "N061011_vs_N61311_difference.txt",sep = "\t",row.names = FALSE,quote =FALSE)
#######
inputFileName <- "N061011_vs_N61311_difference.txt"
# Import DESeq2 data
volcanodata <- read.table(inputFileName, header=TRUE, sep="\t",stringsAsFactors = FALSE)
str(volcanodata)
colnames(volcanodata)
# [1] "baseMean"       "log2FoldChange" "lfcSE"          "stat"           "pvalue"         "padj"          
# [7] "Gene"  
# Set transcript names as row names
row.names(volcanodata)<-make.names(volcanodata[,"Gene"], unique=TRUE)
p <- EnhancedVolcano(volcanodata,
                lab = rownames(volcanodata),
                x = "log2FoldChange",
                y = "pvalue",
                xlim = c(-8, 8),
                #ylim = c(0,6.1),
                title = 'N061011 versus N61311',
                pCutoff = 10e-16,
                FCcutoff = 1.5,
                transcriptPointSize = 1.5,
                transcriptLabSize = 3.0)
p
ggsave(p,file="N061011_vs_N61311_EnhancedVolcano_plot.pdf",width = 9,height = 9)
#############################################
# create custom key-value pairs for 'high', 'low', 'mid' expression by fold-change
# 通过named vector生成自定义颜色
# set the base colour as 'black'
keyvals <- rep('black', nrow(res2))

# set the base name/label as 'Mid'
names(keyvals) <- rep('Mid', nrow(res2))

# modify keyvals for transcripts with fold change > 2.5
keyvals[which(res2$log2FoldChange > 2.5)] <- 'gold'
names(keyvals)[which(res2$log2FoldChange > 2.5)] <- 'high'

# modify keyvals for transcripts with fold change < -2.5
keyvals[which(res2$log2FoldChange < -2.5)] <- 'royalblue'
names(keyvals)[which(res2$log2FoldChange < -2.5)] <- 'low'

unique(names(keyvals))
# define different cell-types that will be shaded
celltype1 <- c('ENSG00000106565', 'ENSG00000002933',
               'ENSG00000165246', 'ENSG00000224114')
celltype2 <- c('ENSG00000230795', 'ENSG00000164530',
               'ENSG00000143153', 'ENSG00000169851',
               'ENSG00000231924', 'ENSG00000145681')

# create custom key-value pairs for different cell-types
# set the base shape as '3'
keyvals.shape <- rep(3, nrow(res2))

# set the base name/label as 'PBC'
names(keyvals.shape) <- rep('PBC', nrow(res2))

# modify the keyvals for cell-type 1
keyvals.shape[which(rownames(res2) %in% celltype1)] <- 17
names(keyvals.shape)[which(rownames(res2) %in% celltype1)] <- 'Cell-type 1'

# modify the keyvals for cell-type 2
keyvals.shape[which(rownames(res2) %in% celltype2)] <- 64
names(keyvals.shape)[which(rownames(res2) %in% celltype2)] <- 'Cell-type 2'

unique(names(keyvals.shape))
p1 <- EnhancedVolcano(res2,
                      lab = rownames(res2),
                      x = 'log2FoldChange',
                      y = 'pvalue',
                      selectLab = rownames(res2)[which(names(keyvals) %in% c('high', 'low'))],
                      xlim = c(-6.5,6.5),
                      xlab = bquote(~Log[2]~ 'fold change'),
                      title = 'Custom shape over-ride',
                      pCutoff = 10e-14,
                      FCcutoff = 1.0,
                      transcriptPointSize = 3.5,
                      transcriptLabSize = 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')

# create custom key-value pairs for 'high', 'low', 'mid' expression by fold-change
# set the base colour as 'black'
keyvals.colour <- rep('black', nrow(res2))

# set the base name/label as 'Mid'
names(keyvals.colour) <- rep('Mid', nrow(res2))

# modify keyvals for transcripts with fold change > 2.5
keyvals.colour[which(res2$log2FoldChange > 2.5)] <- 'gold'
names(keyvals.colour)[which(res2$log2FoldChange > 2.5)] <- 'high'

# modify keyvals for transcripts with fold change < -2.5
keyvals.colour[which(res2$log2FoldChange < -2.5)] <- 'royalblue'
names(keyvals.colour)[which(res2$log2FoldChange < -2.5)] <- 'low'

unique(names(keyvals.colour))

p2 <- EnhancedVolcano(res2,
                      lab = rownames(res2),
                      x = 'log2FoldChange',
                      y = 'pvalue',
                      selectLab = rownames(res2)[which(names(keyvals) %in% c('High', 'Low'))],
                      xlim = c(-6.5,6.5),
                      xlab = bquote(~Log[2]~ 'fold change'),
                      title = 'Custom shape & colour over-ride',
                      pCutoff = 10e-14,
                      FCcutoff = 1.0,
                      transcriptPointSize = 5.5,
                      transcriptLabSize = 0.0,
                      shapeCustom = keyvals.shape,
                      colCustom = keyvals.colour,
                      colAlpha = 1,
                      legendPosition = 'top',
                      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)))
grid.rect(gp=gpar(fill=NA))