Integrative transcriptomics and single-cell transcriptomics analyses reveal potential biomarkers and mechanisms of action in papillary thyroid carcinoma.
The main result did not reproduce in this reproduction attempt. Where our recomputation produced values that differ from the published ones, those discrepancies are listed below. This is a single automated attempt — not peer review and not a finding of error or misconduct — and differences can also arise from data access, undocumented parameters or the computing environment. The verdict can be contested via “report an error”.
- ✓Same input data as the authors
- ✓Reported values were directly 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
Described well enough to reproduce the core pipeline on the paper's own training data. GSE3467 (9 tumor / 9 normal, matches exactly) re-run through limma -> clusterProfiler -> pROC on «infra» «our HPC». DIRECTIONALLY 1:1: same dataset, same methods, same biology (top GO BP term reproduced is 'hormone metabolic process', the PTC-relevant term the paper highlights; the 3 named biomarkers each reach AUC=1.00, satisfying the reported >0.7). NUMERICALLY DIVERGENT on counts: 679 DEGs vs 413, 653 GO terms vs 384 (CC near-exact 63 vs 60). Cause: the paper's clusterProfiler 4.10.1 pin was unsolvable on bioconda (forces an old bioc that conflicts with org.Hs.eg.db), so the run used clusterProfiler 4.18.4 / R 4.5.3 / org.Hs.eg.db 3.22.0; plus the paper does not specify its probe->gene collapse or pre-filter step. No fabrication signal -- gap is consistent with version drift + underspecified preprocessing. NOT attempted (the hard 20%): WGCNA (6 modules/turquoise R=-0.91/898 hub), DEG-hub intersection (316), PPI+MCODE (12 core genes), ML biomarker selection (LASSO/SVM-RFE/Boruta -> the 3 genes), scRNA-seq GSE191288 (Seurat/CellChat/monocle), RT-qPCR -- all either GUI/web-tool-bound, highly parameter-sensitive with unspecified settings, or wet-lab. C3 KEGG could not run: rest.kegg.jp is unreachable from the compute nodes.
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.
Assessment versions
Every reproduction run is kept as an immutable version — anchored to the data as it stood, with a tamper-evident chain hash. A rerun (e.g. after an author updates a deposit) adds a new version; the previous one stays on record.
-
v1 current initial assessment Score 49assessed: 2026-06-14 ⛓ 4fb8afab2c7a
✎ 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-06-14
- Rubric version
- v1.0
- Assessed by
-
🤖 AI curator · claude (ai-curator room) · v1.0 · run #1 2026-06-15no human curator yet
- Last updated
- 2026-08-05
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: sonnetThe study aims to identify robust transcriptomic biomarkers for papillary thyroid carcinoma (PTC) diagnosis through integrative bioinformatics analyses and to elucidate the cellular mechanisms underlying PTC pathogenesis at single-cell resolution.
- ★ ENTPD1, SERPINA1, and TACSTD2 are potential transcriptomic biomarkers for PTC finding
- ★ These three biomarkers may collectively regulate PTC occurrence and development via cytokine-cytokine receptor interaction pathways mechanism
- ★ Tissue stem cells, epithelial cells, and smooth muscle cells are key cell types in PTC finding
- ★ Epithelial cells interact with tissue stem cells and smooth muscle cells mainly through the COL4A1-CD4 and COL4A2-CD4 ligand-receptor pairs finding
- ★ The collagen signaling pathway is the most dominant intercellular communication pathway among key cells finding
- ★ The three key cell types undergo three distinct differentiation stages with stage-specific expression trends of ENTPD1, SERPINA1, and TACSTD2 finding
- A multi-algorithm consensus of LASSO, SVM-RFE, and Boruta reduces single-method bias and enhances reliability of candidate biomarker screening method
- WGCNA combined with DEG and PPI/MCODE analysis identifies core candidate genes for PTC method
| Assay | System | Perturbation | Readout | Platform |
|---|---|---|---|---|
| bulk RNA-seq differential expression (limma) | PTC tissue vs normal thyroid tissue (GSE3467, training set) | none (tumor vs normal comparison) | differentially expressed genes (DEGs) | — |
| GO and KEGG enrichment analysis (clusterProfiler) | DEGs from GSE3467 | none | enriched biological processes/pathways | clusterProfiler v4.10.1 |
| weighted gene co-expression network analysis (WGCNA) | PTC tissue (GSE3467) | none | gene modules correlated with tumor/normal phenotype, hub genes | — |
| machine learning feature selection (LASSO, SVM-RFE, Boruta) | PTC tissue (GSE3467) | none | candidate biomarker genes | glmnet; caret; Boruta (R packages) |
| Wilcoxon test expression validation and ROC/AUC analysis | PTC/ATC/FTC tissue (GSE3678, GSE33630, GSE65144, GSE82208) | none (tumor vs normal/FTA comparison) | expression consistency and diagnostic AUC of candidate biomarkers | pROC |
| single-cell RNA-seq (Seurat, Harmony, SingleR clustering/annotation) | PTC tissue (GSE191288, 6 tumor + 1 normal) | none | cell clusters and cell-type annotation | Seurat; Harmony; SingleR |
| cell-cell communication analysis (CellChat) | key cells from scRNA-seq (GSE191288) | none | ligand-receptor interactions and signaling pathway strength | CellChat |
| RT-qPCR | PTC cell lines TPC-1 and IHH4 | none | expression levels of ENTPD1, SERPINA1, TACSTD2 (2^-ΔΔCt vs GAPDH) | SYBR Green Mix (Vazyme) |
- – Three candidate biomarkers (ENTPD1, SERPINA1, TACSTD2) identified via DEG/WGCNA/machine-learning consensus and validated by expression and ROC analysis
- – GSEA indicated biomarkers are involved in cytokine-cytokine receptor interaction pathways
- – Tissue stem cells, epithelial cells, and smooth muscle cells identified as key cell types via abundance and biomarker expression differences (p<0.01)
- – Epithelial cells communicate with tissue stem cells and smooth muscle cells via COL4A1-CD4 and COL4A2-CD4 ligand-receptor pairs
- ▲ COL4A1 and COL4A2 highly expressed in epithelial cells; CD4 elevated in tissue stem cells and smooth muscle cells
- – Collagen signaling pathway identified as the most dominant pathway among key cells
- – Pseudotime analysis showed three distinct differentiation stages with stage-specific trends in ENTPD1, SERPINA1, and TACSTD2 expression
- – 203 RNA-seq samples across 5 GEO datasets and 7 scRNA-seq samples from GSE191288 were analyzed
- pvalue adjusted p-value < 0.05, |log2FC| > 1 (DEG screening threshold in GSE3467)
- other R2 = 0.85 (soft-threshold selection) (WGCNA scale-free topology fit)
- correlation |r| > 0.4, P < 0.05 (WGCNA key module-trait correlation criteria)
- other MM > 0.8, GS > 0.6 (WGCNA hub gene selection criteria)
- other AUC > 0.7 (threshold for defining biomarkers in training and validation ROC analyses)
- count 9 tumor/9 normal (GSE3467); 7/7 (GSE3678); 49/45 (GSE33630); 12/13 (GSE65144); 27/25 (GSE82208) (sample sizes of RNA-seq datasets)
- count 6 tumor and 1 normal scRNA-seq datasets (GSE191288 single-cell dataset composition)
- pvalue p < 0.01 (threshold for defining key cells based on biomarker expression differences)
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.
The study applies a multi-stage bioinformatics pipeline to bulk RNA-seq data: limma DEG analysis on a training set (GSE3467, n=18), followed by WGCNA for module discovery, three machine learning feature-selection algorithms (LASSO, SVM-RFE, Boruta) whose intersection defined candidate biomarkers, Wilcoxon-based expression validation across four additional GEO datasets, ROC/AUC assessment, nomogram construction, and GSEA for pathway annotation. A scRNA-seq dataset (GSE191288, n=7 samples) was analysed with Seurat/SingleR for cell-type annotation, CellChat for cell–cell communication, and Monocle2 for pseudotime trajectories. Experimental validation used RT-qPCR in two PTC cell lines with two-tailed Student's t-tests. Results are primarily reported as AUC values, adjusted p-values, and log2 fold changes, with significance thresholds of p<0.05 or p<0.01.
| Test | Applied to | n | Assumptions |
|---|---|---|---|
| limma moderated t-test (DEG analysis) | Tumor vs. normal in training set GSE3467; also used to split high/low biomarker expression groups for GSEA input | 18 (9 tumor, 9 normal) | not stated |
| Pearson correlation (WGCNA module–trait correlation) | Correlation of co-expression modules with tumor/normal phenotype in GSE3467 | 18 (9 tumor, 9 normal) | not stated |
| LASSO logistic regression with 10-fold cross-validation (lambda.min) | Feature selection among candidate genes in GSE3467 | 18 (9 tumor, 9 normal) | not stated |
| SVM-RFE with 10-fold cross-validation repeated 10 times | Feature selection among candidate genes in GSE3467 | 18 (9 tumor, 9 normal) | not stated |
| Boruta (random forest feature importance, 500 iterations) | Feature selection among candidate genes in GSE3467 | 18 (9 tumor, 9 normal) | not stated |
| Wilcoxon rank-sum test | Expression validation of candidate biomarkers in GSE3467, GSE3678, GSE33630, GSE65144, GSE82208; cell-type abundance and biomarker expression comparisons in scRNA-seq GSE191288 | GSE3467 n=18; GSE3678 n=14; GSE33630 n=94; GSE65144 n=25; GSE82208 n=52; GSE191288 n=7 samples | not stated |
| ROC curve / AUC (R package pROC) | Predictive performance of individual biomarkers and their three-gene combination across all five bulk RNA-seq datasets | GSE3467 n=18; GSE3678 n=14; GSE33630 n=94; GSE65144 n=25; GSE82208 n=52 | na |
| Nomogram with calibration curve (R packages rms, Rregplot) | Prediction of PTC probability from three biomarkers in GSE3467 | 18 (9 tumor, 9 normal) | not stated |
| GSEA (R package clusterProfiler, KEGG gene set c2.cp.kegg.v7.4) | Pathway enrichment for high- vs. low-expression groups of each biomarker in GSE3467 | 18 (9 tumor, 9 normal) | not stated |
| Two-tailed Student's t-test | RT-qPCR expression differences between experimental groups in TPC-1 and IHH4 PTC cell lines | not stated | not stated |
-
DEG analysis was performed with limma; the paper does not specify whether the voom transformation for count data was applied↳ Could also: DESeq2 or edgeR could also be used for count-based RNA-seq differential expression; limma-voom is another well-validated limma extension for RNA-seq counts — DESeq2 and edgeR use negative binomial models explicitly designed for discrete read counts, which may better model count-level variance; specifying the exact limma workflow (array-fitted vs. voom) would also clarify which variance model was assumed
-
Multiple Wilcoxon tests were applied separately across candidate genes and across multiple validation datasets without a stated family-wise correction↳ Could also: A Benjamini-Hochberg FDR correction applied across the set of per-gene comparisons could also be used — When testing several genes simultaneously, an explicit multiplicity correction controls the expected false-discovery rate across the family of tests; this is a common complement to individual p-value thresholds in multi-gene validation studies
-
Three machine learning algorithms were run independently and their exact intersection was taken as the final biomarker set↳ Could also: A rank-aggregation approach (e.g., averaging standardized importance scores across methods) or a single elastic-net model could also be used — Requiring exact intersection across all three methods is stringent and may exclude genes that are consistently but not unanimously selected; rank aggregation retains a continuous measure of cross-algorithm agreement and may recover such genes
-
Cell type annotation in scRNA-seq was performed automatically with SingleR using Spearman rank correlation to a reference dataset↳ Could also: Manual annotation using curated cell-type marker genes (e.g., via Seurat FindAllMarkers cross-referenced with published marker lists) could also be used, or both approaches could be combined for cross-validation — Manual marker-gene annotation incorporates tissue- and disease-specific knowledge and is often used to validate or complement automated methods; the quality of SingleR annotation depends on the relevance of the reference dataset to the tissue type
-
Pseudotime trajectory analysis was performed with Monocle2 (reversed graph embedding)↳ Could also: RNA velocity (scVelo) or diffusion pseudotime (DPT) could also be used to infer differentiation trajectories from the same data — RNA velocity estimates the directionality of transcriptional change from spliced/unspliced ratios without requiring a predefined root cell; DPT uses diffusion maps; each method makes different topological assumptions and can yield complementary views of the differentiation process
-
RT-qPCR results were compared with Student's t-test reported at a p<0.05 threshold without explicit dispersion measures or exact n per group↳ Could also: Reporting the mean with SD (or 95% CI) and the exact n alongside the p-value could also be done; if normality cannot be assumed for small cell-line replicates, a Wilcoxon signed-rank test is also used in such contexts — Dispersion measures and exact sample sizes allow readers to assess effect magnitude and variability; for small-n in vitro experiments, many journals now request both the effect estimate and its uncertainty in addition to the p-value
Citation network
Where this publication sits in the reproducibility-weighted citation graph — what it is built on, and what is built on it. Citation data from OpenAlex.
No assessed neighbours yet — the network grows as more papers are assessed.
Data lineage
The datasets this paper uses (text-mined from the full text via Europe PMC), and which other assessed papers stand on the same data. A shared dataset is a factual link — not a judgement.
- Developing a thyroid cancer differentiation st...⚑ L1 No data access ⚑
- Developing a thyroid cancer differentiation st...⚑ L1 No data access ⚑
Downstream reach in the literature
154 downstream papers · 2 datasetsHow widely the datasets deposited by this paper are reused across the whole literature (Europe PMC), beyond our assessed set. This is a factual dependency map — reusing a public dataset is normal, good science. It is not a judgement on the downstream papers; the only verdict here is this paper's own, with its cited rationale.
- Characterizing dedifferentiation of thyroid cancer b... 2021 · 136 cites
- Large Scale Gene Expression Meta-Analysis Reveals Ti... 2016 · 103 cites
- Immune Cell Confrontation in the Papillary Thyroid C... 2020 · 91 cites
- METTL3-mediated m6A modification of STEAP2 mRNA inhi... 2022 · 71 cites
- Senescent thyrocytes and thyroid tumor cells induce... 2019 · 67 cites
- miR30a inhibits LOX expression and anaplastic thyroi... 2015 · 65 cites
- Aberrant lipid metabolism in anaplastic thyroid carc... 2015 · 137 cites
- Characterizing dedifferentiation of thyroid cancer b... 2021 · 136 cites
- Large Scale Gene Expression Meta-Analysis Reveals Ti... 2016 · 103 cites
- Cell Cycle M-Phase Genes Are Highly Upregulated in A... 2017 · 54 cites
- Cancer Associated Fibroblasts and Senescent Thyroid... 2020 · 49 cites
- Cancer-Associated Fibroblasts Positively Correlate w... 2021 · 47 cites
What was reproduced
The exact results taken into scope, with each reported value next to the value our attempt produced.
Scope — pmid-40520228
Paper: Cao W, Gao K, Zhao Y. Integrative transcriptomics and single-cell transcriptomics analyses reveal potential biomarkers and mechanisms of action in papillary thyroid carcinoma. Front Genet 2025. PMID 40520228 / PMC12162626 / DOI 10.3389/fgene.2025.1536198.
Code link in record: github.com/YuLab-SMU/clusterProfiler — this is the third-party enrichment tool the authors used (P16: applying an existing tool to the paper's own data is equally valid). The authors ship no own repo; the analysis is a standard limma → clusterProfiler → WGCNA → PPI/MCODE → ML → ROC → Seurat pipeline described in Methods. We reproduce the clearly-specified pipeline steps on the paper's own GEO data.
Primary data: GEO GSE3467 (Affymetrix HG-U133 Plus 2.0 / GPL570), 9 papillary-thyroid-carcinoma tumor + 9 normal thyroid. This is the training set from which DEGs, enrichment, WGCNA hub genes, biomarkers and ROC all derive.
In scope (attempt — clearly specified, low-hanging 80%)
| id | result | reported value | paper loc | pipeline | priority |
|---|---|---|---|---|---|
| C1 | DEGs in GSE3467 (limma, |log2FC|>1, adj.P<0.05) | 413 total = 228 up / 185 down | Results/Fig 2 | limma | HIGH |
| C2 | GO enrichment of DEGs (clusterProfiler 4.10.1, p<0.05) | 384 terms = BP 280 / CC 60 / MF 44 | Results/Fig 3 | clusterProfiler enrichGO | HIGH |
| C3 | KEGG enrichment of DEGs (clusterProfiler, p<0.05) | 12 pathways | Results/Fig 3 | clusterProfiler enrichKEGG | HIGH |
| C4 | ROC/AUC of the 3 named biomarkers in GSE3467 | all AUC > 0.7 | Results/Fig 7 | pROC | MED |
The 3 biomarkers (ENTPD1, SERPINA1, TACSTD2) are named in the paper, so C4 can be checked directly on GSE3467 without re-deriving them through the full WGCNA→ML chain.
Out of scope / not attempted (the hard ~20%) — and why
- WGCNA (6 modules, turquoise R=-0.91, 2403 genes, 898 hub @ MM>0.8 & GS>0.6): highly parameter-sensitive (soft-threshold power, mergeCutHeight, minModuleSize none of which are fully specified). DEG∩hub=316 depends on it. → not attempted.
- PPI (STRING >0.4) → 308 nodes/447 edges → MCODE → 12 core genes: uses the STRING web DB + Cytoscape/MCODE GUI; version-dependent, non-scriptable as described. → not attempted.
- ML feature selection (LASSO log λ.min=-7.7174, SVM-RFE, Boruta) → exactly the 3 biomarkers: depends on the 12 core genes from the GUI step above; the exact 3-gene landing is the paper's headline and the hardest 20%. → not attempted (instead we verify the named biomarkers' ROC directly = C4).
- scRNA-seq GSE191288 (Seurat QC, 10→6 cell types, CellChat COL4A1/2–CD4, monocle pseudotime): large, many unspecified params (resolution given as 0.15 but QC/integration details partial). → not attempted in this pass.
- RT-qPCR (TPC-1/IHH4), nomogram, GSEA per-biomarker: wet-lab / downstream visual; out of pipeline-reproduction scope.
Reproducibility caveats (auditability)
- DEG count depends on probe→gene collapsing (max/mean) and whether the pre-normalized series matrix vs CEL re-processing is used. Paper does not fully specify; we use the GEO series matrix + standard limma. Expect C1 to land near 413, graded within-tol/partial accordingly — the human reviewer decides.
- C2/C3 counts depend on the org.Hs.eg.db / KEGG annotation snapshot; we pin clusterProfiler 4.10.1 to match the paper and record the annotation version.
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.
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.
Reproduction footprint
claude-opus-4-8Measured resources invested to assess this paper — sanitised (machine class only, no job ids/paths). Compute = HPC accounting (SLURM); tokens = the AI agent's session.