Analyzing RNA-seq data with DESeq2：输入数据和差异表达分析

• 基于DESeq2分析RNA-seq数据
• • Abstract
  • 标准流程
  • 快速上手
  • 如何获取DESeq2的帮助
  • 致谢
  • 资金支持
• 输入数据
• • 为何必须输入非标准化（非均一化）的counts值？
  • DESeqDataSet
  • 导入转录本丰度文件以及tximport文件生成的矩阵
  • Tximeta用于自动导入元数据
  • 导入count matrix矩阵
  • 导入htseq-count数据
  • SummarizedExperiment输入
  • 预过滤
  • 关于因子的说明
  • 折叠技术重复
• 差异表达分析
• • Log fold change shrinkage的可视化和排名
  • 加速与并行化思路
  • P 值和调整后的 p 值
  • 独立假设权重
  • 对差异分析结果进行绘图
  • • MA-plot MA图
  • p-value设定为NA的情况：
• 参考资料

基于DESeq2分析RNA-seq数据

Abstract

从 RNA-seq 中分析计数数据的基本任务是检测差异表达的基因。计数数据以表格的形式显示，每个样本报告了分配给每个基因的序列片段的数量。类似的数据也出现在其他分析类型中，包括comparative ChIP-Seq、 HiC、 shRNA screening和质谱法。一个重要的分析问题是，与条件内变异性相比，条件之间系统性变化的量化和推论统计学。软件包 DESeq2提供了用负二项广义线性模型检验差分表达式的方法; dispersion 和logarithmic fold changes的估计包含了数据驱动的先验分布。本文解释了包的使用，并演示了典型的工作流。Bioconductor 网站上的 RNA-seq 工作流程涵盖了类似于本文的内容，但是速度较慢，包括从 FASTQ 文件生成计数矩阵。2包版本: 1.34.0

标准流程

注意: 如果你在已发表的研究中使用DESeq2，请注明:

Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. 10.1186/s13059-014-0550-8

其他具有类似目的的Bioconductor软件包还有edgeR、limma、DSS、EBSeq和baySeq。

快速上手

这里我们展示了差异表达式分析的最基本步骤。在DESeq2的上游有各种步骤，这些步骤会生成每个样本的计数或估计计数，我们将在下面的章节中讨论。这个代码块假设您有一个称为 cts 的 count 矩阵和一个称为 coldata 的样本信息表。design表明如何建模样本，这里，我们想要测量条件的影响，控制批次的差异。两因素变量batch和condition应该是coldata。

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~batch + condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, contrast=c("condition","treated","control"))
# or to shrink log fold changes association with condition:
res <- lfcShrink(dds, coef="condition_trt_vs_untrt", type="apeglm")

• cts是count matrix
• coldata是table of sample information
• design是指明如何对样本进行比较的信息。

上述代码中是为了在控制批次效应的同时研究condition间的效应。

下面将解释以下启动函数:

• 如果您已经执行了转录量化（使用Salmon、kallisto、RSEM等)，您可以使用tximport导入数据，它会生成一个列表，然后您可以使用DESeqDataSetFromTximport()。
• 如果您使用txmeta导入量化数据，它会生成一个含有额外元数据的summarizeexperiment，那么您可以使用DESeqDataSet()。
• 如果您有htseq-count文件，您可以使用DESeqDataSetFromHTSeq()。

如何获取DESeq2的帮助

任何和所有关于DESeq2的问题都应该张贴在 Bioconductor 支持网站上，该网站作为问题和答案的可搜索知识库:Https://support.bioconductor.org
发布一个问题并用“DESeq2”标记将自动向包作者发送一个提醒，以便在支持站点上做出响应。有关如何构建信息性文章的信息，请参阅常见问题(FAQ)列表中的第一个问题。
你不应该把你的问题通过电子邮件发送给软件包的作者，因为我们只是回复这个问题应该发布到 Bioconductor 的支持网站上。

致谢

