首页 > 美文阅读

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

更新时间:2023-06-30 20:48:19 阅读：评论：0

绘制超漂亮的基因差异表达⽕⼭图基因表达，漂亮⽕⼭图

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

# baMean 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_",p = "\t",row.names = FALSE,quote =FALSE) >##

inputFileName <- "N061011_vs_"

# Import DESeq2 data

volcanodata <- read.table(inputFileName, header=TRUE, p="\t",stringsAsFactors = FALSE)

str(volcanodata)

colnames(volcanodata)

# [1] "baMean" "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)

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⽣成⾃定义颜⾊

# t the ba colour as 'black'

keyvals <- rep('black', nrow(res2))

# t the ba 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))

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

# t the ba shape as '3'

keyvals.shape <- rep(3, nrow(res2))

# t the ba 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',

lectLab = 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

# t the ba colour as 'black'

# t the ba name/label as 'Mid'

lour) <- rep('Mid', nrow(res2))

# modify keyvals for transcripts with fold change > 2.5

lour)[which(res2$log2FoldChange > 2.5)] <- 'high'

# modify keyvals for transcripts with fold change < -2.5

lour)[which(res2$log2FoldChange < -2.5)] <- 'low'

unique(lour))

p2 <- EnhancedVolcano(res2,

lab = rownames(res2),

x = 'log2FoldChange',

y = 'pvalue',

lectLab = 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 = lour,

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

<(gp=gpar(fill=NA))

本文发布于:2023-06-30 20:48:19，感谢您对本站的认可！

本文链接：https://www.wtabcd.cn/fanwen/fan/82/1070860.html

上一篇：出国高中成绩证明模板

下一篇：Academiclistening02.

标签：结果表达基因分析赋值差异效应

留言与评论（共有 0 条评论）