occdn.limour.top/source/_posts/火山图美化.md

---

title: 火山图美化 tags: [] id: '2317' categories:

• - 绘图 date: 2022-09-09 21:06:05 ---

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

读入数据

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_limma、f_counts2TMM

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值

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']

绘图保存

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)