Constantin Ahlmann-Eltze 为提高DESeq2的计算性能和构建glmGamPoi包的接口贡献了核心代码。在DESeq2的开发过程中，我们受益于许多人的帮助和反馈，包括但不限于：Bionconductor核心小组，亚历杭德罗•雷耶斯(Alejandro Reyes)、安德烈•奥莱斯(Andrzej Oles)、亚历山德拉•佩科沃斯卡(Aleksandra Pekowska)、费利克斯•克莱因(Felix Klein)、尼古拉斯•伊格纳蒂亚迪斯(Nikolaos Ignatiadis)、,朱(阿佩格尔姆) ，约瑟夫 · 易卜拉欣(阿佩格尔姆) ，文斯 · 凯里，欧文 · 索尔伯格，鲁平 · 太阳，德文 · 赖安，史蒂夫 · 利亚诺格鲁，杰西卡 · 拉尔森，克里斯蒂娜 · 查沃拉波尔，潘杜，理查德 · 伯根，威廉 · 塔洛恩，伊琳 · 维德沃尔，汉内克 · 范德沃尔，托德 · 伯威尔，杰西 · 罗利,伊戈尔 · 多尔加莱夫，斯蒂芬 · 特纳，瑞安 · c · 汤普森，蒂尔 · 维斯纳-汉克斯，康拉德 · 鲁德，大卫 · 罗宾逊，明祥 · 腾，马蒂亚斯 · 莱什，索纳利 · 阿罗拉，乔丹 · 拉米洛夫斯基，伊恩 · 德沃金，比约恩 · 格鲁宁，瑞安 · 麦思维，保罗 · 戈登，莱昂纳多 · 科拉多，恩里科 · 费雷罗，彼得 · 兰费尔德，加文 · 凯利，罗布 · 帕特罗，夏洛特 · 索尼森，科恩 · 范登伯格,范尼 · 佩罗多，大卫 · 里索，斯蒂芬 · 恩格尔布雷希特，尼古拉斯 · 阿尔卡拉，杰里米 · 西蒙，特拉维斯 · 皮塔塞克，罗里 · 基什内尔 · 巴特勒，本 · 基思，丹 · 梁，尼尔 · 艾根，罗里 · 诺兰，迈克尔 · 舒伯特，雨果 · 塔瓦雷斯，埃里克 · 戴维斯，万森 · 穆，张成。

资金支持

Deseq2及其开发商得到了欧洲联盟第七个框架方案通过 RADIANT 项目、 NIH NHGRI r01-hg09937和 CZI EOSS 奖提供的部分资助。

---

输入数据

为何必须输入非标准化（非均一化）的counts值？

作为输入，DESeq2包期望得到的计数数据，例如，从RNA-seq或另一个高通量测序实验，以整数矩阵的形式。该矩阵中第i行和第j列的交点代表在样品j中有多少条reads落在基因i上（就是说矩阵中样品为列，基因为行）。类似地,对于其他类型的分析,矩阵的行可能如对应绑定地区(ChIP-Seq)或肽序列(定量质谱)。我们将在下面的章节中列出获取计数矩阵的方法。

矩阵中的值应该是未归一化计数或测序reads(单端RNA-seq)或片段(成对RNA-seq)的估计计数。RNA-seq工作流程描述了制备这种计数矩阵的多种技术。提供计数矩阵作为DESeq2统计模型(Love, Huber, and Anders 2014)的输入是很重要的，因为只有计数值才能正确评估测量精度。DESeq2的统计模型会自动根据矩阵大小进行校正，因此根据矩阵大小转化后的或者标准化（均一化）的矩阵不能作为DESeq2的输入数据。

DESeqDataSet

在统计分析期间，DESeq2包用来存储read counts和中间估计数量的对象类是DESeqDataSet，在这里的代码中，它通常表示为一个对象dds。

一个技术细节是DESeqDataSet类扩展了SummarizedExperiment包的RangedSummarizedExperiment类。“Ranged”部分指的是分析数据的行(这里是计数)可以与基因组范围(基因的外显子)相关联。这种关联有助于对结果的下游探索，利用其他Bioconductor封装的range-based的功能(例如，找到与差异表达基因最接近的ChIP-seq峰)。

一个DESeqDataSet对象必须关联相应的design公式。design公式指明了要对哪些变量进行统计分析。该公式（上文中的design = ~batch + condition）以短波浪字符开头，中间用加号连接变量。design公式可以修改，但是相应的差异表达分析就需要重新做，因为design公式是用来估计统计模型的分散度以及log2 fold change的。

注意：要将你感兴趣的变量放在design公式的后面，对照变量放在前面。

以下介绍4种不同pipiline下构建DESeqDataSet方法：

1. 从转录本丰度文件以及tximport构建
2. 从count matrix 构建
3. 从htseq-count文件构建
4. 从SummarizedExperiment对象构建

