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

基因表达，漂亮火山图

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. 0. 0. 2. 0.0 0. # ENSG00000000005 0.0000000 NA NA NA NA NA # ENSG00000000419 520. -0.0 0. -0. 0. 0. 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('ENSG00000', 'ENSG00000002933', 'ENSG00000', 'ENSG00000') celltype2 <- c('ENSG00000', 'ENSG00000', 'ENSG00000', 'ENSG00000', 'ENSG00000', 'ENSG00000') # 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))