Data used to produce S10 Fig.

We assess inferential quality in the field of differential expression profiling by high-throughput sequencing (HT-seq) based on analysis of datasets submitted from 2008 to 2020 to the NCBI GEO data repository. We take advantage of the parallel differential expression testing over thousands of genes, whereby each experiment leads to a large set of p-values, the distribution of which can indicate the validity of assumptions behind the test. From a well-behaved p-value set π0, the fraction of genes that are not differentially expressed can be estimated. We found that only 25% of experiments resulted in theoretically expected p-value histogram shapes, although there is a marked improvement over time. Uniform p-value histogram shapes, indicative of <100 actual effects, were extremely few. Furthermore, although many HT-seq workflows assume that most genes are not differentially expressed, 37% of experiments have π0-s of less than 0.5, as if most genes changed their expression level. Most HT-seq experiments have very small sample sizes and are expected to be underpowered. Nevertheless, the estimated π0-s do not have the expected association with N, suggesting widespread problems of experiments with controlling false discovery rate (FDR). Both the fractions of different p-value histogram types and the π0 values are strongly associated with the differential expression analysis program used by the original authors. While we could double the proportion of theoretically expected p-value distributions by removing low-count features from the analysis, this treatment did not remove the association with the analysis program. Taken together, our results indicate widespread bias in the differential expression profiling field and the unreliability of statistical methods used to analyze HT-seq data.

本研究针对高通量测序（high-throughput sequencing, HT-seq）差异表达谱分析领域的推断质量展开评估，分析对象为2008年至2020年间提交至美国国家生物技术信息中心（National Center for Biotechnology Information, NCBI）基因表达综合数据库（Gene Expression Omnibus, GEO）的数据集。研究依托针对数千个基因的并行差异表达检验展开分析：每项实验均可生成大量p值，其分布特征可反映检验背后所依赖假设的合理性。研究发现，仅25%的实验呈现出理论预期的p值直方图形态，不过该比例随时间推移已有显著提升。服从均匀分布的p值直方图形态极少出现，此类形态对应实际效应少于100个的情况。此外，尽管多数HT-seq分析流程默认大多数基因未发生差异表达，但37%的实验的π₀值小于0.5，仿佛大多数基因的表达水平发生了改变。多数HT-seq实验的样本量极小，功效普遍不足。然而，估算得到的π₀值与样本量N并未呈现预期的关联，这表明实验在控制错误发现率（False Discovery Rate, FDR）方面普遍存在问题。不同p值直方图类型的占比以及π₀值均与原始作者所使用的差异表达分析程序显著相关。尽管通过移除分析中的低计数特征可使理论预期的p值分布占比翻倍，但该处理并未消除与分析程序的关联。综合来看，本研究结果表明差异表达谱分析领域普遍存在偏倚，且用于分析HT-seq数据的统计方法可靠性不足。