导入转录本丰度文件以及tximport文件生成的矩阵

该方法是最新的及推荐的方法。 具体是将上游的快速转录本丰度定量软件获得的基因水平的count matrix，然后用tximport包导入，这种方法支持的软件有：

• Salmon
• Sailfish
• kailisto
• RSEM

用以上软件进行转录本丰度估计的优点有：

1. 校正了潜在的不同样本间基因长度导致的错误。
2. 运行速度快，与alignment-based的方法相比需要的内存和硬盘空间少
3. 避免将那些比对到带有同源序列的多个基因上的片段丢弃，从而增加敏感度

以下我们演示如何将Salmon软件导出的定量文件“quant.sf”导入到DESeq2中。

注意，与使用 system.file 定位目录不同，用户通常只需提供一个路径，例如/path/to/quant/files。对于典型的使用，条件信息应该已经显示为示例表示例的一列，而在这里，我们构建人工条件标签以供演示。

library("tximport")
library("readr")
library("tximportData")
dir <- system.file("extdata", package="tximportData")
samples <- read.table(file.path(dir,"samples.txt"), header=TRUE)
samples$condition <- factor(rep(c("A","B"),each=3))
rownames(samples) <- samples$run
samples[,c("pop","center","run","condition")]

接下来，我们使用适当的样本列指定文件的路径，并读取将转录本与该数据集的基因连接起来的表格。

files <- file.path(dir,"salmon", samples$run, "quant.sf.gz")
names(files) <- samples$run
tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v27.csv"))

利用 tximport 函数导入 deseq2所需的量化数据。关于 tximport 使用的更多细节，包括建立 tx2基因表将数据集中的转录本与基因相连接，请参考 tximport 包装 vignette。

txi <- tximport(files, type="salmon", tx2gene=tx2gene)

最后，我们可以从txi对象和样本中的样本信息构造一个DESeqDataSet。

library("DESeq2")
ddsTxi <- DESeqDataSetFormTximport(txi, colData = samples, design =~condition)

这里的ddsTxi对象可以在下面的分析步骤中用作dds。

Tximeta用于自动导入元数据

另一个 Bioconductor 软件包，tximeta (Love et al. 2020) ，扩展了 tximport，提供了相同的功能，还提供了为常用转录组自动添加注释元数据(GENCODE，ensemble bl，RefSeq for human and mouse)的额外好处。有关更多细节，请参见 tximatpackage 插件。使用 DESeqDataSet 函数可以很容易地将一个 summarydexperiment 加载到 deseq2中，在 txmeta 包 vignette 中有一个示例，如下:

coldata <- samples
coldata$files <- files
coldata$names <- coldata$run

library("tximeta")
se <- tximeta(coldata)
ddsTxi <- DESeqDataSet(se, design = ~ condition)

这里的 ddsTxi 对象可以在以下分析步骤中用作 dds。如果 tximate 将引用转录组识别为具有预先计算的散列校验和的转录组之一，那么 dds 对象的 rowrangs 将被预先填充。同样，请参阅 tximetavignette 以获得完整的细节。

导入count matrix矩阵

用DESeqDataSetFromMatrix命令导入count matrix。除count matrix外，使用者还需要将count matrix的列名单独建立一个DataFrame，以及一个design formula。 以下我们演示DESeqDataSetFromMatrix的用法，主要过程是从pasilla包中加载count data，读取的count matrix命名为cts，样品信息表命名为coldata。下边演示的方法是从featureCounts的输出结果中提取以上信息。

library("pasilla")
pasCts <- system.file("extdata", "pasilla_gene_counts.tsv",package="pasilla",mustWork=TRUE)
pasAnno <- system.file("extdata", "pasilla_sample_annotation.csv", packages="pasilla", mustWork=TRUE)
cts <- as.matrix(read.csv(pasCts, sep="\t", row.names="gene_id"))
coldata <- read.csv(pasAnno,row.names=1)
coldata <- coldata[,c("condition", "type")]

coldata$condition <- factor(coldata$condition)
coldata$type <- factor(coldata$type)

我们检查计数矩阵和列数据，看看它们在样本顺序方面是否一致。

head(cts,2)

coldata

注意 count matirx 的列与coldata的行的顺序一致，DESeq2默认这两者间的顺序是一致的。
上面的结果中可以看出cts和coldata中列和行的顺序并不一致，另外coldata中的行名也多了“fb”两个字母，需要调整顺序并且去掉多余的字母。

