A comparative study of techniques for differential expression analysis on RNA-Seq data.
The main results reproduced, with only marginal, non-material deviations.
Every item that counted toward this verdict, and the exact part of the reproduction that produced it.
- ✓Same input data as the authors
- 🔴Reported values were only indirectly comparable
- 🟡A deviation arose in the data or preprocessing
- 🔴A deviation was attributed to the published material
- 🟡Reported values were not (fully) derivable from the shared data
- 🟡The deviation was non-trivial in magnitude
- 🟡The central claim did not (fully) hold under reproduction
- 🟡Overall, the reproduction showed a material discrepancy
A 0–100 reproducibility-quality score from the per-question grades, shown as a z-score: standard deviations above (+) or below (−) the mean of comparable assessments.
▸Reproduction agent’s raw note
Reproduced the core DE pipeline (TopHat -> HTseq-count -> DESeq/edgeR, plus TopHat -> Cuffdiff2) of Zhang et al. 2014 (PMID 25119138) on the paper's assigned no-replicate MAQC dataset (GSE24284/GSE24283, GSM597210 HBR + GSM597211 UHR, 'Wu' technical replicate). Raw read counts exactly match the paper's own supplementary Table S1 (53,238,798 hbr / 59,461,348 uhr), validating correct dataset selection. TopHat alignment (79.9%/82.2% mapping rates), htseq-count (cross-validated against TopHat's multi-mapping counts), DESeq (blind/fit-only/local dispersion, 0 genes padj<0.05, reflecting the method's known conservatism without replicates), and edgeR (fixed BCV=0.4 per the paper's stated value, 4,317/63,677 genes FDR<0.05) all completed successfully. Cuffdiff2 on pooled BAMs found 193/63,652 genes significant (q<0.05). Benchmarked all three methods against a TaqMan qRT-PCR gold standard (596 genes, sourced via the bioconductor-seqc package as a documented substitute for raw GSE5350 files) using AUC1 (full ROC) and AUC2 (FPR<=0.05 restricted ROC, the paper's own definitions, independently confirmed against its supplementary Table S4): AUC1 = 0.95/0.92/0.92 and AUC2 = 0.59/0.61/0.59 for DESeq/edgeR/Cuffdiff2 respectively. This reproduces the paper's qualitative finding that edgeR performs comparatively best at the restricted, decision-relevant FPR range, though DESeq narrowly leads on full-ROC AUC1 in our run. Exact published MAQC-specific AUC1/AUC2 numbers were not extractable from the paper's text, figures-as-images, or its five supplementary tables (S1-S5, all individually checked), so several quantitative claims are graded 'partial' rather than 'exact' pending human reviewer access to Figure 2/Figure 7 source data. No drops; all assigned pipeline components were run to completion on real data.
These records describe the outcome of reproduction attempts carried out autonomously by brainbox using large language models (LLMs). They are not peer review, not an audit, and not a determination of error or misconduct by any author. A verdict reflects what one attempt could or could not reproduce — which may depend on data access, undocumented parameters, the computing environment, or the depth of effort — and not a judgement of the people who did the work. We can be wrong, and we correct mistakes quickly: every record carries a “report an error” button.
✎ I am an author of this paper
Updated or fixed a deposit, or is there an erratum? Ask us to re-run the metrics. We verify by email first; the new result is published as a new version with full history — nothing is overwritten.
Provenance — full disclosure
When this reproduction was carried out, which methodology version was used, and by whom — so the record can be audited and checked independently.
- Reproduced
- 2026-08-02
- Rubric version
- v1.0
- Assessed by
-
🤖 AI curator · claude (ai-curator room) · v1.0 · run #1 2026-08-02no human curator yet
- Last updated
- 2026-08-02
Provisional, curator- or AI-assessed, and independently checkable. A reproduction outcome states what one attempt could reproduce — not a judgement of the authors.
Deep full-text extraction
Model: opusThere is no consensus on which software tool is optimal for detecting differentially expressed genes from RNA-Seq data, nor on how to design an efficient RNA-Seq experiment. This study asks how Cufflinks-Cuffdiff2, DESeq and edgeR perform for differential expression analysis when the number of biological replicates, sequencing depth, and balanced vs. unbalanced sequencing depth within and between groups are varied.
- ★ edgeR performs slightly better than DESeq and Cuffdiff2 in terms of the ability to uncover true positives. finding
- ★ DESeq, or taking the intersection of DEGs from two or more tools, is recommended when the number of false positives is a major concern. finding
- ★ edgeR is slightly preferable for differential expression analysis in other circumstances, at the expense of potentially introducing more false positives. finding
- ★ Cufflinks-Cuffdiff2, which takes BAM files rather than count matrices as input, had been omitted from previous comparison studies; this study includes it alongside the count-matrix-based DESeq and edgeR. method
- ★ Performance is benchmarked against sets of DEGs identified independently by quantitative RT-PCR (MAQC data) or by microarray (mouse neurosphere data). method
- The Poisson distribution underestimates variation when biological replicates are included (overdispersion), so the negative binomial distribution, with a scaling factor for the variance, has achieved a dominant position for modelling RNA-Seq feature counts. mechanism
- ★ A new in-house mouse neurosphere RNA-Seq dataset (K_N, KCl vs. norepinephrine) plus matched microarray data are released as a benchmarking resource (NCBI SRA SRX516577; GEO GSE57440). resource
- Analysis pipelines and custom PERL scripts for subsampling datasets and generating BAM files from simulated count values are publicly available at https://github.com/Qiongyi/RNA-Seq-comparison/. resource
| Assay | System | Perturbation | Readout | Platform |
|---|---|---|---|---|
| Single-end RNA-Seq (MAQC dataset I and II) | Human brain reference RNA (hbr) and universal human reference RNA (uhr) biological samples; no biological replicates | none (hbr vs. uhr sample comparison) | Mapped read counts per gene / differentially expressed genes between hbr and uhr | Illumina GAII (35 bp reads, Dudoit lab; 50 bp reads, Wu lab) |
| Paired-end RNA-Seq (K_N dataset) | Primary neurospheres from hippocampal tissue of adult (8–12 week) male C57BL/6J mice, cultured in complete neurosphere medium with EGF and bFGF; 4 biological replicates per group | Drug/ion treatment: norepinephrine (10 µM) vs. potassium chloride (15 mM), neurospheres collected on day 14 | Mapped read counts per gene / DEGs between K and N groups | Illumina HiSeq2000; Illumina TruSeq RNA Sample Preparation Kit v2 (cat#RS-122-2001); TRIzol Reagent; Ambion DNA-free kit; QuBit RNA assay kit (Q32852); Agilent RNA 6000 Pico Kit on Agilent 2100 Bioanalyser (RIN>8) |
| RNA-Seq (LCL1 dataset, real data used for false positive rate estimation) | Lymphoblastoid cell lines from Yoruban individuals (HapMap); 40 samples with >8 M reads sequenced at Yale, out of 69 cell lines | none (samples randomly assigned to two hypothetical treatment groups of N = 20, 14, 8, 6, 5, 4, 3, 2, 1) | Number of DEGs called / false positive rate, averaged over 10 independent simulations | Illumina GAII |
| Simulated RNA-Seq (LCL2 dataset) | In silico simulation based on the LCL1 lymphoblastoid cell line dataset; BAM files generated by randomly adding/removing aligned reads with custom PERL scripts and SAMtools | DEGs simulated in a random 10% of total genes; counts scaled by exp{(−1)^i δ_j}, δ_j from a two-component normal distribution | True positive rate and false positive rate against known DEG/non-DEG sets, averaged over 10 simulations | — |
| Subsampled RNA-Seq depth/balance experiments (K_N subsets and LCL3 dataset) | Subsets of K_N_full (K_N_30M, 20M, 10M, 5M; unbalanced between groups 1 and 2) and of LCL2 (S1_8M_balanced, S2_5M_balanced, S3_1M_balanced, S4_5M_1M_btw, S5_5M_1M_within) | Varying sequencing depth; unbalanced sequencing depth between and within groups | Number of DEGs, TPR, FPR, AUC1 (0≤FPR≤1) and AUC2 (0≤FPR≤0.05) per tool | — |
| Quantitative RT-PCR (benchmark/gold standard) | MAQC uhr and human brain samples, four technical replicates each (GEO GSE5350, GSM129638–GSM129645) | none | log2 fold change between hbr and uhr; positive set |log2FC|>2, negative set |log2FC|<0.2 | — |
| cDNA microarray (benchmark) | Mouse neurospheres grown with norepinephrine or KCl, n = 3 each (two biological replicates in common with the RNA-Seq data) | Norepinephrine (N) vs. potassium chloride (K) treatment | log2 fold change and P-value; positive set |log2FC|>1 and P<0.05, negative set |log2FC|<0.1 and P>0.1 | Affymetrix GeneChip Mouse Gene 1.0 ST arrays, Affymetrix GeneChip Scanner; Applause WT-Amp ST kit (NuGEN), Encore Biotin module (NuGEN); Partek Genomics Suite |
| Read alignment and count quantification pipeline | Human reference genome hg19 (MAQC, LCL) and mouse genome mm10 (K_N) | none | Aligned BAM files for Cuffdiff2 and raw count tables for DESeq/edgeR; DEG lists; AUC of ROC curves | Tophat v2.0.8; HTseq-count v0.5.4p2; Cufflinks-Cuffdiff2 v2.1.1; DESeq v1.10.1; edgeR v3.0.8 (Bioconductor 2.12, R 2.15.1) |
- ▲ edgeR uncovers more true positives than DESeq and Cuffdiff2
- – edgeR's higher sensitivity comes at the expense of potentially introducing more false positives, so DESeq or the intersection of two or more tools is preferred when false positives are a major concern
- – MAQC dataset I contains 137.77 M single-end 35 bp reads and dataset II contains 112.70 M single-end 50 bp reads, both without biological replicates 137.77 M and 112.70 M reads
- – The in-house K_N dataset yielded 476.38 M paired-end 101 bp x2 reads across 4 biological replicates per group 476.38 M reads (231.55 M K, 244.83 M N)
- – qRT-PCR benchmark yielded 410 genes in the positive set and 86 genes in the negative set out of 1044 genes 410 positive / 86 negative of 1044
- – Microarray benchmark yielded 77 genes in the positive set and 9072 genes in the negative set out of 35557 genes 77 positive / 9072 negative of 35557
- – No differentially expressed genes are expected between the two hypothetical LCL1 groups because samples were randomly drawn from the same population, allowing false positive rate estimation
- – Performance was quantified by AUC1 (area under the ROC over 0≤FPR≤1, maximum 1) reflecting overall DEG identification and AUC2 (0≤FPR≤0.05, maximum 0.05) reflecting performance in the discovery-relevant FPR range AUC1 max = 1; AUC2 max = 0.05
- count 137.77 million RNA-Seq reads (read length 35 bp) (MAQC dataset I, Dudoit group, Illumina GAII)
- count 112.70 million RNA-Seq reads (read length 50 bp) (MAQC dataset II, Wu group, Illumina GAII)
- count 476.38 M paired-end reads (231.55 M for K, 244.83 M for N) (K_N_full mouse neurosphere dataset, 4 biological replicates per group)
- count 410 genes in the positive set and 86 genes in the negative set among 1044 genes (qRT-PCR benchmark; positive |log2FC|>2, negative |log2FC|<0.2)
- count 77 genes in the positive set and 9072 genes in the negative set among 35557 genes (K_N microarray benchmark; positive |log2FC|>1 and P<0.05, negative |log2FC|<0.1 and P>0.1)
- count 40 samples selected (>8 M reads each) from 69 sequenced Yoruban cell lines (LCL1 dataset used for false positive rate estimation)
- other 10% of total genes simulated as DEGs; counts scaled by exp{(−1)^i δ_j}, δ_j ~ two-component normal with μ = (−0.5, 0.5) and σ = (0.7, 0.7) (LCL2 simulation design)
- count 10 independent simulations, averaged for number of DEGs, FPR and TPR (LCL1, LCL2 and LCL3 analyses)
Statistical methods review
Model: sonnetA neutral, descriptive read of the statistical approach — what was done, and (for shared learning, not as criticism) what could also have been done.
This is a methods-comparison paper that benchmarks three RNA-Seq differential expression (DE) analysis tools (Cuffdiff2, DESeq, edgeR) using both real datasets (MAQC, an in-house K_N mouse neurosphere dataset, and a Yoruban lymphoblastoid cell line population dataset) and simulated data with known ground-truth DEGs. Each tool's built-in statistical model (beta negative binomial for Cuffdiff2; negative binomial with local regression for DESeq; negative binomial with empirical Bayes dispersion moderation for edgeR) was used to call DEGs, and performance was evaluated against qRT-PCR- and microarray-derived positive/negative gene sets using ROC curves and area-under-curve (AUC1 full range, AUC2 restricted to FPR≤0.05) rather than a single hypothesis test. Simulation-based results (LCL2/LCL3) were averaged over 10 independent simulation runs, and effects of replicate number, sequencing depth, and balanced vs. unbalanced depth were assessed via true/false positive rates.
| Test | Applied to | n | Assumptions |
|---|---|---|---|
| Negative binomial Wald/likelihood-based DE test (DESeq) | DEG calling for MAQC, K_N, and LCL datasets | varies by dataset (e.g., 0 biological replicates for MAQC technical data; 4 biological replicates for K_N; 1-20 hypothetical replicates for LCL1) | not stated |
| Negative binomial DE test with empirical Bayes dispersion moderation (edgeR) | DEG calling for MAQC, K_N, and LCL datasets | varies by dataset (see above) | not stated |
| Beta negative binomial model-based DE test (Cuffdiff2) | DEG calling for MAQC, K_N, and LCL datasets | varies by dataset (see above) | not stated |
| ROC/AUC (AUC1 full range, AUC2 restricted to FPR≤0.05) | Comparing each tool's DEG calls to qRT-PCR (MAQC) or microarray (K_N) gold-standard positive/negative gene sets | 410 positive / 86 negative genes (qRT-PCR benchmark); 77 positive / 9072 negative genes (microarray benchmark) | not stated |
| Fixed log2 fold-change (and, for microarray, P-value) thresholds to define positive/negative gold-standard gene sets | Construction of benchmark sets from qRT-PCR and microarray data | 1044 genes (qRT-PCR); 35557 genes (microarray) | not stated |
-
Tool performance was evaluated using ROC curves and AUC computed over the full FPR range and a restricted low-FPR range (AUC2), against benchmark sets that are highly imbalanced (e.g., 77 positive vs. 9072 negative genes in the microarray benchmark).↳ Could also: Precision-recall (PR) curves and area under the PR curve — PR curves are often preferred to ROC/AUC when the positive class is much rarer than the negative class, since they can be more sensitive to differences in performance on the minority (positive/DEG) class.
-
Gold-standard positive and negative gene sets were defined using hard log2 fold-change (and, for microarray, P-value) cutoffs rather than using the qRT-PCR/microarray measurements as continuous quantities.↳ Could also: Rank or linear correlation (e.g., Spearman or Pearson correlation) between RNA-Seq-derived and qRT-PCR/microarray-derived fold changes across all genes — A continuous correlation-based comparison uses information from all genes rather than only those passing a threshold, and avoids sensitivity to the specific cutoff values chosen to define the positive/negative sets.
-
All three evaluated tools model RNA-Seq counts with a negative binomial (or beta negative binomial) distribution and test for differential expression directly on counts.↳ Could also: A linear-model-based approach such as limma-voom, which transforms counts and associated precision weights for use in linear models with empirical Bayes moderated t-statistics — Linear-model-based approaches can be convenient for more complex experimental designs (e.g., covariates, multi-factor models) and provide moderated t/F-statistics within a familiar linear-modeling framework.
-
Simulation-based results (LCL2/LCL3 datasets) were summarized as average values (e.g., mean TPR/FPR/AUC) across 10 independent simulations.↳ Could also: Reporting the spread across simulation runs (e.g., SD or a range/CI of AUC values across the 10 runs) alongside the means — Showing variability across simulation replicates would convey how consistent each tool's performance was from run to run, in addition to the average performance.
-
The number of biological replicates in the LCL1 false-positive-rate analysis was varied descriptively (N = 1 to 20) to observe how false positive rate changes with replicate number.↳ Could also: A formal statistical power/sample-size analysis for RNA-Seq (e.g., using tools such as RNASeqPower or PROPER) to estimate the replicate number needed for a target power and false discovery rate — A power-analysis-based approach can translate the same underlying variance/dispersion information into a specific replicate-number recommendation for a desired detection power, complementing the descriptive replicate sweep used here.
-
Genome-wide DEG calls from DESeq, edgeR, and Cuffdiff2 were compared to benchmark sets without an explicitly stated multiple-testing correction step in the analyzed text.↳ Could also: Explicitly reporting a Benjamini-Hochberg false discovery rate (FDR) threshold applied across all tested genes — Because DE analysis tests thousands of genes simultaneously, explicitly stating and reporting an FDR-controlling procedure (commonly built into DESeq/edgeR output) can make the multiple-testing handling transparent to readers comparing tools.
What was reproduced
The exact results taken into scope, with each reported value next to the value our attempt produced.
Assessments & scoring basis
Each contributor’s verdict, the per-question basis, and the auditable, itemised worksheet behind it.
An automated assessment. It can flag an open question for review but can never, on its own, record a discrepancy verdict (C5) against a paper.
Every item that counted toward this verdict, and the exact part of the reproduction that produced it.
Input side is exemplary: the merged 7-lane MAQC read counts (hbr=53,238,798 / uhr=59,461,348) match the paper's own Table S1 'Wu' column exactly, and htseq-count's __alignment_not_unique cross-validates against TopHat's multi-alignment counts, so dataset identity and alignment are not in question. The deviation sits at the endpoint, on the authors' reporting side: the paper publishes its MAQC AUC1/AUC2 comparison only inside Figure 2/Figure 7 images with a prose summary, and none of Tables S1-S5 carry the numbers — so our AUC1 (DESeq 0.9535 / edgeR 0.9229 / Cuffdiff2 0.9233) and AUC2 (0.5877 / 0.6067 / 0.5884) can be computed but never checked. Severity is moderate, not critical: the paper's 'edgeR slightly better' trend is reproduced on AUC2 but reversed on AUC1, within a narrow 0.92-0.95 band, and the claim was stated as pooled across three datasets. Two self-chosen steps (DESeq fitType='local' after non-convergence; TaqMan gold standard via bioconductor-seqc instead of raw GSE5350) are legitimate and documented but add reproducer-side degrees of freedom — no fabrication signal anywhere.
Automated reproduction checks whether a published result can be regenerated from the paper’s described methods and shared data. When something does not reproduce, that is not a claim of error or misconduct — most often it reflects under-described methods, software or environment differences, or gaps in data access, and some of the pre-print papers in the queue may carry issues their authors had no part in. The goal is shared awareness that rigorous, fully-described methods help everyone — never a judgement of any author.
Are you an author? We would genuinely like to hear from you — to clarify the record, add data or code, re-run the pipeline after an accession update, and publish your response right next to the assessment. Everything here is open and auditable.
🚩 Report an error in this record
Spotted something wrong — a verdict you’d contest, a data or value error, or a private detail that slipped through? Tell us, with a short justification. Authors and readers are equally welcome to write in; we review every report.
Prefer email, or the form below not working? Contact us at support@doesitreproduce.com.