PowerBacGWAS: a computational pipeline to perform power calculations for bacterial genome-wide association studies.
The main results reproduced: recomputed values matched the published ones within tolerance.
Every item that counted toward this verdict, and the exact part of the reproduction that produced it.
- ✓Same input data as the authors
- ✓No authors-side cause for any deviation
- ✓Reported values are derivable from the shared data
- ✓Any deviation was negligible
- ✓The central claim held under reproduction
- 🟡Reported values were only indirectly comparable
- 🟡A deviation arose in the data or preprocessing
- 🟡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 (all 4 in-scope results graded within-tol). Software/methods paper (PowerBacGWAS) reproduced 1:1 by running the pipeline on the authors' own shipped clade-scale example data via the repo's OWN scripts (github:francesccoll/powerbacgwas@693db85 / zenodo 5950535), recovering the documented power-vs-sample-size surfaces. BOTH pipeline approaches covered, BOTH input types (VCF variants + Roary pangenome), 3 species: (1) SUB-SAMPLING - R1 E. faecium pangenome streptomycin/aadE («job», 1680 LMM runs) and R2 M. tuberculosis burden isoniazid/katG-gene1945 («job», 2497 burden runs); (2) PHENOTYPE-SIMULATION - R4 K. pneumoniae pangenome («job», 30000 runs) and R3 K. pneumoniae individual-variant («job» pyseer + 2214689 aggregate/plot, 30000 runs). Every reproduced power surface matches its reference image structure-for-structure (side-by-side compares in artifacts/{efm_R1,mtb_R2,R4_kpn_pg,R3_kpn_var}_compare.png): same axes, AF/effect-size curve families, legend, and 80%-power-line crossings; causal genes significant in powered combos; observed AF/OR & variants_tested match the wiki. R3's pyseer step hit the 12h association walltime DURING aggregation only (pyseer outputs were complete: 30000 files/27482 nonempty), so aggregation+plot were re-run as a light job on the existing outputs - no recompute of pyseer. Grades are within-tol because phensim/sub-sampling are stochastic (random phenotype/sample draws): per-point power% varies run-to-run but the surface shape and required-N structure reproduce. Env: conda env310 python 3.10.20 + pyseer/plink/gcta/pastml/R4.5, close to paper's pins (paper pinned py>=3.6.9). NOT attempted (and why): upstream raw-read->VCF/pangenome assembly (out of scope - pipeline takes these as shipped inputs); cluster CPU-hour/runtime figures (hardware-dependent, not a reproducible numeric claim); species-wide Table 2/3 full-data values (heavy; kpn_spe VCF lives on FigShare, clade-scale tutorials reproduced instead). All grades PROVISIONAL - a human reviewer signs off.
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 85assessed: 2026-06-22 ⛓ 7dc28cc6e16e
✎ 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-22
- Rubric version
- v1.0
- Assessed by
-
🤖 AI curator · claude (ai-curator headless) · v1.0 · run #1 2026-06-22no 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: sonnetExisting collections of whole-genome sequenced bacterial isolates can be used to perform power calculations for bacterial GWAS, and the paper tests whether two novel approaches (sub-sampling and phenotype-simulation) can determine the sample sizes required to detect genotype-phenotype associations of varying allele frequency, effect size, heritability and homoplasy.
- ★ Two computational approaches (sub-sampling and phenotype-simulation) can be implemented to perform power calculations for bacterial GWAS using existing genome collections, packaged as the PowerBacGWAS pipeline method
- ★ Larger effect size and allele frequency of causal variants require smaller sample sizes to detect via GWAS finding
- ★ Burden testing has more power than SNP-level (variant) GWAS, detecting mutated genes down to 2.5% MAF not detectable by SNP GWAS finding
- ★ Higher phenotype heritability reduces sample sizes needed to detect common variants but has little/no effect on detection of rarer variants finding
- ★ Higher degree of homoplasy of causal SNPs increases the power of GWAS to detect them finding
- ★ Lower bacterial population diversity (single-clade vs species-wide) reduces the sample size required to detect causal variants of the same MAF and effect size finding
- PowerBacGWAS pipeline and user documentation are made publicly available for application to other bacterial populations resource
- Power calculation approaches developed for human GWAS cannot be applied to bacteria due to clonal reproduction, strong population structure and uneven recombination mechanism
| Assay | System | Perturbation | Readout | Platform |
|---|---|---|---|---|
| Pan-genome GWAS (sub-sampling approach) | Enterococcus faecium (species-wide n=1432, single-clade n=761) | sequential sub-sampling of known AMR phenotype labels (kanamycin/streptomycin resistance) to reduce sample size, AF and effect size | power (proportion of replicates above Bonferroni-corrected significance threshold) vs sample size | — |
| Pan-genome GWAS (sub-sampling approach) | Klebsiella pneumoniae (species-wide n=2628) | sub-sampling of meropenem resistance phenotype labels | power vs sample size to detect blaKPC | — |
| Burden test GWAS (sub-sampling approach) | Mycobacterium tuberculosis (species-wide n=2655, single-clade n=1139) | sub-sampling of isoniazid resistance phenotype labels | power vs sample size to detect katG mutations | — |
| Pan-genome GWAS (phenotype-simulation approach) | E. faecium, K. pneumoniae, M. tuberculosis (species-wide and single-clade collections) | simulated binary phenotypes from randomly selected genes at defined MAF, effect size (OR) and heritability | minimum sample size required for 80% power | — |
| Variant (SNP) GWAS (phenotype-simulation approach) | E. faecium, K. pneumoniae, M. tuberculosis (species-wide and single-clade collections) | simulated phenotypes from randomly selected SNPs at defined MAF and effect size | sample size required for 80% power | — |
| Burden GWAS (phenotype-simulation approach) | E. faecium, K. pneumoniae, M. tuberculosis (species-wide and single-clade collections) | simulated phenotypes from mutated genes at defined MAF and effect size | sample size required for 80% power | — |
| Homoplasy-stratified variant GWAS | E. faecium and other bacterial populations | causal SNPs selected with varying degrees of homoplasy (1-5 vs 50-100 independent phylogenetic occurrences) | sample size/power to detect SNPs of given MAF and effect size | — |
| Whole-genome sequence/population diversity characterization | E. faecium, K. pneumoniae, M. tuberculosis genome collections (species-wide vs single-clade) | none | pan-genome size, number of SNP sites, average pairwise genetic distance (SNPs/kb), linkage disequilibrium | — |
- – 500-700 genomes sufficient to detect pan-genome AMR genes of moderate (OR=5) to very large (OR=100) effect size present in ≥10% of population 500-700 genomes
- – Genes of small effect size (OR=1.5) could not be detected using the maximum sample sizes available in the collections
- ▲ Burden testing detected mutated genes down to 2.5% MAF, which was not detectable by SNP-level GWAS
- – Increasing heritability sharply decreased sample sizes needed to detect common genes (25% frequency) but had little/no effect on rarer genes (2.5% frequency)
- ▼ Highly homoplasic SNPs (acquired 50-100 times) in E. faecium at 10% MAF detectable with half the sample size needed for low-homoplasy SNPs (acquired 1-5 times) of same MAF ~2-fold reduction in required sample size
- ▼ Lower sample sizes required to detect causal mutated genes in M. tuberculosis single-clade (lower diversity) population than species-wide population, for genes of same MAF and effect size
- – Known AMR causal variants used for sub-sampling approach showed very large effect sizes and highly significant GWAS associations (e.g. aph(3')-IIIa OR=1083, p=8.25×10^-145; ant(6)-Ia OR=8986, p=1.61×10^-51; blaKPC OR=180, p=8.90×10^-110; katG OR=220, p=2.54×10^-101) OR up to 8986
- – K. pneumoniae single-clade ST288 population could not be used for meropenem resistance power calculation due to unbalanced resistant (95.4%) vs susceptible (1.3%) cases
- fold_change OR = 1083 (aph(3')-IIIa gene and kanamycin resistance, E. faecium species-wide)
- pvalue 8.25 × 10^-145 (GWAS p-value for aph(3')-IIIa association with kanamycin resistance)
- fold_change OR = 8986 (ant(6)-Ia/aad(6) gene and streptomycin resistance, E. faecium single-clade)
- fold_change OR = 180 (blaKPC and meropenem resistance, K. pneumoniae species-wide)
- fold_change OR = 220 (katG mutations and isoniazid resistance, M. tuberculosis species-wide)
- count 500-700 genomes (sample size to detect pan-genome AMR genes OR=5-100 at ≥10% frequency)
- other half the sample size (highly homoplasic (50-100x) vs low homoplasy (1-5x) SNPs at 10% MAF in E. faecium)
- mean 5.6 SNPs/kb (average pairwise genetic diversity, E. faecium species-wide population)
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/tool paper describing PowerBacGWAS, a computational pipeline for performing power calculations for bacterial GWAS using two simulation-based approaches (sub-sampling of real phenotype-genotype data, and phenotype simulation from existing genetic variants). Statistical significance in each simulated GWAS replicate was assessed via association test p-values compared against a Bonferroni-corrected genome-wide significance threshold, and power was reported as the proportion of replicates in which the causal variant exceeded that threshold across varying sample sizes, minor allele frequencies (MAF), effect sizes (odds ratios), heritability, and homoplasy levels. Results are presented primarily as tables and plots of required sample sizes for 80% power, alongside descriptive summaries (AF, OR, p-values) of the real genotype-phenotype relationships used as the basis for simulations.
| Test | Applied to | n | Assumptions |
|---|---|---|---|
| Genome-wide association test (pan-genome GWAS) with Bonferroni-corrected significance threshold | Detection of acquired AMR genes in E. faecium and K. pneumoniae (Fig. 2a, b; Table 2) | n=1432 (E. faecium species-wide), n=761 (single-clade); n=2628 (K. pneumoniae species-wide), n=1193 (single-clade) | not stated |
| Burden test GWAS with Bonferroni-corrected significance threshold | Detection of mutated causal genes in M. tuberculosis and other species (Fig. 2c; Table 3) | n=2655 (species-wide), n=1139 (single-clade) | not stated |
| Variant/SNP GWAS with Bonferroni-corrected significance threshold | Detection of individual causal SNPs across all three species (Supplementary Fig. 2; Table 3) | varies by species/population as listed in Table 1 | not stated |
| Power estimation via replicate proportion exceeding significance threshold | All power calculations (sub-sampling and phenotype-simulation approaches) | varies per simulated sample size scenario | not stated |
-
Multiple testing across variants within each simulated GWAS was controlled using a Bonferroni correction.↳ Could also: A Benjamini-Hochberg false discovery rate (FDR) correction could also be used — FDR-based methods are less conservative than Bonferroni and are often preferred in GWAS with large numbers of correlated variants, potentially increasing power to detect true associations while still controlling for multiplicity.
-
Effect sizes for real genotype-phenotype relationships are reported as point-estimate odds ratios without accompanying confidence intervals.↳ Could also: Reporting odds ratios alongside 95% confidence intervals could also be used — Confidence intervals convey the precision of the effect size estimate in addition to its magnitude, which can be informative for interpreting how variable the OR might be, especially for rarer variants.
-
Power calculations rely on simulation-based replicate counting (proportion of replicates exceeding a significance threshold) rather than closed-form analytic power formulas.↳ Could also: Analytic power calculation formulas (where applicable, e.g., for simple case-control designs) could also be used as a complementary check — Analytic formulas can offer faster, formula-based estimates for simpler scenarios and may serve as a cross-validation of simulation-based results, though bacterial population structure often necessitates simulation-based approaches as the paper notes.
-
Genetic diversity and linkage disequilibrium are summarized using median R2 with interquartile range.↳ Could also: Reporting the full distribution (e.g., via a violin or histogram plot) or additional dispersion measures such as the full range could also be used — This can give readers a fuller picture of the spread and any skewness in pairwise LD values beyond the median and IQR.
-
The pipeline uses a fixed genome-wide Bonferroni threshold across pan-genome, burden, and variant GWAS types.↳ Could also: Permutation-based empirical significance thresholds could also be used — Permutation approaches can account for the specific correlation structure (e.g., due to population structure or linkage) of a given dataset, which may offer an alternative to a fixed Bonferroni threshold in some GWAS contexts.
-
Population structure effects were assessed by comparing species-wide versus single-clade collections rather than through a model-based correction term.↳ Could also: Linear mixed models (LMMs) with a kinship/relatedness matrix could also be used, as mentioned as a methodological advance in the Introduction — LMMs can explicitly model fine-scale population structure and lineage effects within a single analysis, which the paper notes as an alternative approach used elsewhere in bacterial GWAS methodology.
What was reproduced
The exact results taken into scope, with each reported value next to the value our attempt produced.
Scope — pmid-35338232 (PowerBacGWAS)
Paper: Coll F et al. (2022) PowerBacGWAS: a computational pipeline to perform power calculations for bacterial genome-wide association studies. Commun Biol. DOI 10.1038/s42003-022-03194-2 · PMCID PMC8956664.
Nature of the paper: This is a software/methods paper. It presents a pipeline (PowerBacGWAS) and demonstrates it on three bacterial species (E. faecium, K. pneumoniae, M. tuberculosis), at species-wide and single-clade scales. The "results" are power curves / required-sample-size tables produced BY the pipeline. So reproducing the paper = running the pipeline on the shipped data and recovering the same power outputs.
Code: https://github.com/francesccoll/powerbacgwas (own code; Snakemake-style python+R scripts + a Nextflow wrapper + Docker image). Zenodo 5950535 is just a v1.0.0 snapshot (34.2 MB zip = the repo). Wiki has step-by-step tutorials.
Data shipped in repo (data/): multi-sample VCFs (.vcf.gz), Roary
pan-genome tables (.gene_presence_absence.Rtab.gz), Newick trees, phenotype
.phen files, causal-loci lists, and pre-computed power-parameter tables, for:
efm_clade/efm_spe(E. faecium clade A1 / species-wide)kpn_clade/kpn_spe(K. pneumoniae CC258 / species-wide)mtb_clade/mtb_spe(M. tuberculosis lineage4.3 / species-wide) Pre-computed reference output PNGs live inimages/(the wiki tutorial outputs). Note:kpn_spe.vcf.gz(63 MB) is NOT in the repo — it is on FigShare (doi:10.6084/m9.figshare.13353236.v1). The other 5 datasets are complete in-repo.
In scope (pipeline-derived, reproducible)
The wiki tutorials run the full pipeline on the SHIPPED example data and produce specific power outputs. These are the clean, low-hanging reproductions:
| # | Result | Pipeline | Data | Reference output |
|---|---|---|---|---|
| R1 | Sub-sampling power curve, E. faecium pan-genome (streptomycin/aadE) | prepare_gwas_runs_subsampling_roary.py → pyseer → process_gwas_runs → plot | efm_clade pan-genome | images/efm_clade.pg.subsampling.wiki.power.png |
| R2 | Sub-sampling power curve, M. tuberculosis variant GWAS (isoniazid) | subsampling (vcf) → pyseer → process → plot | mtb_clade VCF | images/mtb_clade.var.subsampling.wiki.power.png |
| R3 | Phenotype-simulation power, K. pneumoniae variant GWAS (effect-size) | prepare_gwas_runs.py (GCTA sim) → pyseer → process → plot | kpn_clade VCF | images/kpn_clade.var.phensim.plot_type1.efs.power.png |
| R4 | Phenotype-simulation power, K. pneumoniae pan-genome (effect-size) | roary phensim → pyseer → process → plot | kpn_clade pan-genome | images/kpn_clade.pg.phensim.plot_type1.efs.power.png |
| R5 | Burden-test GWAS runs, K. pneumoniae | burden pipeline → plot | kpn_clade | images/kpn_clade.burden.gwas_runs.results.plot.png |
| R6 | Ancestral-state-reconstruction / homoplasy power (Fig 3 concept) | ancestral_state_reconstruction.py (PastML) + sampling-by-homoplasy | efm_clade tree+vcf | (Fig 3 in paper) |
Primary reproduction targets: R1–R4 (well-documented tutorials, shipped data, self-contained). The agreement criterion is reproducing the power-vs-sample-size relationship and the key required-N for the documented effect (not bit-identical PNGs — GWAS power sim is stochastic; tolerance on power %).
Headline paper numbers (partial / hard)
- Table 2 (samples for 80% power, pan-genome): species-wide values (E. faecium 200, K. pneumoniae 300, M. tuberculosis 500 @10% freq OR=100; and 1100/1200/1000 @2.5%). These use the FULL species-wide datasets (N=1432/2628/2655). efm_spe & mtb_spe full data ARE shipped → attemptable but heavy (many GWAS runs). kpn_spe needs the 63 MB FigShare VCF.
- Table 3 (variant vs burden, 5% MAF OR=5): same, full data, heavy.
- Table 1 (descriptive: AF 56.3%, OR 180, p=2.54e-101): derived from GWAS on full data; treat as PARTIAL/external — cross-check against shipped phenotype + causal-loci files whe
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.
This is a software/methods paper: the in-scope targets were reproduced 1:1 on the authors' own shipped clade-scale data (zenodo 5950535) using the repo's own scripts, recovering all four power-vs-sample-size surfaces structure-for-structure, with causal genes significant in powered combos and observed AF/OR/variants_tested matching the wiki. The deviations are entirely on our/technical side — stochastic simulation noise, version pins, a chosen Bonferroni denominator, and param grids that are supersets of the published figures — not authors' defects. The values are derivable from shared data and the core claim holds fully, so q5/q7 are green. Overall is yellow rather than green because the comparison is qualitative/structural (no exact-numeric 1:1) and species-wide Table 2/3 values were not attempted (FigShare dependency).
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.