(DESeq2) Why are some p values set to NA?
引入
在上一期奇怪的转录组差异表达矩阵之实验分组中，我们谈到DESeq2输出NA的问题，这周我们仍使用上周 GSE126548-分组差异并不大，这个数据集来进行分析
本文主要参考bioconductor文档：Analyzing RNA-seq data with DESeq2 地址：http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
这里先贴出一般使用DESeq2分析的代码：
counts <- read.table("../11.一些奇怪的转录组差异表达矩阵/1.GSE126548-分组差异并不大//GSE126548_HTSeq_count.txt",sep = '\t',header = TRUE,row.names = 1)
range(counts)
# 0 71407700
# 1根据表达量过滤
keep <- rowSums(counts>0) >= floor(0.75*ncol(counts))
filter_count <- counts[keep,] #获得filter_count矩阵
dim(counts);dim(filter_count)
# [1] 60160 6
# [1] 15065 6
# 2根据基因过滤
# 一般行名不能一样，我们现在行名形式就是gene ENSG 应该没问题
table(duplicated(rownames(filter_count)))
# FALSE
# 15065
# 分组
group_list <- factor(c(rep("withBM",3),rep("noBM",3)),levels = c("withBM","noBM"))
table(group_list)
# group_list
# withBM noBM
# 3 3
# 所以这里其实相当于只对低表达量过滤
# R package DESeq, which utilizes a generalized linear model (GLM) and is applied directly to raw read counts.
# 默认counts所以要求整数 内部有log2等 不需手动
#1 查看分组信息和表达矩阵数据
exprSet <- filter_count
dim(exprSet)
# 加载包
library(DESeq2)
#2 第一步，构建DESeq2的DESeq对象dds(建立DEseq数据矩阵)
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list)
colData
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~ group_list)
dds
#3 第二步，进行差异表达分析
dds2 <- DESeq(dds)
#4 提取差异分析结果，trt组对untrt组的差异分析结果
tmp <- results(dds2,contrast=c("group_list","withBM","noBM"))
DEG_DESeq2_raw <- as.data.frame(tmp[order(tmp$padj),])
head(DEG_DESeq2_raw) # 联系之前的内容DEseq2和edge会存在p值为0的情况
# baseMean log2FoldChange lfcSE stat pvalue padj
# ENSG00000257671 2937.01756 -5.662520 1.0831750 -5.227706 1.716262e-07 0.001160022
# ENSG00000121552 261.23701 6.711915 1.3209047 5.081301 3.748577e-07 0.001266831
# ENSG00000228113 449.33229 -6.395500 1.2846913 -4.978239 6.416533e-07 0.001445645
# ENSG00000233271 70.50139 -5.201629 1.0669793 -4.875098 1.087542e-06 0.001837674
# ENSG00000126005 359.50238 -5.362662 1.1296268 -4.747286 2.061642e-06 0.002438831
# ENSG00000269834 55.09172 -4.616794 0.9745453 -4.737382 2.164964e-06 0.002438831
# 查看结果
# 去除差异分析结果中包含NA值的行
DEG_DESeq2 = na.omit(DEG_DESeq2_raw)
# 为什么会出现NA？
# https://blog.****.net/linkequa/article/details/83116789
# 筛选上下调差异基因，设定阈值
fc_cutoff <-2.0
fdr <- 0.05#p值
DEG_DESeq2$regulated <- "normal"
loc_up <- intersect(which(DEG_DESeq2$log2FoldChange>fc_cutoff),
which(DEG_DESeq2$pvalue<fdr))
loc_down <- intersect(which(DEG_DESeq2$log2FoldChange< -fc_cutoff),
which(DEG_DESeq2$pvalue<fdr))
DEG_DESeq2$regulated[loc_up] <- "up"
DEG_DESeq2$regulated[loc_down] <- "down"
table(DEG_DESeq2$regulated)
# down normal up
# 175 6389 195
可以发现去掉包含NA的行后，影响很大
在这个过程中我们使用了DESeq2的默认设置，也就意味着除了手动根据低表达量和基因名进行过滤，我们还使用了DESeq2内部的自动过滤
- Pre-filtering
在上面的代码中，我们使用 keep <- rowSums(counts>0) >= floor(0.75*ncol(counts))进行过滤，也就是保留那些样本基因表达量不为0不少于样本数四分之三的基因
在官方的文档中对于Pre-filtering，有这么一段描述：
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pre-filtering 虽然在运行DESeq2函数之前没有必要对低计数基因进行预过滤，但有两个原因使预过滤变得有用：1减少了dds数据对象的内存大小，并提高了DESeq2中转换和检测函数的速度；2改善可视化效果。
这是因为DESeq2本身有Independent filtering（更严格）
同时，官方也给出了“a minimal pre-filtering”rowSums(counts(dds)) >=10 , 也就是保留所有样本基因表达量计数和大于10的基因
- Independent filtering
Filtering criteria：
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#indfilttheory
Independent filtering of results：
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#indfilt
DESeq2包的 results函数默认情况下使用归一化计数的平均值作为过滤统计信息来执行独立过滤，找到过滤统计量的阈值，该阈值优化了低于显著性水平α的调整后的p值的数量，未通过过滤阈值的基因的调整后的p值被设置为NA
默认的独立过滤是使用genefilter包的 filtered_p函数执行的，filtered_p的所有参数都可以传递给 results函数。过滤阈值和过滤统计量的每个分位数处的拒绝次数可用作结果返回的对象的元数据metadata
例如，我们可以通过绘制results对象的 filterNumRej属性来可视化优化。results函数在过滤统计量的分位数（归一化计数的平均值）上最大化拒绝次数（调整后的p值小于显著性水平）。所选择的阈值（垂直线）是过滤的最低分位数，对于该分位数，拒绝次数在拟合过滤分位数上拒绝次数的曲线峰值的1个残差标准偏差内：
metadata(tmp)$alpha
# [1] 0.1
metadata(tmp)$filterThreshold
# 54.28571%
# 14.34643
plot(metadata(tmp)$filterNumRej,
type="b", ylab="number of rejections",
xlab="quantiles of filter")
lines(metadata(tmp)$lo.fit, col="red")
abline(v=metadata(tmp)$filterTheta)
可以通过设置 independentFiltering 为 FALSE关闭独立过滤
# Independent filtering can be turned off by setting independentFiltering to FALSE.
tmpNoFilt <- results(dds2, independentFiltering=FALSE)
addmargins(table(filtering=(tmp$padj < .1),
noFiltering=(tmpNoFilt$padj < .1)))
# noFiltering
# filtering FALSE TRUE Sum
# FALSE 6719 0 6719
# TRUE 17 23 40
# Sum 6736 23 6759
查看关闭过滤后的结果：
DEG_DESeq2_NoFilt <- as.data.frame(tmpNoFilt[order(tmpNoFilt$padj),])
# > dim(DEG_DESeq2_NoFilt)
# [1] 15065 6
# > dim(na.omit(DEG_DESeq2_NoFilt))
# [1] 14934 6
# 为什么还是有NA？
# 发现仍为NA的 pvalue就已为NA了
可以发现仍存在一些基因p值为NA，所有这些基因和之前大部分的区别在于，它们的pvalue就已经为NA
- How can I get unfiltered DESeq2 results?
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#how-can-i-get-unfiltered-deseq2-results
Users can obtain unfiltered GLM results, i.e. without outlier removal or independent filtering with the following call:
dds2.2 <- DESeq(dds, minReplicatesForReplace=Inf)
tmpUnFilt <- results(dds2.2, cooksCutoff=FALSE, independentFiltering=FALSE)
DEG_DESeq2_UnFilt <- as.data.frame(tmpUnFilt[order(tmpUnFilt$padj),])
dim(DEG_DESeq2_UnFilt)
# [1] 15065 6
dim(na.omit(DEG_DESeq2_UnFilt))
# [1] 15065 6
这时完全属于unfiltered
可以发现DESeq2内部过滤不仅仅有independent filtering还存在outlier removal
Why are some p values set to NA?
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pvaluesNA
- 如果在一行中，所有样本的计数都为零，则基础平均值（baseMean）列将为零，log2 FC、p值和调整后的p值都将被设置为NA
- 如果一行平均归一化计数较低，会被自动独立过滤掉，只有调整后的p值将被设置为NA
上述两条都很好理解，我们往期推文无论是使用DESeq2、edgeR还是limma，都或多或少考虑到了这些
我们将重点看看outlier removal
- 如果一行包含一个具有极端计数异常值的样本，则p值和调整后的p值将被设置为NA。这些异常计数值由Cook距离检测到。自定义离群值过滤和替换离群值计数并进行重新拟合的功能描述如下
Approach to count outliers http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#outlier
RNA测序数据有时会包含与实验或研究设计明显无关的非常大的计数的孤立实例，这些计数可能被视为异常值。异常值可能产生的原因有很多，包括罕见的技术或实验人工制品，在遗传上不同的样品中读取映射问题，以及真实但罕见的生物事件。
在很多情况下，用户主要关注表现一致的基因，这就是为什么默认情况下，DESeq2会过滤受这些异常值影响的基因，而如果有足够的样本，异常值计数将被替换以进行模型拟合，这两种方式将在下面进行介绍：
DESeq函数对每个基因和每个样本进行计算，用一种叫做Cook距离的异常值诊断检测。Cook距离是衡量单个样本对基因拟合系数影响程度的指标，而较大的Cook距离则表示存在异常值。Cook距离的矩阵存储在 assays(dds)[["cooks"]]中。
results函数会自动标记那些在具有3个或更多重复样本的情况下，包含高于Cooks距离截止值的基因。这些基因的p值和调整后的p值将被设置为NA。至少需要3个重复样本才能进行标记，因为仅有2个重复样本时很难判断哪个样本是异常值。可以使用 results(dds, cooksCutoff=FALSE)命令关闭这个过滤功能。
当*度很大——即样本数远大于要估计的参数数时，完全因为一个计数异常值而从分析中移除整个基因是不可取的。当给定样本的重复次数为7次或更多次时，DESeq函数将自动用所有样本的修剪均值来替换大的Cook距离值，该平均值经过该样本的尺寸因子或正则化因子进行缩放。这种方法是谨慎的，它不会导致假阳性，因为它用空假设预测的值替换异常值。这种异常值替换仅在有7个或更多个重复时发生，并且可以使用 DESeq(dds，minReplicatesForReplace = Inf)命令关闭。
上述行文提到的两种方式的默认Cooks距离截止值取决于样本大小和要估计的参数数量。默认值是使用F(p，m-p)分布的99%分位数（其中p是参数数量，包括截距，m是样本数）。
基因标记的默认值可以使用 results函数的 cooksCutoff参数进行修改。对于异常值替换，在 DESeq中保留原始计数，并将替换计数保存为矩阵，命名为 assays(dds)中的 replaceCounts。请注意，如果在设计中存在连续自变量，则不会自动执行异常值检测和替换，因为我们当前的方法涉及对组内方差进行鲁棒估计，难以简单地扩展到连续协变量。不过，用户可以通过 assays(dds)[["cooks"]]检查Cooks距离，以便在必要时进行手动可视化和过滤。
基因标记 "gene flagging"是指DESeq2在RNA测序数据分析中，针对每个基因对所有样本进行异常值检测将存在异常值的样本标记出来。当一个样本的Cooks距离超过F(p,m-p)分布的0.99分位数时，DESeq2会将其标记为异常值。
Note on many outliers: 如果summary（res）（我们上述代码中应是tmp变量）报告了非常多的异常值（例如数百个或数千个），则可以考虑进一步探索，以查看是否应由于质量不佳而去除一个或几个样本。自动异常值过滤/替换在仅有少量异常值的情况下最为有用。当报告的异常值数量有数千个时，可能更有意义地关闭异常值过滤/替换（使用 DESeq函数中的 minReplicatesForReplace = Inf和 results函数中的 cooksCutoff = FALSE），并进行手动检查：
首先可以制作PCA图来检测单个样本的异常值;其次，可以制作Cook’s距离的箱型图，以查看某个样本是否始终比其他样本高
######Approach to count outliers######
par(mar=c(8,5,2,2))
boxplot(log10(assays(dds2)[["cooks"]]), range=0, las=2)
当然，这里不存在这种情况
小结
在本文中，我们介绍了三种DESeq2结果输出NA的情况：
- 如果在一行中，所有样本的计数都为零，则基础平均值（baseMean）列将为零，log2 FC、p值和调整后的p值都将被设置为NA
- 如果一行平均归一化计数较低，会被自动独立过滤掉，只有调整后的p值将被设置为NA
- 如果一行包含一个具有极端计数异常值的样本，则p值和调整后的p值将被设置为NA。这些异常计数值由Cook距离检测到。自定义离群值过滤和替换离群值计数并进行重新拟合的功能描述如下
大家可以联系自己的表达矩阵和差异分析结果对感兴趣的基因进行解读
同时，我们着重介绍了基因计数异常值的处理，包括小样本（但大于3）中的直接过滤和大样本（大于7，大*度）中类似处理缺失值的拟合替换
关于官方文档，很多内容都非常详尽，感兴趣的同学可以挑自己需要的部分进行阅读
http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
下一篇： 假设检验和 p 值数据