---
title: 火山图美化
tags: []
id: '2016'
categories:
  - - uncategorized
date: 2022-09-09 13:10:49
---

用[DESeq2](https://occdn.limour.top/2315.html)计算的结果，有时候绘图时，会发现有些差异基因的P值为0，这时候可以通过使用[其他方法](https://occdn.limour.top/2171.html)来计算P值，代替为0的p值，在不改变结果的前提下让[火山图](https://occdn.limour.top/1568.html)更好看。

## 读入数据

```R
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](https://occdn.limour.top/2171.html)、[f\_counts2TMM](https://occdn.limour.top/2159.html)

```R
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值

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

## 绘图保存

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