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

基因表达，漂亮火山图

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

绘制超漂亮的基因差异表达火山图相关推荐

sangerbox平台使用（三）绘制火山图
目前绘制火山图有很多工具,我们最常用的就是R语言,本次想向大家介绍的是一个在线工具直接绘制火山图-sangerbox 官方链接为:http://sangerbox.com/AllTools?tool_ ...
R语言绘制火山图（volcano plot）实战：为差异表达基因（DEGs）添加颜色、基于显著性阈值进行点的颜色美化、为选定基因添加标签
R语言绘制火山图(volcano plot)实战:为差异表达基因(DEGs)添加颜色.基于显著性阈值进行点的颜色美化.为选定基因添加标签目录 R语言绘制火山图(volcano plot)实战 #导入 ...
ggplot2绘制差异表达基因火山图
一.前置环境 1.1 R 语言下载对应系统的R软件 R: The R Project for Statistical Computing (r-project.org) 以win11为演示 http ...
差异表达基因-火山图和聚类图解释
想研究某现象的分子机制,老板豪气的来一句,先测个转录组吧,看下差异表达基因. 是否在心里窃喜,制个样就完事了,太easy有木有.等大堆数据回来的时候,是不是傻眼了? 从何下手挑选差异表达基因呢? 今天 ...
差异表达基因火山图（ggplot函数）
1. 读入数据差异表达基因来自limma分析结果. # read the file data <- read.csv("diff_expr_genes.csv",row.n ...
BIC无代码绘制差异基因火山图
无代码绘制差异基因火山图 Volcano plot | 别再问我这为什么是火山图一文解释了火山图如何解读.不太难看懂,而一旦看懂了,图也就知道怎么绘制了. 假设我们已经有了一个差异基因鉴定后的表格文 ...
ImageGP/BIC无代码绘制差异基因火山图
无代码绘制差异基因火山图 Volcano plot | 别再问我这为什么是火山图一文解释了火山图如何解读.不太难看懂,而一旦看懂了,图也就知道怎么绘制了. 假设我们已经有了一个差异基因鉴定后的表格文 ...
再肝一个R包！一行代码绘制精美火山图！
一行代码绘制火山图的R包诞生了!在过去的一年中,师兄先后生信绘图系列和高分SCI复现系列中更新了多种不同的火山图的绘制方法,包括普通的火山图.渐变火山图.以及包含GO通路信息的火山图!但是很多小伙伴反 ...
R绘图笔记 | 火山图的绘制
参考前文:R绘图笔记 | R语言绘图系统与常见绘图函数及参数关于绘图,前面介绍了一些: R绘图笔记 | 一般的散点图绘制 R绘图笔记 | 柱状图绘制 R绘图笔记 | 直方图和核密度估计图的绘制 R绘 ...

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

绘制超漂亮的基因差异表达火山图相关推荐

最新文章

热门文章