火山图美化

DESeq2计算的结果,有时候绘图时,会发现有些差异基因的P值为0,这时候可以通过使用其他方法来计算P值,代替为0的p值,在不改变结果的前提下让火山图更好看。

读入数据

1
2
3
4
5
6
7
8
9
10
11
12
DEG <- readRDS('../../../DEG/CELL/C42vsLNCaP_EtOH.rds')
DEG <- subset(DEG, !grepl('pseudogene', type) & baseMean > quantile(DEG$baseMean)['25%'])
DEG <- f_dedup_IQR(df = DEG, rowNn = DEG$symbol, select_func = maxbaseMean)
FC_up <- with(DEG, log2FoldChange[log2FoldChange>0])
FC_up <- mean(FC_up) + 0*sd(FC_up)
FC_up
FC_down <- with(DEG, log2FoldChange[log2FoldChange<0])
FC_down <- mean(FC_down) - 0*sd(FC_down)
FC_down

tmp <- load('../../../DEG/CELL/GSE.rdata')
tmp # 'cts_b''geneInfo'

重新计算P值

f_DE_limmaf_counts2TMM

1
2
3
4
TMM <- f_counts2TMM(cts_b)
keep <- rowSums(cts_b) > ncol(cts_b)
r2 <- f_DE_limma(TMM[keep,], geneInfo[keep,], Ct2, Tt2, F)
rownames(r2) <- r2$ID

调整P值

1
2
3
4
5
6
7
DEG[DEG$padj ==0, 'padj'] = min(DEG$padj[DEG$padj > 0])
ID <- DEG[DEG$padj < 1e-200, 'ID']
tmp <- r2[ID, 'P.Value'] / max(r2[ID, 'P.Value'])
tmp <- tmp^6
tmp <- tmp / min(tmp)
# DEG[DEG$padj < 1e-200, 'padj'] <- tmp * ((DEG[DEG$padj < 1e-200, 'padj']*1e200)^0.5)*1e-200
DEG[DEG$padj < 1e-200, 'padj'] <- tmp * DEG[DEG$padj < 1e-200, 'padj']

绘图保存

1
2
3
4
5
6
7
8
9
10
11
require(EnhancedVolcano)
options(repr.plot.width=6, repr.plot.height=4)
p <- EnhancedVolcano(DEG,
lab = rownames(DEG),
selectLab = NA,
x = 'log2FoldChange',
y = 'padj',
pCutoff = 0.01,
FCcutoff = min(abs(FC_down),FC_up)) + theme_bw() + theme(panel.grid=element_blank())
p
ggsave(plot = p, filename = 'C42vsLNCaP_EtOH.pdf', width = 6, height = 4)

火山图美化
https://b.limour.top/2016.html
Author
Limour
Posted on
September 9, 2022
Licensed under