rownames(coldata) <- sub("fb", '''', rownames(coldata)) #去掉fb
all(rownames(coldata)%in%colnames(cts)) 
[1] TRUE  #此行为运行结果
all(rowname(coldata)==colname(cts))
[1]FALSE #此行为运行结果
cts <- cts[, rownames(coldata)] #直接把列名提换了吗？该列的数据提换了吗？
all(rowname(coldata)==colname(cts))

如果在 Rsubread 包中使用了 featureCounts 函数(Liao、 Smyth 和 Shi 2013) ，那么可以直接从列表输出中的“ counts”元素提供读取计数矩阵。通常可以使用基本的 r 函数(如 read.csv 或 read.delim)将 count 矩阵和列数据从平面文件读入 r。对于 htseq-count 文件，请参阅下面的专用输入函数。

有了count matrix（cts）和sample information（coldata），我们就可以构建DESeqDataSet矩阵了

library("DESeq2")
dds <- DESeqDataSetFromMatrix(countData=cts, colData = coldata, design = ~condition)
dds

如果您有其他特性数据，可以通过添加到新构建的对象的元数据列，将其添加到 DESeqDataSet。(这里我们添加冗余数据只是为了演示，因为基因名称已经是 dds 的行名。)

featureData <- data.frame(gene=rownames(cts))
mcols(dds) <- DataFrame(mcols(dds), featureData)
mcols(dds)

导入htseq-count数据

通过HTSeq软件进行read count的结果数据，可以用DESeqDataSetFromHTSeqCount导入到DESeq2中。首先要将htseq-count软件的结果文件的位置信息保存在一个变量中。下文中我们用pasilla包中的htseq-count结果文件作为例子。

directory <- "/path/to/your/files/"

directory <- systme.file("extdata", package = "pasilla", mustWork =TRUE)
sampleFiles <- grep("treated",list.files(directory), value=TRUE) #用list.file指明which file to read，再用grep选择包含"treated"的文件
sampleCondition <- sub("(.*treated).*". "\\1", sampleFiles)      #用sub函数从sample filename得到condition status
sampleTable <- data.frame(sampleName =sampleFiles, fileName =sampleFiles, condition =sampleCondition)
#build the DESeqDataSet
library("DESeq2")
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory= directory, design=~condition)
ddsHTseq #显示该矩阵

SummarizedExperiment输入

如果已经创建或获得了SummarizedExperiment，那么可以很容易地将它输入到DESeq2中，如下所示。首先，我们加载包含气道数据集的包。

library("airway")
data("airway")
se <- airway

下面的构造函数显示了从RangedSummarizedExperiment se生成一个DESeqDataSet。

library("DESeq2")
ddsSE <- DESeqDataSet(se, design = ~ cell + dex)
ddsSE

预过滤

虽然提前去掉low count基因不是运行DESeq进行差异表达分析必须的，但是提前过滤掉low count基因有两个好处：减少dds矩阵的大小，提高DESeq运行的速度。 这里我们演示如何去掉read count在10以下的基因。 注意：更严格的过滤会在DESeq的results功能中自动实现

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

关于因子的说明

默认情况下，R会安装字母顺序选择一个level作为参考level。就是说，如果没有给DEseq2函数指定要与哪个level进行比较（即指定哪个level是对照组），默认会安装字母顺序选择排在最前面的level作为对照组。

如果要指定哪个level是对照组，可以直接在函数#results#的参数contrast中指明对照level（下文中会具体介绍），或者直接在factor level中指明。

注意，更改design公式中变量的factor level只能在运行DESeq2分析前，分析中或分析后都不能再进行更改。更改dds矩阵中的factor level可以用以下两个方法：

#用factor
dds$condition <- factor(dds$condition, levels = c("untreated", "treated"))
#或者用relevel，直接指定参考level
dds$condition <- relevel(dds$condition, ref = "untreated")

如果需要从dds中提取一个子集，比如说从dds中去掉部分样品，这时有可能会同时去掉所有样品的一个或者多个level。这种情况下，可以用droplevles函数去掉没有对应sample的level：

dds$condition <- droplevels(dds$condition)

折叠技术重复

DESeq2包中的collapseReplicates函数可以用于将技术重复的定量数据折叠成一个样品的。注意不能把这个函数用到生物学重复上，更多信息请查看collapseReplicates的帮助页面。

差异表达分析

标准的差异表达分析已经被打包整合到DESeq这一个函数中。其具体分析过程可以通过“?DESeq”在DESeq包的帮助页面查看。

差异表达分析的结果可以用results函数进行提取，该函数导出一个结果表格，**结果包括log2 fold change值、p值（Wald检验）以及校正后的p值。**如果不指定额外的参数，结果表格会输出disign公式中最后一个变量的log2 fold change以及wald检验的p值，如果该变量是一个因子，则将该变量的最后一个水平与第一个水平做比较。但是如果用results函数中的name和contrast参数指定了哪一组和哪一组做比较（参见results的帮助页），design公式中变量的顺序就不重要了。 我们看一下结果表：

dds <- DESeq(dds)  #进行差异表达分析
res <- results(dds)     #提取分析结果并赋予变量res
res  #显示差异分析结果表 

结果如下：

如上图所示，在表格上方有log2 fold change以及p-value值是以哪一组比哪一组得到的说明。
注意，我们可以使用下面的等价命令来指定我们想要建立的结果表的系数或对比度:

res <- results(dds, name="condition_treated_vs_untreated")
res <- results(dds, contrast=c("condition","treated","untreated"))

这两个命令等价的一个例外是，使用对比将在两组比较中额外将估计的LFC设置为0，其中两组中的所有计数都等于0(而其他组的计数为正)。由于这可能是将这些情况下的LFC设置为0的理想特性，因此可以使用对比来构建这些结果表。关于从拟合的DESeqDataSet对象提取特定系数的更多信息可以在帮助页面的结果中找到。下面还将进一步讨论对比论证的使用。

Log fold change shrinkage的可视化和排名

effect size 的缩小(LFC 估计值)对于基因的可视化和排序是有用的。为了缩小 LFC，我们将 dds 对象传递给函数 lfcShrink。下面我们指定使用 apeglm 方法计算effect size 收缩(Zhu，Ibrahim，and Love 2018) ，这是对以前估计方法的改进。

我们提供 dds 对象和要shrink的系数的名称或数字，其中的数字指的是在 resultsNames (dds)中显示的系数的顺序。

resultsNames(dds)

resLFC <- lfcShrink(dds, coef="condition_treated_vs_untreated", type="apeglm")
resLFC

加速与并行化思路

对于大多数分析来说，上述步骤的执行时间应该不超过30秒。对于具有复杂设计和许多样本(例如几十个系数，大约100个样本)的实验，人们可能希望获得比DESeq默认运行提供的更快的计算速度。我们有两个建议:

1. 通过使用参数 fitType = “ glmGamPoi”，可以利用 Constantin Ahlmann-Eltze 编写的更快的 NB GLM 引擎。请注意，在 deseq2中，glmGamPoi 的接口需要使用 test = “ LRT”和简化设计的规范。
2. 可以利用并行计算的优势。通过加载 BiocParallel 包，然后设置以下参数，可以很容易地实现 DESeq、 results 和 lfcShrink 的并行化: parallel = TRUE 和 BPPARAM = MulticoreParam (4) ，例如，在4个核上分割作业。然而，一些关于并行化的建议: 首先，建议过滤所有样本计数较低的基因，以避免向子进程发送不必要的数据，当这些基因功耗较低并且无论如何都会被独立过滤; 其次，由于向子进程发送数据的开销，往往需要增加更多的报酬递减，因此我建议首先从少量额外的核开始。注意，获得 resultsNames (dds)中列出的系数或对比结果的速度很快，不需要并行化。作为 BPPARAM 的替代方案，可以在分析开始时注册核心，然后在调用时指定 parallel = TRUE。

library("BiocParallel")
register(MulticoreParam(4))

P 值和调整后的 p 值

我们可以根据最小的p值对结果表进行排序:

resOrdered <- res[order(res$pvalue),]

我们可以使用summary函数总结一些基本的计数。

summary(res)

查看padj小于0.1的基因有多少个

sum(res$padj < 0.1, na.rm=TRUE)

results函数有很多可以调整的参数，可以通过？results了解每个参数的意义。

需要注意的是，results函数自动基于每个基因的标准化的count的平均数进行indenpendent filtering（下文会详细介绍）。对那些padj值低于给定FDR阈值（参数alpha）的基因的数量进行优化。参数alpha的默认值是0.1。如果padj值的阈值要取0.1以外的值，参数alpha也要调整成对应的值。

res05 <- results(dds, alpha=0.05)
summary(res05)

sum(res05$padj < 0.05, na.rm=TRUE)

独立假设权重

p值滤波思想的一个推广是用权值假设来优化功率。Bioconductor软件包IHW可实现独立假设加权法(Ignatiadis等人，2016)。在这里，我们展示了使用IHW来调整DESeq2结果的p值。有关更多详情，请参阅IHW包的简介。IHW结果对象存储在元数据中。

注:如果发表的研究中使用独立假设加权的结果，请引用:

Ignatiadis, N., Klaus, B., Zaugg, J.B., Huber, W. (2016) Data-driven hypothesis weighting increases detection power in genome-scale multiple testing. Nature Methods, 13:7. 10.1038/nmeth.3885

# (unevaluated code chunk)
library("IHW")
resIHW <- results(dds, filterFun=ihw)
summary(resIHW)
sum(resIHW$padj < 0.1, na.rm=TRUE)
metadata(resIHW)$ihwResult

对于高级用户，请注意，deseq2包计算的所有值都存储在 DESeqDataSet 对象或 DESeqResults 对象中，下面将讨论对这些值的访问。

对差异分析结果进行绘图

MA-plot MA图

在DESeq2中，plotMA函数显示了DESeqDataSet中所有样本归一化计数的平均值上由给定变量引起的log2倍变化。如果调整后的p值小于0.1，点将被涂成红色。落在窗口外的点被绘制为指向上或向下的开放三角形。

plotMA(res, ylim=c(-2,2))

MA图常用在芯片数据的分析中。M是纵轴，代表log2 fold change。A是横轴，代表标准化后的count值的平均数。红色的点代表padj值小于0.1，落在窗口外的点用小三角表示。

对log2 fold change进行收缩一下，得到的MA图会好看一些。

resLFC <- lfcShrink(dds, coef=2)
plotMA(resLFC, ylim=c(-2,2))

lfcShrik是DESeq中对log2 fold change进行收缩的函数，它有三个参数：dds矩阵，coef协同系数，contrast与results函数中的contrast的用法一致。可以参考以下命令：

res.shrink <- lfcShrink(dds, contrast = c("condition","KD","control"), res=res)

上面的MA图收缩后就会变成这个样子：
绘制MA图后，还可以通过identify函数获得图中感兴趣的点所代表的基因。点击图中感兴趣的点，然后运行以下代码：

idx <- identify(res$baseMean, res$log2FoldChange)
rownames(res)[idx]

关于收缩的算法，lfcShrik提供了三个参数：normal，apeglm和ashr

• normal 是DESeq2包原始的收缩估计量（shrikage estimator），自适应正态先验分布（adaptive normal prior）
• apeglm是apeglm包中的收缩估计量，自适应t先验分布（adaptive t prior）
• ashr是ashr包中的收缩估计量。

我们比较一下这三种收缩算法：

resApe <-lfcShrink(dds, coef=2,type="apeglm")
resAsh <-lfcShrink(dds, coef=2,type="ashr")
par(mfrow=c(1,3), mar=c(4,4,2,1))
xlim <- c(1,1e5); ylim<-c(-3,3)
plotMA(resLFC, xlim=xlim, ylim=ylim, main="normal")
plotMA(resApe, xlim=xlim, ylim=ylim, main="apeglm")
plotMA(resAsh, xlim=xlim, ylim=ylim, main="ashr")

p-value设定为NA的情况：

• 在同一行中，所有样品的count值都为0，这时baseMean为0，log2 fold change、p-value以及adjusted p-value都可设为NA
• 如果一行中包含了一个落在可信区域之外的极值，p-value 以及 adjusted p-value可设定为NA。极值的确定可根据Cook距离，也可自行设定，具体方法下文会介绍。
• 如果一行因为含有过低的mean normalized count而被自动独立过滤（automatic independent filtering）过滤掉了，这时可把adjusted p-value设定为NA

参考资料

1. https://zhuanlan.zhihu.com/p/31952724
2. https://zhuanlan.zhihu.com/p/32668252
3. https://blog.csdn.net/xiaomotong123/article/details/106900481
4. https://www.jianshu.com/p/9ac7f6aa6027
5. https://djnb.gitee.io/2020/12/14/P019_Bioconductor_DESeq2_Package/