Reproducing Wang et al. (2024) from Dryad dataset tdz08kq5d¶
Paper. Wang J, Wang X, Yang K, Lu C, Fields B, Xu Y, Shen Q, Wei Z, Friman V-P (2024). Phage selection drives resistance-virulence trade-offs in Ralstonia solanacearum plant-pathogenic bacterium irrespective of the growth temperature. Evolution Letters 8(2), 253-266. doi:10.1093/evlett/qrad056. The claims file was written from the version of record (Europe PMC full text, PMC10959482) with its supplement (Figures S1 to S5, Tables S1 to S6).
Dataset. doi:10.5061/dryad.tdz08kq5d, Dryad version id
263142 (dryad:tdz08kq5d/263142). The paper cites the dataset's DOI and no version. Dryad's API lists one
version, 263142 (version number 5, published 2023-11-03), so it is the most recent published version and the
one pinned here. It holds two files: README.md (SHA-256
79a23ead0ade6ec77d36e7d1dc34ed51508417d827175829093609d8740c051a) and
temperature_phage_data_allfigures2.xlsx (SHA-256
6bc2c64abc9e1c65bcdc4080fcaa81c669c39f56a6f63708483e47c563c5c0dd), both read here. Dryad needs a token to
download a file, so the registry's workflow stages the files before the run and the helper reads them from
there, checking each SHA-256 against Dryad's file list (docs/archives.md, "Dryad").
What the dataset holds. One sheet per figure: bacterial OD600 of the 40 selection lines (Figure 1B), phage densities (1C), colony counts of 16 lines at day 12 (1D), the per-clone phage resistance, maximum growth rate and carrying capacity of 128 clones (2A to 2C), twitching, biofilm and disease incidence of 28 clones and the ancestor (2D to 2F), the variant calls of the 28 sequenced clones aggregated by group (Table 1, Table S4), the DESeq2 table and gene lists of Figure 4, the KEGG enrichment table of Figure S5, the adsorption assay of Figure S3, and the curves of Figures S2 and S4. The growth curves the growth parameters were fitted to, the per-clone variant calls, the raw reads and the pilB/pilM data of Table S6 are not in it.
Outside data. The paper names NCBI BioProjects PRJNA1028400 and PRJNA1032726. Their public metadata are read in the data check to say what they hold; no data from them enter a result.
Answered questions. The first run of this attempt asked four questions (questions/dryad-tdz08kq5d.yaml),
and a reader answered them in the claims file: the OD600 time points of Figure 1B are days 2 to 12; the 0.24%
of Figure 1D is read both as a mean of replicate frequencies and from pooled counts; the unit of the Figure 2A
to 2C ANOVAs is the mean of the clones of each colony type within each replicate population, with each clone
as a variant; and the sequencing and RNA-seq claims are data: dataset, tested on the deposited tables. The
claims' sections below follow the criteria as they now stand.
No random numbers are drawn anywhere in this notebook.
| Claim | Data | Status |
|---|---|---|
| fig1b-temperature-density | dataset | tested |
| fig1b-phage-reduces-density | dataset | tested |
| fig1b-phage-effect-diminishes | dataset | tested |
| auc-phage-reduction-by-temperature | dataset | tested |
| fig1c-phage-density-temperature | dataset | tested |
| fig1d-mucoid-lost-with-phage | dataset | tested |
| fig1d-phage-raises-scv-frequency | dataset | tested |
| fig2a-phage-exposed-more-resistant | dataset | tested |
| fig2a-temperature-resistance | dataset | tested |
| fig2a-no-temperature-phage-interaction | dataset | tested |
| fig2a-scv-more-resistant | dataset | tested |
| fig2a-mucoid-25p-replicate-resistant | dataset | tested |
| fig2b-scv-higher-growth-rate | dataset | tested on the deposited fitted values |
| fig2c-scv-lower-carrying-capacity | dataset | tested on the deposited fitted values |
| fig2d-scv-twitching-no-exposure-effect | dataset | tested |
| fig2d-mucoid-25p-reduced-twitching | dataset | tested |
| fig2e-mucoid-biofilm-increased | dataset | tested |
| fig2e-biofilm-no-exposure-effect | dataset | tested |
| fig2f-phage-exposure-reduces-virulence | dataset | tested |
| fig2f-scv-avirulent | dataset | tested |
| fig2f-mucoid-25p-avirulent | dataset | tested |
| sequencing-mutation-count | dataset | tested on the dataset's variant table |
| table1-phcs-in-scvs | dataset | tested on the dataset's variant table |
| table1-pild-in-25p-mucoid | dataset | tested on the dataset's variant table |
| table1-t2ss-in-35p-scv | dataset | blocked: the variant table does not identify clones |
| tables6-pil-mutants-resistant | dataset | blocked: the data are not in the dataset |
| figs3-adsorption | dataset | tested |
| rnaseq-differential-expression-count | dataset | tested on the dataset's DESeq2 table |
| rnaseq-kegg-enrichment | dataset | blocked: the criterion's readings differ on the dataset's KEGG table |
| rnaseq-mutated-genes-expression | dataset | tested on the dataset's DESeq2 table |
| rnaseq-defense-systems | dataset | tested on the dataset's DESeq2 table |
| rnaseq-virulence-genes-down | dataset | tested on the dataset's DESeq2 table |
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
from statsmodels.stats.multicomp import pairwise_tukeyhsd
import phage_temp as pt
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 300)
pd.set_option("display.max_colwidth", 70)
warnings.filterwarnings("ignore", category=FutureWarning)
started = time.time()
wb = pt.load_workbook()
readme = pt.load_readme()
print(f"Dryad {pt.DOI}, version {pt.VERSION}: workbook {pt.WORKBOOK} read, SHA-256 checked; "
f"{len(wb)} sheets, in {time.time() - started:.1f} s")
Dryad 10.5061/dryad.tdz08kq5d, version 263142: workbook temperature_phage_data_allfigures2.xlsx read, SHA-256 checked; 16 sheets, in 0.6 s
Data check¶
Each check of the description against the values is stated under its output.
listing = requests.get(pt.FILES_API, timeout=60)
listing.raise_for_status()
files = listing.json()["_embedded"]["stash:files"]
print("Dryad file list of version", pt.VERSION)
for f in files:
print(f" {f['path']}: {f['size']} bytes, sha-256 {f['digest']}")
assert {f["path"]: f["digest"] for f in files} == {"README.md": pt.README_SHA256, pt.WORKBOOK: pt.WORKBOOK_SHA256}
print("\nEvery sheet of the workbook, as read with its first row as header:")
for name, sheet in wb.items():
print(f"\n[{name}] {sheet.shape[0]} rows x {sheet.shape[1]} columns")
print(" columns:", [str(c) for c in sheet.columns])
Dryad file list of version 263142 README.md: 3770 bytes, sha-256 79a23ead0ade6ec77d36e7d1dc34ed51508417d827175829093609d8740c051a temperature_phage_data_allfigures2.xlsx: 518639 bytes, sha-256 6bc2c64abc9e1c65bcdc4080fcaa81c669c39f56a6f63708483e47c563c5c0dd Every sheet of the workbook, as read with its first row as header: [transfer_OD_FIG 1B] 280 rows x 5 columns columns: ['replicate', 'temperature', 'phage', 'time(days)', 'OD'] [phage_density_FIG 1C] 60 rows x 4 columns columns: ['ID', 'temperature', 'time(days)', 'phagedensity (log10 PFU)'] [scv frequency_FIG 1D] 16 rows x 8 columns columns: ['Temperature', 'phage', 'replicate', 'SCV(CFU/mL)', 'mucoid(CFU/mL)', 'Total(CFU/mL)', 'proportion of SCV', 'proportion of m'] [resistance to phage_FIG 2A] 128 rows x 8 columns columns: ['treatment', 'strain', 'morphotype', 'resistance to phage', 'Unnamed: 4', 'Unnamed: 5', 'strian', 'resistance to phage.1'] [maximum growth rate_FIG 2B] 128 rows x 8 columns columns: ['treatment', 'strain', 'morphotype', 'maxgrowthrate(OD600)', 'Unnamed: 4', 'Unnamed: 5', 'strian', 'maxgrowthrate(OD600).1'] [carrying capacity FIG2C] 128 rows x 7 columns columns: ['treatment', 'strain', 'morphotype', 'maxdensity(OD600)', 'Unnamed: 4', 'strian', 'maxgrowthrate(OD600)'] [twitching_FIG2D] 28 rows x 7 columns columns: ['treatment', 'strain', 'time(h)', 'twitching(mm)', 'Unnamed: 4', 'strain.1', 'twitching(mm).1'] [biofilm_FIG 2E] 28 rows x 8 columns columns: ['treatment', 'morphotype', 'strain', 'biofilm(OD570)', 'Unnamed: 4', 'Unnamed: 5', 'strain.1', 'biofilm(OD570).1'] [disease incidence_FIG 2F] 28 rows x 7 columns columns: ['treatment', 'strain', 'morphology', 'Disease incidence', 'Unnamed: 4', 'strain.1', 'Disease incidence.1'] [snp_indel] 138 rows x 19 columns columns: ['CHROM', 'POS', 'TYPE', 'REF', 'ALT', 'EVIDENCE', 'FTYPE', 'STRAND', 'NT_POS', 'AA_POS', 'EFFECT', 'LOCUS_TAG', 'GENE', 'PRODUCT', 'temperature', 'phage', 'colonial_morphology', 'treatment', 'num_samples_with_snp'] [gene_blast in uniprot] 23 rows x 6 columns columns: ['chrome', 'LOCUS_TAG', 'sequence', 'uniprotKB', 'gene name', 'protein'] [transcriptome] 5147 rows x 18 columns columns: ['Fig 4A', 'Unnamed: 1', 'Unnamed: 2', 'Unnamed: 3', 'Unnamed: 4', 'Unnamed: 5', 'Unnamed: 6', 'Unnamed: 7', 'fig 4B', 'Unnamed: 9', 'Unnamed: 10', 'Unnamed: 11', 'Unnamed: 12', 'fig 4C', 'Unnamed: 14', 'Unnamed: 15', 'Unnamed: 16', 'Unnamed: 17'] [anc RS growth at 30°C_fig.S2] 48 rows x 4 columns columns: ['P', 'replicate', 'time(hour)', 'OD600'] [adsorption assay_fig. S3] 39 rows x 3 columns columns: ['group', 'replicate', 'adsorption rate'] [one-step growth curve_fig. S4] 33 rows x 2 columns columns: ['time(mins)', 'phage density(log10)'] [kegg enrichment_fig. S5] 30 rows x 12 columns columns: ['layer1', 'layer2', 'layer3', 'ID', 'Input.number', 'Background.number', 'Gene.ratio', 'Pvalue', 'class', 'FDR', 'Input', 'Hyperlink']
The workbook has the 16 sheets the README lists (it names them 1 to 16 in the same order). Several sheets
hold a second small table to the right of the clone values: the ancestor's value (strian/strain = anc).
In the carrying capacity sheet that ancestor column is headed maxgrowthrate(OD600) though it holds the
ancestor's carrying capacity (1.44, against a growth rate of 0.05 in the growth rate sheet). The
transcriptome sheet holds three tables side by side, headed Fig 4A, fig 4B and fig 4C.
print(readme)
print("Ancestor block of each Figure 2 sheet (the columns to the right of the clone table):")
for sheet in ("resistance to phage_FIG 2A", "maximum growth rate_FIG 2B", "carrying capacity FIG2C",
"twitching_FIG2D", "biofilm_FIG 2E", "disease incidence_FIG 2F"):
print(f"--- {sheet}")
print(wb[sheet].iloc[:, 4:].dropna(how="all").to_string())
# Phage selection drives resistance-virulence trade-offs in Ralstonia solanacearum plant pathogenic bacterium irrespective of the growth temperature
[https://doi.org/10.5061/dryad.tdz08kq5d](https://doi.org/10.5061/dryad.tdz08kq5d)
This is an experiment focused on the coevolution of Ralstonia solanacearum and phages. We gathered information concerning the growth of R. solanacearum and the presence of phages at various stages of the experiment. Subsequent to isolating the evolved strains, we conducted assessments of their fundamental characteristics, encompassing growth parameters, motility, biofilm formation capability, and pathogenicity. Furthermore, we performed genomic analysis on these evolved strains, identified pivotal genes through SNP analysis, and corroborated the functionality of these genes via experimental assays.
## Description of the data and file structure
Datasets included:
1. transfer OD
• every two days’ OD600 value of each replicate in R. solanacearum -phage selection experiment
* in the row of “phage”, “1” represents in the presence of phage, and “0” represents in the absence of phage
2. phage density
* every four days’ phage density of each replicate of treatments with the presence of phage in R. solanacearum -phage selection experiment
3. scv frequency
* the densities and frequencies of SCV and mucoid ones in each replicate of the end of R. solanacearum -phage selection experiment. 4 out of 10 replicates of selection experiment were randomly chosen here.
* in the row of “phage”, “1” represents in the presence of phage, and “0” represents in the absence of phage
4. resistance to phage
* the resistance to ancestral phage of all isolates (N=128) from all treatments, and together with the phage resistance of ancestral RS-N
* the resistance to ancestral phage of each tested clone were replicated three times, and only average resistance to phage were shown here
* the “treatment” row refers to the evolution trajectory of isolates, in detail, “35\_1” represents evolved with phage at 35°C, “35\_0” represents evolved without phage at 35°C, “25\_1” represents evolved with phage at 25°C, and “25\_0” represents evolved without phage at 25°C
5. maximum growth rate and
6. carrying capacity have the same dataset structure with dataset 4)
7. twitching,
8. biofilm and
9. disease incidence were the twitching ability, biofilm formation ability and ability of causing tomato wilt disease of a subset of all isolates (N=28)
10. snp\_indel
* variant calling of a subset isolates from all treatments comparing with RS-N using Snippy
* Missing values are denoted by “n/a”
11. gene blast in uniprot
* blast variant gene sequence in UniProt (https://www.uniprot.org) to ensure the gene name and related protein
12. transcriptome
* changes in RS-N gene expression in response to phage infection
13. anc RS growth at 30°C
* RS-N density dynamics in the presence of phage and absence of phage at 30°C for 24 hours
14. adsorption assay
* adsorption rate of NNP42 phage to RS-N (anc) and phcS, gspE, gspF and pilD mutants isolated at the end of the selection experiment
15. one-step growth cuve
* One-step growth curve of NNP42 phage infecting RS-N.
16. kegg enrichment
* KEGG enrichment analysis for differentially expressed genes in the presence of phage
## Sharing/Access information
## Code/Software
The R code is employed for various tasks, including generating line plots, boxplots, and more. Additionally, it encompasses code for computing bacterial growth parameters like maximum growth rate and maximum density. Furthermore, it incorporates scripts for performing SNP analysis and creating circular genome plots.
Ancestor block of each Figure 2 sheet (the columns to the right of the clone table):
--- resistance to phage_FIG 2A
Unnamed: 4 Unnamed: 5 strian resistance to phage.1
0 NaN NaN anc 0.320821
--- maximum growth rate_FIG 2B
Unnamed: 4 Unnamed: 5 strian maxgrowthrate(OD600).1
0 NaN NaN anc 0.050785
--- carrying capacity FIG2C
Unnamed: 4 strian maxgrowthrate(OD600)
0 NaN anc 1.435515
--- twitching_FIG2D
Unnamed: 4 strain.1 twitching(mm).1
0 NaN anc 9.383333
--- biofilm_FIG 2E
Unnamed: 4 Unnamed: 5 strain.1 biofilm(OD570).1
0 NaN NaN anc 0.488667
--- disease incidence_FIG 2F
Unnamed: 4 strain.1 Disease incidence.1
0 NaN anc 0.722222
Each Figure 2 sheet holds one ancestor row, so the ancestor's value used below is that row (the helper checks there is exactly one). The README says the resistance values are means of three replicate assays and lists the sheets; it says nothing about the clones of the twitching sheet, the units of the growth rate or the blank correction of the OD600 values.
od = pt.transfer_od(wb)
print("Figure 1B sheet: OD600 values per treatment (temperature, phage) and day")
print(od.pivot_table(index=["temperature", "phage"], columns="day", values="od", aggfunc="size").to_string())
print("\nOD600 range:", od.od.min(), "to", od.od.max(), "; values below 0:", int((od.od < 0).sum()),
"(", od[od.od < 0].groupby(["temperature", "phage", "day"]).size().to_dict(), ")")
print("Day 0 OD600 by treatment (mean):", od[od.day == 0].groupby(["temperature", "phage"]).od.mean().round(4).to_dict())
ph = pt.phage_density(wb)
print("\nFigure 1C sheet: log10 phage density, populations per temperature and day")
print(ph.pivot_table(index="temperature", columns="day", values="ID", aggfunc="count").to_string())
print("IDs per temperature:", ph.groupby("temperature").ID.agg(["min", "max", "nunique"]).to_dict("index"))
print("log10 density range:", round(ph.log10_pfu.min(), 3), "to", round(ph.log10_pfu.max(), 3))
Figure 1B sheet: OD600 values per treatment (temperature, phage) and day
day 0 2 4 6 8 10 12
temperature phage
25 0 10 10 10 10 10 10 10
1 10 10 10 10 10 10 10
35 0 10 10 10 10 10 10 10
1 10 10 10 10 10 10 10
OD600 range: -0.0286 to 1.1579 ; values below 0: 8 ( {(25, 1, 4): 7, (25, 1, 6): 1} )
Day 0 OD600 by treatment (mean): {(25, 0): 0.0639, (25, 1): 0.0637, (35, 0): 0.0639, (35, 1): 0.0629}
Figure 1C sheet: log10 phage density, populations per temperature and day
day 4 8 12
temperature
25 10 10 10
35 10 10 10
IDs per temperature: {25: {'min': 11, 'max': 20, 'nunique': 10}, 35: {'min': 1, 'max': 10, 'nunique': 10}}
log10 density range: 5.051 to 11.686
The Methods describe 40 populations (four treatments, ten replicates) with OD600 measured before every serial transfer, and phage densities at every second transfer in the 20 phage populations; the README says OD every two days and phage density every four days. The data agree: 40 populations at seven days (day 0, the inoculation, then days 2 to 12), and 20 phage populations at days 4, 8 and 12. The ANOVA of Table S1a has six time points (time df 5); the claims file names days 2 to 12 (fig1b-temperature-density). Eight OD values at 25 degrees C with phage (days 4 and 6) are below 0 (down to -0.029), so the values are blank-corrected, which neither the README nor the Methods says. Phage IDs 1 to 20 cannot be matched to the OD sheet's replicates 1 to 10; no analysis here needs that. The phage density column is in log10 units as its header says (PFU per ml in the Methods); the Methods give the inoculation as 10^4 PFU/ml.
cc = pt.colony_counts(wb)
print("Figure 1D sheet (colony counts at day 12, CFU/ml):")
print(cc.to_string(index=False))
print("\nTotal = SCV + mucoid in every row:", bool((cc.scv + cc.mucoid == cc.total).all()))
print("Stored proportions (cached formula values) equal the counts' ratios:",
bool(np.allclose(cc.p_scv, cc.scv / cc.total) and np.allclose(cc.p_m, 1 - cc.p_scv)))
Figure 1D sheet (colony counts at day 12, CFU/ml): temp phage replicate scv mucoid total p_scv p_m 35 1 1 2870000 0 2870000 1.000000 0.000000 35 1 2 1032000 0 1032000 1.000000 0.000000 35 1 3 3370000 0 3370000 1.000000 0.000000 35 1 4 2120000 0 2120000 1.000000 0.000000 35 0 1 0 35000 35000 0.000000 1.000000 35 0 2 1700000 0 1700000 1.000000 0.000000 35 0 3 0 570000 570000 0.000000 1.000000 35 0 4 1200000 300000 1500000 0.800000 0.200000 25 1 1 23000 0 23000 1.000000 0.000000 25 1 2 71000 700 71700 0.990237 0.009763 25 1 3 9050 0 9050 1.000000 0.000000 25 1 4 21000 0 21000 1.000000 0.000000 25 0 1 1720000 150000 1870000 0.919786 0.080214 25 0 2 1450000 350000 1800000 0.805556 0.194444 25 0 3 1640000 250000 1890000 0.867725 0.132275 25 0 4 1880000 490000 2370000 0.793249 0.206751 Total = SCV + mucoid in every row: True Stored proportions (cached formula values) equal the counts' ratios: True
clone_sheets = {"R": "resistance to phage_FIG 2A", "mu": "maximum growth rate_FIG 2B", "A": "carrying capacity FIG2C"}
clones, ancestor = {}, {}
for v, sheet in clone_sheets.items():
clones[v], ancestor[v] = pt.clone_table(wb, sheet, v)
c128 = clones["R"][["treatment", "strain", "ctype", "temp", "phage", "name_temp", "name_rep", "name_type", "unit"]].copy()
for v in clone_sheets:
assert list(clones[v].strain) == list(c128.strain) and list(clones[v].ctype) == list(c128.ctype)
c128[v] = clones[v][v].values
print("Figure 2A to 2C sheets: the same", len(c128), "strains in the same order with the same morphotype in all three")
print("Ancestor values:", {k: round(x, 4) for k, x in ancestor.items()})
print("\nClones per treatment and colony type:")
print(pd.crosstab(c128.treatment, c128.ctype).to_string())
print("\nStrain name agrees with the treatment's temperature:", bool((c128.temp == c128.name_temp).all()),
"; with the morphotype column:", bool((c128.ctype == c128.name_type).all()))
tubes = c128.groupby(["treatment", "name_rep"]).ctype.agg(lambda s: "+".join(sorted(set(s)))).unstack()
print("\nColony types among the clones of each replicate tube (from the strain names):")
print(tubes.to_string())
print("\nUnits (replicate tube and colony type):", c128.unit.nunique())
for v in clone_sheets:
print(f"{v}: range {c128[v].min():.4f} to {c128[v].max():.4f}")
Figure 2A to 2C sheets: the same 128 strains in the same order with the same morphotype in all three
Ancestor values: {'R': 0.3208, 'mu': 0.0508, 'A': 1.4355}
Clones per treatment and colony type:
ctype m scv
treatment
25_0 16 16
25_1 4 28
35_0 20 12
35_1 0 32
Strain name agrees with the treatment's temperature: True ; with the morphotype column: True
Colony types among the clones of each replicate tube (from the strain names):
name_rep 1 2 3 4 5 6 7 8
treatment
25_0 NaN NaN NaN NaN m+scv m+scv m+scv m+scv
25_1 scv m+scv scv scv NaN NaN NaN NaN
35_0 NaN NaN NaN NaN m scv m m+scv
35_1 scv scv scv scv NaN NaN NaN NaN
Units (replicate tube and colony type): 22
R: range 0.1833 to 1.1187
mu: range 0.0367 to 0.0717
A: range 0.6097 to 1.5112
The 128 clones of Figures 2A to 2C are the same strains, in the same order and with the same colony type, in
all three sheets: 32 per treatment, with mucoid and SCV clones as Figure S1 describes (four of each from a
replicate holding both types, eight of one type otherwise). The tubes are numbered 1 to 4 (35 degrees C with
phage, 25 degrees C with phage) and 5 to 8 (without phage). The colony types per tube agree with the colony
counts of Figure 1D taken in order: at 35 degrees C without phage, replicates 1 and 3 are all mucoid, 2 all
SCV and 4 mixed, as tubes 5 and 7 (mucoid), 6 (SCV) and 8 (both); at 25 degrees C with phage, only replicate 2
has mucoid colonies, as only tube 2 has mucoid clones. This link of tubes to the Figure 1D replicates is
consistent with the data, not shown by an identifier. The README says resistance is the mean of three
technical replicates; R runs from 0.18 to 1.12. The growth rate column is headed maxgrowthrate(OD600); its
values (0.04 to 0.07) are consistent with OD600 per hour. Clones grouped by tube and colony type make 22
units, which with the seven terms of Table S3a leaves the 15 residual degrees of freedom of that table.
tw, anc_tw = pt.twitching(wb)
bf, anc_bf = pt.clone_table(wb, "biofilm_FIG 2E", "biofilm")
dz, anc_dz = pt.disease(wb)
print("Twitching sheet (no morphotype column):")
print(tw[["treatment", "strain", "hours", "diameter_mm", "type_as_named", "type_design"]].to_string())
print("\nAncestor: twitching", round(anc_tw, 3), "mm; biofilm", round(anc_bf, 4), "OD570; disease", round(anc_dz, 4))
print("\nStrains named in the 35_1 rows of the twitching sheet that are also its 35_0 strains:",
sorted(set(tw[tw.treatment == "35_1"].strain) & set(tw[tw.treatment == "35_0"].strain)))
print("35_1 strains in the biofilm and disease sheets:", sorted(bf[bf.treatment == "35_1"].strain),
sorted(dz[dz.treatment == "35_1"].strain))
print("25_0 strains and colony types in the biofilm sheet:", dict(zip(bf[bf.treatment == "25_0"].strain, bf[bf.treatment == "25_0"].ctype)))
print("25_0 strains and colony types in the disease sheet:", dict(zip(dz[dz.treatment == "25_0"].strain, dz[dz.treatment == "25_0"].ctype)))
print("Strain names repeated within a treatment in the twitching sheet:",
tw[tw.duplicated(["treatment", "strain"], keep=False)].groupby("treatment").strain.unique().to_dict())
print("\nBiofilm and disease sheets list the same 28 strains:", sorted(bf.strain) == sorted(dz.strain),
"; with the same colony types:",
bool((bf.set_index("strain").ctype.sort_index() == dz.set_index("strain").ctype.sort_index()).all()))
print("Biofilm/disease strains with a twitching row of the same name and treatment:",
len(set(zip(bf.treatment, bf.strain)) & set(zip(tw.treatment, tw.strain.str.replace("^([bd])", r"6\1", regex=True)))), "of 28")
print("\nClones per treatment and colony type, biofilm:")
print(pd.crosstab(bf.treatment, bf.ctype).to_string())
print("twitching, colony type by design reading:")
print(pd.crosstab(tw.treatment, tw.type_design).to_string())
print("\nRanges: twitching", tw.diameter_mm.min(), "to", round(tw.diameter_mm.max(), 3), "mm (all at",
tw.hours.unique().tolist(), "h); biofilm", bf.biofilm.min(), "to", bf.biofilm.max(),
"; disease", dz.incidence.min(), "to", dz.incidence.max())
k = dz.incidence * 18
print("Disease incidence x 18 plants (3 blocks x 6) is a whole number for every clone:",
bool(np.allclose(k, k.round(), atol=0.01)), "; wilted plants:", sorted(k.round().astype(int).tolist()))
Twitching sheet (no morphotype column):
treatment strain hours diameter_mm type_as_named type_design
0 35_0 b5a2 72 7.033333 m m
1 35_0 b7a2 72 9.166667 m m
2 35_0 b8a1 72 9.700000 m m
3 35_0 b8a2 72 8.350000 m m
4 35_0 6b6s1 72 7.183333 scv scv
5 35_0 6b6s2 72 11.166667 scv scv
6 35_0 6b8s1 72 9.666667 scv scv
7 35_0 6b8s2 72 11.283333 scv scv
8 35_1 6b6s1 72 9.400000 scv scv
9 35_1 6b6s2 72 9.316667 scv scv
10 35_1 6b8s1 72 8.083333 scv scv
11 35_1 6b8s2 72 10.716667 scv scv
12 25_0 6d5a2 72 11.333333 m m
13 25_0 6d7a2 72 8.750000 m m
14 25_0 6d6a2 72 8.750000 m m
15 25_0 6d8a2 72 9.400000 m m
16 25_0 6d5a2 72 8.000000 m scv
17 25_0 6d7a2 72 9.650000 m scv
18 25_0 6d6a2 72 7.466667 m scv
19 25_0 6d8a2 72 8.750000 m scv
20 25_1 6d2a1 72 5.783333 m m
21 25_1 6d2a2 72 6.633333 m m
22 25_1 6d2a3 72 6.300000 m m
23 25_1 6d2a4 72 6.250000 m m
24 25_1 6d1s1 72 9.066667 scv scv
25 25_1 6d2s2 72 10.516667 scv scv
26 25_1 6d2s3 72 10.783333 scv scv
27 25_1 6d3s2 72 11.033333 scv scv
Ancestor: twitching 9.383 mm; biofilm 0.4887 OD570; disease 0.7222
Strains named in the 35_1 rows of the twitching sheet that are also its 35_0 strains: ['6b6s1', '6b6s2', '6b8s1', '6b8s2']
35_1 strains in the biofilm and disease sheets: ['6b1s1', '6b1s2', '6b2s3', '6b2s4'] ['6b1s1', '6b1s2', '6b2s3', '6b2s4']
25_0 strains and colony types in the biofilm sheet: {'6d5a2': 'm', '6d7a2': 'm', '6d6a2': 'm', '6d8a2': 'm', '6d5s2': 'scv', '6d6s2': 'scv', '6d7s2': 'scv', '6d8s2': 'scv'}
25_0 strains and colony types in the disease sheet: {'6d5a2': 'm', '6d7a2': 'm', '6d6a2': 'm', '6d8a2': 'm', '6d5s2': 'scv', '6d6s2': 'scv', '6d7s2': 'scv', '6d8s2': 'scv'}
Strain names repeated within a treatment in the twitching sheet: {'25_0': <StringArray>
['6d5a2', '6d7a2', '6d6a2', '6d8a2']
Length: 4, dtype: str}
Biofilm and disease sheets list the same 28 strains: True ; with the same colony types: True
Biofilm/disease strains with a twitching row of the same name and treatment: 20 of 28
Clones per treatment and colony type, biofilm:
ctype m scv treatment 25_0 4 4 25_1 4 4 35_0 4 4 35_1 0 4 twitching, colony type by design reading: type_design m scv treatment 25_0 4 4 25_1 4 4 35_0 4 4 35_1 0 4 Ranges: twitching 5.78333333333333 to 11.333 mm (all at [72] h); biofilm 0.576333333333333 to 2.909 ; disease 0.0 to 0.5556 Disease incidence x 18 plants (3 blocks x 6) is a whole number for every clone: True ; wilted plants: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 5, 7, 8, 8, 10, 10]
The Methods describe 29 clones for twitching, biofilm and virulence: four mucoid and four SCV clones from each
treatment, except four SCV clones only at 35 degrees C with phage, and the ancestor. The biofilm and disease
sheets hold that design, with the same 28 strain names and colony types. The twitching sheet does not
identify its clones: it has no morphotype column; its four 35_1 rows carry the names of four 35_0 SCV clones
(6b6s1, 6b6s2, 6b8s1, 6b8s2), where the other sheets name 6b1s1, 6b1s2, 6b2s3 and 6b2s4, and their values
differ from those of the 35_0 rows; and in 25_0 the four mucoid names appear twice with different values,
where the other sheets have four SCVs (6d5s2 to 6d8s2). The treatment column is taken as right, so the 35_1
rows are read as the 35_1 SCVs. The second 25_0 block is read as the four SCVs (type_design, the
reading that gives the design and the 21 residual degrees of freedom of Table S3d); reading every row by its
name (type_as_named) is shown as a variant. The twitching values (5.8 to 11.3) are in mm as the header
says. Disease incidence is a proportion (0 to 0.72), and each value times 18 is a whole number of wilted
plants, consistent with three blocks of six plants pooled per clone.
var = pt.variants(wb)
names = pt.gene_names(wb)
print("snp_indel sheet:", len(var), "rows;", var.LOCUS_TAG.notna().sum(), "in a gene")
print("Groups and the largest number of clones carrying a variant:")
print(var.groupby("group").num_samples_with_snp.agg(["size", "max"]).to_string())
print("\nEFFECT kinds:", var.EFFECT.str.split(" ").str[0].value_counts(dropna=False).to_dict())
print("FTYPE of rows with an EFFECT:", var[var.EFFECT.notna()].FTYPE.value_counts().to_dict())
print("\nVariants summed over clones, by group, replicon and type (Table S4 counts deletions, insertions and SNPs):")
s4 = var.pivot_table(index=["group", "CHROM"], columns="TYPE", values="num_samples_with_snp", aggfunc="sum", fill_value=0)
print(s4.to_string())
print("\nGene names: the GENE column (PROKKA) against the gene_blast sheet for the genes of Table 1:")
tagged = var.dropna(subset=["LOCUS_TAG"]).drop_duplicates("LOCUS_TAG")
print(pd.DataFrame({"locus": tagged.LOCUS_TAG, "GENE column": tagged.GENE, "gene_blast name": tagged.LOCUS_TAG.map(names)}).to_string(index=False))
print("\nVariants in the genes of Table 1 named in the claims and in RSc3091, with their positions:")
t1 = var.assign(name=var.LOCUS_TAG.map(names))
t1 = t1[t1.name.isin(["phcS", "phcR", "pilD", "gspE", "gspF", "RSc3091"])]
print(t1[["group", "CHROM", "POS", "TYPE", "REF", "ALT", "LOCUS_TAG", "name", "num_samples_with_snp"]].to_string())
snp_indel sheet: 138 rows; 81 in a gene
Groups and the largest number of clones carrying a variant:
size max
group
25,+P,mucoid 20 4
25,+P,scv 21 4
25,-P,mucoid 19 4
25,-P,scv 22 4
35,+P,scv 19 4
35,-P,mucoid 18 4
35,-P,scv 19 4
EFFECT kinds: {nan: 57, 'frameshift_variant': 46, 'conservative_inframe_insertion': 13, 'disruptive_inframe_insertion': 8, 'frameshift_variant&start_lost': 6, 'missense_variant': 4, 'disruptive_inframe_deletion': 4}
FTYPE of rows with an EFFECT: {'CDS': 81}
Variants summed over clones, by group, replicon and type (Table S4 counts deletions, insertions and SNPs):
TYPE complex del ins snp
group CHROM
25,+P,mucoid RS-N.Chr1 0 5 34 5
RS-N.Chr2 0 1 24 0
25,+P,scv RS-N.Chr1 0 11 35 4
RS-N.Chr2 0 0 24 0
25,-P,mucoid RS-N.Chr1 0 8 34 4
RS-N.Chr2 0 0 24 0
25,-P,scv RS-N.Chr1 1 8 32 5
RS-N.Chr2 0 0 22 0
35,+P,scv RS-N.Chr1 0 6 30 7
RS-N.Chr2 0 0 17 0
35,-P,mucoid RS-N.Chr1 0 4 32 4
RS-N.Chr2 0 1 24 0
35,-P,scv RS-N.Chr1 0 8 32 5
RS-N.Chr2 0 0 24 0
Gene names: the GENE column (PROKKA) against the gene_blast sheet for the genes of Table 1:
locus GENE column gene_blast name
RS-N_01894 mdtD RSc1365
RS-N_02118 dhmA RSc1770
RS-N_02189 NaN RSc1701
RS-N_00287 cdiA RSc3188
RS-N_03355 yccS RSc2889
RS-N_03435 gltB gltB
RS-N_05118 NaN RSc0104
RS-N_03706 NaN RSp0796
RS-N_03478 nanT_2 tISRso1
RS-N_04190 NaN RSc1775
RS-N_03191 sasA_8 phcS
RS-N_03192 sasA_9 phcR
RS-N_00358 epsF_1 gspF
RS-N_00359 hxcR gspE
RS-N_02234 NaN RSc1658
RS-N_01695 yejB RSc1174
RS-N_03289 outO pilD
RS-N_03094 bprV_1 RSc2653
RS-N_00383 NaN RSc3091
Variants in the genes of Table 1 named in the claims and in RSc3091, with their positions:
group CHROM POS TYPE REF ALT LOCUS_TAG name num_samples_with_snp
24 35,-P,scv RS-N.Chr1 3415554 snp G A RS-N_03191 phcS 1
25 35,-P,scv RS-N.Chr1 3415879 del ATTCCAGCTC A RS-N_03191 phcS 3
42 35,+P,scv RS-N.Chr1 3415879 del ATTCCAGCTC A RS-N_03191 phcS 2
43 35,+P,scv RS-N.Chr1 3417086 snp G A RS-N_03192 phcR 1
48 35,+P,scv RS-N.Chr1 405903 del AG A RS-N_00358 gspF 1
49 35,+P,scv RS-N.Chr1 406950 snp C T RS-N_00359 gspE 2
85 25,-P,scv RS-N.Chr1 3415879 del ATTCCAGCTC A RS-N_03191 phcS 3
105 25,+P,mucoid RS-N.Chr1 3523802 ins A ATGGGCGCGC RS-N_03289 pilD 2
125 25,+P,scv RS-N.Chr1 3415879 del ATTCCAGCTC A RS-N_03191 phcS 4
130 25,+P,scv RS-N.Chr1 431461 del AT A RS-N_00383 RSc3091 1
The variant sheet is not per clone: each row is a variant in one of the seven groups of four sequenced clones,
with the number of the group's clones that carry it, so which clone carries which variant cannot be read.
The largest count in each group is 4, consistent with four clones per group (28 in all, as in the Methods).
Summed over clones, deletions, insertions and SNPs per group and replicon are those of Table S4 for every
group (for example 25 degrees C without phage, mucoid: chromosome 8, 34 and 4; megaplasmid 24 insertions);
one complex variant in the 25 degrees C -P SCV group is not in Table S4. The GENE column holds PROKKA's names,
which for the genes of Table 1 are not the paper's: phcS is sasA_8, phcR sasA_9, pilD outO, gspE
hxcR, gspF epsF_1. The dataset's gene_blast sheet gives the paper's names by locus tag, and those are used.
The positions of these variants are those of Table 1 (phcS C-3415879, phcR C-3417086, pilD C-3523802, gspE
C-406950, gspF C-405903), except the RSc3091 (RS-N_00383) deletion, which Table 1 places at C-432461 and the
variant table at 431461.
full, fig4b, fig4c = pt.transcriptome(wb)
print("Figure 4A table:", len(full), "genes, ids", full.index[0], "to", full.index[-1], "; unique:", full.index.is_unique)
print("Columns:", full.columns.tolist())
print("Missing values:", full.isna().sum().to_dict())
print("log2FoldChange range:", round(full.log2FoldChange.min(), 3), "to", round(full.log2FoldChange.max(), 3),
"; padj range:", full.padj.min(), "to", round(full.padj.max(), 4))
print("\nFigure 4B list in the sheet:", len(fig4b), "genes")
print(fig4b.to_string())
print("\nFigure 4C list in the sheet:", len(fig4c), "genes")
print(fig4c[["gene", "label", "log2FoldChange", "padj", "regulation"]].to_string())
both = pd.concat([fig4b.set_index("gene")[["log2FoldChange", "padj"]], fig4c.set_index("gene")[["log2FoldChange", "padj"]]])
both = both.join(full[["log2FoldChange", "padj"]].add_suffix("_4A"))
print("\nGenes of 4B and 4C whose log2FoldChange equals the 4A table's (to 1e-6):",
int(np.isclose(both.log2FoldChange, both.log2FoldChange_4A, atol=1e-6).sum()), "of", len(both))
kegg = wb["kegg enrichment_fig. S5"]
print("\nKEGG sheet:", len(kegg), "pathways; columns", kegg.columns.tolist())
ads = wb["adsorption assay_fig. S3"].rename(columns={"adsorption rate": "rate"})
print("\nAdsorption sheet: groups and replicates")
print(ads.groupby("group").rate.agg(["size", "min", "max"]).to_string())
Figure 4A table: 5146 genes, ids RS-N000001 to RS-N005150 ; unique: True
Columns: ['baseMean', 'log2FoldChange', 'lfcSE', 'stat', 'pvalue', 'padj']
Missing values: {'baseMean': 0, 'log2FoldChange': 1, 'lfcSE': 1, 'stat': 1, 'pvalue': 32, 'padj': 32}
log2FoldChange range: -5.942 to 4.352 ; padj range: 0.0 to 0.999
Figure 4B list in the sheet: 11 genes
gene label log2FoldChange padj
0 RS-N003355 0_RSc2889 0.291628 7.636127e-02
1 RS-N002118 1_RSc1770 -0.296235 1.238303e-01
2 RS-N002189 2_RSc1771 0.317011 4.940984e-01
3 RS-N001894 3_RSc1365 1.089598 7.891401e-10
4 RS-N003191 4_phcS 0.411209 3.043844e-03
5 RS-N003289 5_pilD -1.348558 6.525955e-36
6 RS-N003094 6_RSc2653 -0.451678 5.067883e-03
7 RS-N000359 7_gspE -0.048754 8.598386e-01
8 RS-N000383 10_RSc3091 0.057861 7.856951e-01
9 RS-N000358 11_gspF 0.488300 8.604830e-03
10 RS-N003192 12_phcR 1.834769 8.035238e-28
Figure 4C list in the sheet: 20 genes
gene label log2FoldChange padj regulation
0 RS-N004717 RS-N004717_JetC_II -4.209067 6.918707e-65 DOWN
1 RS-N000213 RS-N000213_Type_I_MTases 2.302917 2.167891e-54 UP
2 RS-N000216 RS-N000216_Effector_1TM_DUF1043 2.324671 1.498133e-35 UP
3 RS-N000215 RS-N000215_Type_I_REases 2.600479 2.436394e-26 UP
4 RS-N000217 RS-N000217_Cyclase_SMODS 2.640740 1.763284e-44 UP
5 RS-N003916 type II and III secretion system protein family protein _RS-N003916_cpaC, rcpA -4.405639 5.800000e-186 DOWN
6 RS-N003912 type II secretion system F family protein_RS-N003912_tadC -4.266422 1.790000e-130 DOWN
7 RS-N003622 type III secretion protein_RS-N003622_hrpB2 -2.490343 1.790000e-17 DOWN
8 RS-N003623 type III secretion system apparatus protein_RS-N003623_hrpB1 -2.371315 2.050000e-10 DOWN
9 RS-N003619 type III secretion system protein_RS-N003619_yscL, sctL -2.279596 1.430000e-10 DOWN
10 RS-N003629 type III secretion system export apparatus protein EscS/YscS/HrcS family_RS-N003629_yscS, sctS, hrcS -2.065348 4.660000e-05 DOWN
11 RS-N003625 type III secretion system export apparatus protein EscV/YscV/HrcV family_RS-N003625_yscV, sctV, hrcV -2.055090 1.210000e-19 DOWN
12 RS-N003628 type III secretion system export apparatus protein EscR/YscR/HrcR family _RS-N003628_yscR, sctR, hrcR -2.048767 1.190000e-07 DOWN
13 RS-N003627 type III secretion apparatus protein YscQ/HrcQ family_RS-N003627_yscQ, sctQ, hrcQ -2.001040 9.110000e-09 DOWN
14 RS-N004330 type II secretion system protein_RS-N004330_gspG 2.037027 1.527727e-03 UP
15 RS-N003919 prepilin peptidase_RS-N003919_cpaA, tadV -5.942366 0.000000e+00 DOWN
16 RS-N003920 Flp family type IVb pilin_RS-N003920_flp, pilA -5.625140 8.380000e-273 DOWN
17 RS-N003915 pilus assembly protein_RS-N003915_cpaE, tadZ -4.738090 8.540000e-126 DOWN
18 RS-N003917 Flp pilus assembly protein_RS-N003917_cpaB, rcpC -4.396537 1.730000e-124 DOWN
19 RS-N004603 Sigma-fimbriae chaperone protein_RS-N004603_fimC -2.100763 3.900000e-18 DOWN
Genes of 4B and 4C whose log2FoldChange equals the 4A table's (to 1e-6): 31 of 31
KEGG sheet: 30 pathways; columns ['layer1', 'layer2', 'layer3', 'ID', 'Input.number', 'Background.number', 'Gene.ratio', 'Pvalue', 'class', 'FDR', 'Input', 'Hyperlink']
Adsorption sheet: groups and replicates
size min max
group
anc 3 0.750 0.815
gspE_B+S4 3 0.625 0.840
gspF_B+S1 3 0.630 0.710
phcS_B+S2 3 0.665 0.915
phcS_B+S3 3 0.695 0.775
phcS_B-S1 3 0.750 0.825
phcS_B-S2 3 0.795 0.860
phcS_D+S2 3 0.745 0.905
phcS_D+S3 3 0.805 0.900
phcS_D-S2 3 0.810 0.875
phcS_D-S4 3 0.775 0.800
pilD_D+a1 3 0.175 0.395
pilD_D+a2 3 0.395 0.425
The DESeq2 table has the columns DESeq2 writes (baseMean, log2FoldChange, lfcSE, stat, pvalue, padj) for 5146 genes; padj is missing for 32 of them. The Figure 4B and 4C lists beside it carry the same fold changes as the main table. The Figure 4B list holds 11 genes, numbered 0 to 7 and 10 to 12, where the paper speaks of 13 mutated genes. The adsorption sheet holds the ancestor and twelve mutants, three replicates each, as proportions (0.18 to 0.92), where Figure S3 speaks of percentages.
outside = {}
r = requests.get("https://www.ebi.ac.uk/ena/portal/api/filereport", timeout=120, params={
"accession": "PRJNA1032726", "result": "read_run", "fields": "run_accession,sample_alias,library_strategy,fastq_bytes"})
r.raise_for_status()
runs = pd.read_csv(pd.io.common.StringIO(r.text), sep="\t").sort_values("run_accession")
runs["gb"] = runs.fastq_bytes.astype(str).str.split(";").apply(lambda x: sum(int(b) for b in x)) / 1e9
print("PRJNA1032726 (ENA read runs):")
print(runs[["run_accession", "sample_alias", "library_strategy", "gb"]].round(2).to_string(index=False))
print(f"total {runs.gb.sum():.1f} GB of gzipped FASTQ")
for db in ("sra", "assembly", "biosample"):
r = requests.get("https://eutils.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi", timeout=120,
params={"dbfrom": "bioproject", "db": db, "id": "1028400", "retmode": "json"})
r.raise_for_status()
sets = r.json()["linksets"][0].get("linksetdbs", [])
outside[db] = len(sets[0]["links"]) if sets else 0
print("\nPRJNA1028400 (NCBI BioProject 'Genome sequences of evolved Ralstonia solanacearum'): linked records", outside)
PRJNA1032726 (ENA read runs): run_accession sample_alias library_strategy gb SRR26541645 P5 RNA-Seq 0.84 SRR26541646 P4 RNA-Seq 0.94 SRR26541647 P3 RNA-Seq 1.00 SRR26541648 P2 RNA-Seq 1.06 SRR26541649 P1 RNA-Seq 0.93 SRR26541650 R5 RNA-Seq 0.93 SRR26541651 R4 RNA-Seq 0.91 SRR26541652 R3 RNA-Seq 0.63 SRR26541653 R2 RNA-Seq 0.62 SRR26541654 R1 RNA-Seq 1.27 total 9.1 GB of gzipped FASTQ
PRJNA1028400 (NCBI BioProject 'Genome sequences of evolved Ralstonia solanacearum'): linked records {'sra': 0, 'assembly': 14, 'biosample': 28}
PRJNA1032726 holds the RNA-seq: ten paired-end runs, five phage (P1 to P5) and five control (R1 to R5)
samples, as in the Methods. PRJNA1028400 holds the genome project: 28 BioSamples and 14 assemblies, and no
sequencing runs, so the reads of the 28 clones are not public. The paper's genome and transcriptome claims are
data: dataset in the claims file and are tested on the processed tables the dataset holds (the variant calls, the DESeq2 table and the KEGG table);
the pipelines from reads were not re-run (see each claim).
Claims¶
fig1b-temperature-density¶
R. solanacearum reached higher population densities at 35 degrees C than at 25 degrees C during the selection experiment. (Table S1a: temperature F(1, 36) = 73.989, p < .001)
Criterion: a repeated-measures ANOVA of OD600 with population as subject, temperature, phage and time (a factor, df 5) and all their interactions, on days 2, 4, 6, 8, 10 and 12, the densities measured before each transfer; day 0 is the inoculation. Reproduced when the temperature effect has p < 0.05 and mean OD600 is higher at 35 degrees C. The data hold seven time points; the claims file named days 2 to 12 in answer to this attempt's question. Days 0 to 10 and all seven days were the other readings the first run showed; they are kept as variants and decide nothing.
readings_days = {"criterion, days 2-12": [2, 4, 6, 8, 10, 12], "variant, days 0-10": [0, 2, 4, 6, 8, 10],
"variant, all seven days": [0, 2, 4, 6, 8, 10, 12]}
rm = {}
for label, days in readings_days.items():
rm[label] = pt.rm_anova(od[od.day.isin(days)], "od", "subject", ["temperature", "phage"], "day")
print(f"--- {label}")
print(rm[label].round(4).to_string(index=False))
print("\nTable S1a: temperature SS 5.681 F 73.989; phage SS 3.747 F 48.800; residual (subjects) SS 2.764 df 36;"
"\ntime SS 0.740 F 7.316; temperature x phage x time SS 0.391 F 3.870; residual (within) SS 3.639 df 180")
--- criterion, days 2-12 stratum term df ss ms F df_resid p subject C(temperature) 1.0 5.6807 5.6807 73.9893 36.0 0.0000 subject C(phage) 1.0 3.7467 3.7467 48.7997 36.0 0.0000 subject C(temperature):C(phage) 1.0 0.0139 0.0139 0.1816 36.0 0.6726 subject Residual 36.0 2.7640 0.0768 NaN NaN NaN within C(day) 5.0 0.7395 0.1479 7.3155 180.0 0.0000 within C(day):C(temperature) 5.0 0.2886 0.0577 2.8553 180.0 0.0166 within C(day):C(phage) 5.0 0.9429 0.1886 9.3268 180.0 0.0000 within C(day):C(temperature):C(phage) 5.0 0.3912 0.0782 3.8696 180.0 0.0024 within Residual 180.0 3.6393 0.0202 NaN NaN NaN --- variant, days 0-10 stratum term df ss ms F df_resid p subject C(temperature) 1.0 3.6872 3.6872 70.9107 36.0 0.0000 subject C(phage) 1.0 3.0061 3.0061 57.8115 36.0 0.0000 subject C(temperature):C(phage) 1.0 0.0091 0.0091 0.1747 36.0 0.6785 subject Residual 36.0 1.8719 0.0520 NaN NaN NaN within C(day) 5.0 4.4277 0.8855 52.2313 180.0 0.0000 within C(day):C(temperature) 5.0 0.9978 0.1996 11.7709 180.0 0.0000 within C(day):C(phage) 5.0 1.4372 0.2874 16.9540 180.0 0.0000 within C(day):C(temperature):C(phage) 5.0 0.1241 0.0248 1.4638 180.0 0.2039 within Residual 180.0 3.0517 0.0170 NaN NaN NaN --- variant, all seven days stratum term df ss ms F df_resid p subject C(temperature) 1.0 4.8669 4.8669 74.0189 36.0 0.0000 subject C(phage) 1.0 3.2141 3.2141 48.8810 36.0 0.0000 subject C(temperature):C(phage) 1.0 0.0119 0.0119 0.1803 36.0 0.6737 subject Residual 36.0 2.3671 0.0658 NaN NaN NaN within C(day) 6.0 4.8276 0.8046 43.0559 216.0 0.0000 within C(day):C(temperature) 6.0 1.1025 0.1837 9.8324 216.0 0.0000 within C(day):C(phage) 6.0 1.4755 0.2459 13.1599 216.0 0.0000 within C(day):C(temperature):C(phage) 6.0 0.3933 0.0655 3.5074 216.0 0.0025 within Residual 216.0 4.0365 0.0187 NaN NaN NaN Table S1a: temperature SS 5.681 F 73.989; phage SS 3.747 F 48.800; residual (subjects) SS 2.764 df 36; time SS 0.740 F 7.316; temperature x phage x time SS 0.391 F 3.870; residual (within) SS 3.639 df 180
for label, days in readings_days.items():
sub = od[od.day.isin(days)]
m = sub.groupby("subject").agg(temperature=("temperature", "first"), phage=("phage", "first"), od=("od", "mean"))
print(f"{label}: mean OD600 over populations and time points, 25 C {m[m.temperature == 25].od.mean():.4f}, "
f"35 C {m[m.temperature == 35].od.mean():.4f}; by phage at 25 C",
m[m.temperature == 25].groupby("phage").od.mean().round(4).to_dict(), "at 35 C",
m[m.temperature == 35].groupby("phage").od.mean().round(4).to_dict())
t = rm[label].set_index("term")
print(f" temperature F = {t.loc['C(temperature)', 'F']:.3f}, p = {t.loc['C(temperature)', 'p']:.2e}")
fig, ax = plt.subplots(figsize=(6, 3.5))
for (temp, phage), g in od.groupby(["temperature", "phage"]):
s = g.groupby("day").od.agg(["mean", "sem"])
ax.errorbar(s.index, s["mean"], s["sem"], marker="o", label=f"{temp} C, {'+P' if phage else '-P'}")
ax.set_xlabel("day"); ax.set_ylabel("OD600"); ax.legend(fontsize=8); ax.set_title("Figure 1B, mean +- SEM, n = 10")
plt.show()
print("Plotted means range from", round(od.groupby(["temperature", "phage", "day"]).od.mean().min(), 4), "to",
round(od.groupby(["temperature", "phage", "day"]).od.mean().max(), 4), "; no point falls outside the axes")
criterion, days 2-12: mean OD600 over populations and time points, 25 C 0.2550, 35 C 0.5627; by phage at 25 C {0: 0.3876, 1: 0.1225} at 35 C {0: 0.6801, 1: 0.4454}
temperature F = 73.989, p = 2.96e-10
variant, days 0-10: mean OD600 over populations and time points, 25 C 0.2202, 35 C 0.4681; by phage at 25 C {0: 0.326, 1: 0.1144} at 35 C {0: 0.5862, 1: 0.35}
temperature F = 70.911, p = 4.96e-10
variant, all seven days: mean OD600 over populations and time points, 25 C 0.2277, 35 C 0.4914; by phage at 25 C {0: 0.3414, 1: 0.1141} at 35 C {0: 0.592, 1: 0.3908}
temperature F = 74.019, p = 2.94e-10
Plotted means range from -0.0082 to 0.7346 ; no point falls outside the axes
Result: reproduced. On days 2 to 12 the ANOVA reproduces every sum of squares and F of Table S1a to the printed digits (temperature F(1, 36) = 73.989, p = 2.96e-10); mean OD600 over populations and time points is 0.5627 at 35 degrees C against 0.2550 at 25 degrees C. The variants, days 0 to 10 (F = 70.911) and all seven days (F = 74.019), give the same result.
fig1b-phage-reduces-density¶
Phage presence reduced R. solanacearum densities at both temperatures. (phage F(1, 36) = 48.80, p < .001)
Criterion: the ANOVA above; reproduced when the phage effect has p < 0.05 and, at each temperature, mean OD600 over populations and time points is lower with phage.
for label, days in readings_days.items():
t = rm[label].set_index("term")
sub = od[od.day.isin(days)]
m = sub.groupby(["temperature", "phage"]).od.mean().unstack()
print(f"{label}: phage F = {t.loc['C(phage)', 'F']:.3f}, p = {t.loc['C(phage)', 'p']:.2e}; mean OD600 -P / +P:",
{T: (round(m.loc[T, 0], 4), round(m.loc[T, 1], 4), f"{100 * (1 - m.loc[T, 1] / m.loc[T, 0]):.1f}% lower") for T in (25, 35)})
criterion, days 2-12: phage F = 48.800, p = 3.43e-08; mean OD600 -P / +P: {25: (np.float64(0.3876), np.float64(0.1225), '68.4% lower'), 35: (np.float64(0.6801), np.float64(0.4454), '34.5% lower')}
variant, days 0-10: phage F = 57.811, p = 5.39e-09; mean OD600 -P / +P: {25: (np.float64(0.326), np.float64(0.1144), '64.9% lower'), 35: (np.float64(0.5862), np.float64(0.35), '40.3% lower')}
variant, all seven days: phage F = 48.881, p = 3.37e-08; mean OD600 -P / +P: {25: (np.float64(0.3414), np.float64(0.1141), '66.6% lower'), 35: (np.float64(0.592), np.float64(0.3908), '34.0% lower')}
Result: reproduced. Phage effect F(1, 36) = 48.800, p = 3.43e-08 on days 2 to 12 (variants: 57.811 on days 0 to 10, 48.881 with all seven days). Mean OD600 on days 2 to 12 is lower with phage by 68.4% at 25 degrees C (0.1225 against 0.3876) and by 34.5% at 35 degrees C (0.4454 against 0.6801), and lower with phage under both variants too.
fig1b-phage-effect-diminishes¶
The reduction of bacterial density by phage diminished over time, with a significant temperature by phage by time interaction. (F(5, 180) = 3.870, p < .01)
Criterion: the ANOVA above; reproduced when the three-way interaction has p < 0.05 and, at both temperatures, the difference in mean OD600 between -P and +P at the last time point is smaller than its largest value at an earlier time point; partially reproduced when only one holds.
diff = od.groupby(["temperature", "day", "phage"]).od.mean().unstack()
diff["minus_P_minus_plus_P"] = diff[0] - diff[1]
print("Mean OD600 by temperature and day, and the -P minus +P difference:")
print(diff.round(4).to_string())
for label, days in readings_days.items():
t = rm[label].set_index("term")
term = "C(day):C(temperature):C(phage)"
out = []
for T in (25, 35):
d = diff.loc[T, "minus_P_minus_plus_P"].loc[days]
out.append(f"{T} C: last {d.iloc[-1]:.4f} vs largest earlier {d.iloc[:-1].max():.4f} (day {d.iloc[:-1].idxmax()})")
print(f"{label}: three-way F({t.loc[term, 'df']:.0f}, {t.loc[term, 'df_resid']:.0f}) = {t.loc[term, 'F']:.3f}, "
f"p = {t.loc[term, 'p']:.4f}; " + "; ".join(out))
Mean OD600 by temperature and day, and the -P minus +P difference:
phage 0 1 minus_P_minus_plus_P
temperature day
25 0 0.0639 0.0637 0.0002
2 0.4628 0.0266 0.4362
4 0.3683 -0.0082 0.3765
6 0.3664 0.0617 0.3046
8 0.3249 0.2756 0.0493
10 0.3693 0.2670 0.1023
12 0.4339 0.1120 0.3219
35 0 0.0639 0.0629 0.0010
2 0.6575 0.1529 0.5046
4 0.6746 0.3721 0.3026
6 0.6970 0.4619 0.2352
8 0.7346 0.4917 0.2429
10 0.6893 0.5587 0.1306
12 0.6274 0.6353 -0.0080
criterion, days 2-12: three-way F(5, 180) = 3.870, p = 0.0024; 25 C: last 0.3219 vs largest earlier 0.4362 (day 2); 35 C: last -0.0080 vs largest earlier 0.5046 (day 2)
variant, days 0-10: three-way F(5, 180) = 1.464, p = 0.2039; 25 C: last 0.1023 vs largest earlier 0.4362 (day 2); 35 C: last 0.1306 vs largest earlier 0.5046 (day 2)
variant, all seven days: three-way F(6, 216) = 3.507, p = 0.0025; 25 C: last 0.3219 vs largest earlier 0.4362 (day 2); 35 C: last -0.0080 vs largest earlier 0.5046 (day 2)
Result: reproduced. On days 2 to 12 the three-way interaction is F(5, 180) = 3.870, p = 0.0024, as in Table S1a. The -P minus +P difference at day 12 (0.3219 at 25 degrees C, -0.0080 at 35 degrees C) is smaller than its largest earlier value (0.4362 and 0.5046, both on day 2). Of the variants, days 0 to 10 give an interaction that is not significant (F = 1.464, p = 0.2039), while the differences still decline (0.1023 and 0.1306 on day 10), and all seven days give p = 0.0025; the criterion names days 2 to 12, so the variants decide nothing. The criterion compares only the last day with the largest earlier value. At 25 degrees C the difference does not fall steadily: it drops to 0.049 on day 8 and rises again to 0.322 on day 12 (table above).
auc-phage-reduction-by-temperature¶
Over the whole experiment, phages reduced bacterial density (area under the density curve) more at 25 degrees C (66.8%) than at 35 degrees C (37.18%), with a significant phage effect at each temperature.
Criterion: the area under each population's OD600 curve, compared between +P and -P by one-way ANOVA at each temperature; the reduction is 1 minus the mean +P area over the mean -P area. Readings, as the claims file now names them by day: (a) trapezoidal area over days 2 to 12, (b) trapezoidal area over days 0 to 12, with the day 0 inoculation OD600, (c) the sum of OD600 over days 2 to 12. Reproduced when under every reading both reductions are within tolerance (66.8 +- 6.68, 37.18 +- 3.72), the 25 degrees C reduction is the larger and both ANOVAs give p < 0.05; blocked as methods underspecified when the readings differ. The first run's measures over days 0 to 10 and the sum over days 0 to 12 are kept as variants.
wide = od.pivot_table(index=["temperature", "phage", "replicate"], columns="day", values="od")
def trapezoid(days):
return pd.Series(np.trapezoid(wide[days].values, x=np.array(days, float), axis=1), index=wide.index)
auc_readings = {
"(a) trapezoid, days 2-12": trapezoid([2, 4, 6, 8, 10, 12]),
"(b) trapezoid, days 0-12": trapezoid([0, 2, 4, 6, 8, 10, 12]),
"(c) sum, days 2-12": wide[[2, 4, 6, 8, 10, 12]].sum(axis=1),
"variant: sum, days 0-12": wide[[0, 2, 4, 6, 8, 10, 12]].sum(axis=1),
"variant: trapezoid, days 0-10": trapezoid([0, 2, 4, 6, 8, 10]),
"variant: sum, days 0-10": wide[[0, 2, 4, 6, 8, 10]].sum(axis=1),
}
rows = []
for label, area in auc_readings.items():
row = {"reading": label}
for T in (25, 35):
plus, minus = area.loc[(T, 1)], area.loc[(T, 0)]
f, p, df1, df2 = pt.one_way(plus, minus)
row[f"reduction_{T}C_%"] = 100 * (1 - plus.mean() / minus.mean())
row[f"F_{T}C"] = f
row[f"df_{T}C"] = f"({df1}, {df2})"
row[f"p_{T}C"] = p
row["within_tolerance"] = abs(row["reduction_25C_%"] - 66.8) <= 6.68 and abs(row["reduction_35C_%"] - 37.18) <= 3.72
row["25C_larger"] = row["reduction_25C_%"] > row["reduction_35C_%"]
rows.append(row)
auc = pd.DataFrame(rows)
print(auc.to_string(index=False, float_format=lambda x: f"{x:.4g}"))
print("\nPaper: 66.8% and 37.18%; 25 C F1,19 = 88.78, 35 C F1,19 = 15.66, both p < .001")
reading reduction_25C_% F_25C df_25C p_25C reduction_35C_% F_35C df_35C p_35C within_tolerance 25C_larger
(a) trapezoid, days 2-12 64.55 63.16 (1, 18) 2.692e-07 33.73 11.26 (1, 18) 0.003516 True True
(b) trapezoid, days 0-12 66.8 88.77 (1, 18) 2.216e-08 37.18 15.66 (1, 18) 0.0009238 True True
(c) sum, days 2-12 68.4 109.8 (1, 18) 4.33e-09 34.5 12.29 (1, 18) 0.002521 True True
variant: sum, days 0-12 66.58 110.5 (1, 18) 4.113e-09 34 12.31 (1, 18) 0.002506 True True
variant: trapezoid, days 0-10 70.03 76.9 (1, 18) 6.474e-08 43.02 21.57 (1, 18) 0.0002018 False True
variant: sum, days 0-10 64.9 68.07 (1, 18) 1.576e-07 40.29 19.85 (1, 18) 0.0003057 True True
Paper: 66.8% and 37.18%; 25 C F1,19 = 88.78, 35 C F1,19 = 15.66, both p < .001
Result: reproduced. Under reading (a), the trapezoidal area over days 2 to 12, the reductions are 64.55% at 25 degrees C and 33.73% at 35 degrees C (F = 63.16 and 11.26 at df (1, 18)); under (b), the trapezoid over days 0 to 12, 66.80% and 37.18%, the paper's values, with the paper's F of 88.77 and 15.66 (the paper prints df (1, 19)); under (c), the sum over days 2 to 12, 68.40% and 34.50%. All six reductions are within tolerance, the 25 degrees C reduction is the larger under each reading, and every ANOVA has p < 0.004. Of the variants, the trapezoid over days 0 to 10 gives 43.02% at 35 degrees C, outside the tolerance; the variants decide nothing.
fig1c-phage-density-temperature¶
Phage densities were higher at 25 than at 35 degrees C, and at both temperatures rose rapidly at first and then declined over time. (Table S1b: temperature F(1, 18) = 32.75; time F(2, 36) = 86.963)
Criterion: a repeated-measures ANOVA of log10 phage density with population as subject, temperature, time (df 2) and their interaction. The data hold three phage time points, so readings (a) and (b) of the criterion are the same. Reproduced when the temperature and time effects both have p < 0.05, mean log10 density is higher at 25 degrees C, and at each temperature the mean rises from the inoculation density (10^4 PFU/ml, log10 = 4) to its highest value before the last time point and is lower at the last time point than at that highest value.
rm_phage = pt.rm_anova(ph, "log10_pfu", "ID", ["temperature"], "day")
print(rm_phage.round(4).to_string(index=False))
print("Table S1b: temperature SS 28.61 F 32.75; residual 15.73 (18); time SS 84.88 F 86.963; temp x time F 0.386; residual 17.57 (36)")
means = ph.groupby(["temperature", "day"]).log10_pfu.mean().unstack()
print("\nMean log10 PFU/ml by temperature and day (inoculation 4.0 on day 0):")
print(means.round(3).to_string())
print("Mean over populations and days: 25 C", round(ph[ph.temperature == 25].log10_pfu.mean(), 3),
"; 35 C", round(ph[ph.temperature == 35].log10_pfu.mean(), 3))
for T in (25, 35):
s = means.loc[T]
print(f"{T} C: highest {s.max():.3f} on day {s.idxmax()} (before the last day: {s.idxmax() < s.index[-1]}), "
f"above 4.0: {s.max() > 4}; last day {s.iloc[-1]:.3f}, lower than the highest: {s.iloc[-1] < s.max()}")
stratum term df ss ms F df_resid p subject C(temperature) 1.0 28.6123 28.6123 32.7459 18.0 0.0000 subject Residual 18.0 15.7278 0.8738 NaN NaN NaN within C(day) 2.0 84.8848 42.4424 86.9627 36.0 0.0000 within C(day):C(temperature) 2.0 0.3764 0.1882 0.3856 36.0 0.6828 within Residual 36.0 17.5699 0.4881 NaN NaN NaN Table S1b: temperature SS 28.61 F 32.75; residual 15.73 (18); time SS 84.88 F 86.963; temp x time F 0.386; residual 17.57 (36) Mean log10 PFU/ml by temperature and day (inoculation 4.0 on day 0): day 4 8 12 temperature 25 10.511 7.719 8.037 35 8.947 6.318 6.859 Mean over populations and days: 25 C 8.756 ; 35 C 7.375 25 C: highest 10.511 on day 4 (before the last day: True), above 4.0: True; last day 8.037, lower than the highest: True 35 C: highest 8.947 on day 4 (before the last day: True), above 4.0: True; last day 6.859, lower than the highest: True
Result: reproduced. The repeated-measures ANOVA reproduces Table S1b (temperature F(1, 18) = 32.746, time F(2, 36) = 86.963, both p < 1e-4). Mean log10 density is 8.756 at 25 degrees C and 7.375 at 35 degrees C. At both temperatures it is highest at the first sample, day 4 (10.511 and 8.947, from 4.0 at inoculation), and lower on day 12 (8.037 and 6.859). The minimum is on day 8 at both temperatures, so the decline is not monotonic.
fig1d-mucoid-lost-with-phage¶
At day 12, ancestral mucoid colonies were lost in the presence of phage at both temperatures and replaced by small colony variants (SCVs), except in one phage replicate at 25 degrees C. (0.24%)
Criterion: reproduced when no mucoid colonies are found in any 35 degrees C +P replicate, mucoid colonies are found in exactly one 25 degrees C +P replicate, SCVs are the majority in the other three, and 0.24% (+- 0.024) matches under one reading: (a) mucoid colonies as a share of all colonies of the 25 degrees C +P treatment, taken (a1) as the mean of the four replicate frequencies or (a2) from pooled counts, as the claims file now names them; (b) SCVs as a share of the colonies of the exceptional replicate.
plus = cc[cc.phage == 1]
print(plus[["temp", "replicate", "scv", "mucoid", "total", "p_scv", "p_m"]].to_string(index=False))
p35, p25 = plus[plus.temp == 35], plus[plus.temp == 25]
exc = p25[p25.mucoid > 0]
print("\n35 C +P replicates with mucoid colonies:", int((p35.mucoid > 0).sum()), "of", len(p35))
print("25 C +P replicates with mucoid colonies:", int((p25.mucoid > 0).sum()), "of", len(p25), "(replicate",
exc.replicate.tolist(), ")")
print("SCV share in the other three 25 C +P replicates:", p25[p25.mucoid == 0].p_scv.round(4).tolist())
reads = {"(a1) mucoid share of the 25 C +P treatment, mean of replicate frequencies": 100 * p25.p_m.mean(),
"(a2) mucoid share of the 25 C +P treatment, pooled counts": 100 * p25.mucoid.sum() / p25.total.sum(),
"(b) SCV share of the exceptional replicate": 100 * exc.p_scv.iloc[0]}
for k, v in reads.items():
print(f"{k}: {v:.3f}% (0.24 +- 0.024: {abs(v - 0.24) <= 0.024})")
temp replicate scv mucoid total p_scv p_m 35 1 2870000 0 2870000 1.000000 0.000000 35 2 1032000 0 1032000 1.000000 0.000000 35 3 3370000 0 3370000 1.000000 0.000000 35 4 2120000 0 2120000 1.000000 0.000000 25 1 23000 0 23000 1.000000 0.000000 25 2 71000 700 71700 0.990237 0.009763 25 3 9050 0 9050 1.000000 0.000000 25 4 21000 0 21000 1.000000 0.000000 35 C +P replicates with mucoid colonies: 0 of 4 25 C +P replicates with mucoid colonies: 1 of 4 (replicate [2] ) SCV share in the other three 25 C +P replicates: [1.0, 1.0, 1.0] (a1) mucoid share of the 25 C +P treatment, mean of replicate frequencies: 0.244% (0.24 +- 0.024: True) (a2) mucoid share of the 25 C +P treatment, pooled counts: 0.561% (0.24 +- 0.024: False) (b) SCV share of the exceptional replicate: 99.024% (0.24 +- 0.024: False)
Result: reproduced. No 35 degrees C +P replicate has mucoid colonies (0 of 4); one 25 degrees C +P replicate (replicate 2) has them, and the other three are all SCV. Under reading (a1), the mean of the four replicate frequencies, mucoid colonies are 0.244% of the 25 degrees C +P treatment, the paper's 0.24%. Under (a2), pooled counts, the share is 0.561%, and reading (b) gives 99.024%, neither within tolerance; one matching reading is what the criterion asks. The text calls the 0.24% the frequency of SCVs; in the data it is the frequency of mucoid colonies.
fig1d-phage-raises-scv-frequency¶
SCVs also increased without phage, but around 20% (25 degrees C) and 55% (35 degrees C) of colonies stayed mucoid in the -P treatments, and phage presence raised the SCV frequency. (Table S2: phage F(1, 12) = 7.016)
Criterion: a two-way ANOVA of each replicate's colony type frequency on temperature, phage and their interaction, untransformed. The treatment frequency is read (a) as the mean of the replicate frequencies, (b) from pooled counts. Reproduced when the phage effect has p < 0.05 with more SCVs under phage and both frequencies (20 +- 2 and 55 +- 5.5) are within tolerance under one reading.
from statsmodels.formula.api import ols
import statsmodels.api as sm
a2 = sm.stats.anova_lm(ols("p_scv ~ C(temp) * C(phage)", cc).fit(), typ=1)
print(a2.round(4).to_string())
print("Table S2: temperature SS 0.155 F 2.218; phage SS 0.491 F 7.016 p 0.021; interaction F 2.273; residual SS 0.840 (12)")
g = cc.groupby(["temp", "phage"])
freq = pd.DataFrame({"mean SCV %": 100 * g.p_scv.mean(), "mean mucoid %": 100 * g.p_m.mean(),
"pooled mucoid %": 100 * g.mucoid.sum() / g.total.sum()})
print("\n", freq.round(2).to_string())
for col in ("mean mucoid %", "pooled mucoid %"):
v25, v35 = freq.loc[(25, 0), col], freq.loc[(35, 0), col]
print(f"{col}: 25 C -P {v25:.2f} (20 +- 2: {abs(v25 - 20) <= 2}); 35 C -P {v35:.2f} (55 +- 5.5: {abs(v35 - 55) <= 5.5})")
print("Mean SCV share without phage", round(100 * cc[cc.phage == 0].p_scv.mean(), 2), "% and with phage",
round(100 * cc[cc.phage == 1].p_scv.mean(), 2), "%")
df sum_sq mean_sq F PR(>F)
C(temp) 1.0 0.1553 0.1553 2.2181 0.1622
C(phage) 1.0 0.4914 0.4914 7.0162 0.0212
C(temp):C(phage) 1.0 0.1592 0.1592 2.2734 0.1575
Residual 12.0 0.8404 0.0700 NaN NaN
Table S2: temperature SS 0.155 F 2.218; phage SS 0.491 F 7.016 p 0.021; interaction F 2.273; residual SS 0.840 (12)
mean SCV % mean mucoid % pooled mucoid %
temp phage
25 0 84.66 15.34 15.64
1 99.76 0.24 0.56
35 0 45.00 55.00 23.78
1 100.00 0.00 0.00
mean mucoid %: 25 C -P 15.34 (20 +- 2: False); 35 C -P 55.00 (55 +- 5.5: True)
pooled mucoid %: 25 C -P 15.64 (20 +- 2: False); 35 C -P 23.78 (55 +- 5.5: False)
Mean SCV share without phage 64.83 % and with phage 99.88 %
Result: partially reproduced. The two-way ANOVA reproduces Table S2 (phage F(1, 12) = 7.016, p = 0.0212), and the mean SCV share is 99.88% with phage against 64.83% without. At 35 degrees C without phage 55.00% of colonies are mucoid (mean of replicate frequencies), the paper's value; at 25 degrees C 15.34%, outside 20 +- 2. From pooled counts the shares are 15.64% and 23.78%, both outside.
Figure 2A to 2C: the ANOVA readings¶
The claims for Figures 2A to 2C use one ANOVA with temperature, phage, colony type and the three two-way interactions. The unit is the mean of the clones of each colony type within each replicate population (22 units), as the claims file now fixes in answer to this attempt's question; each clone as the unit is a variant. Sums of squares are read (i) sequential in the order of Table S3, (ii) type II. The cell below runs both readings on unit means and, as the variant, on clones, for the three measures; the claims read from it.
anova2 = {}
for v in clone_sheets:
units = pt.unit_means(c128, v)
for unit_label, data in (("unit means", units), ("variant: clones", c128)):
for typ in (1, 2):
anova2[(v, unit_label, typ)] = pt.clone_anova(data, v, typ)
print(f"--- {v}: unit means n = {len(units)}, clones n = {len(c128)}")
for unit_label in ("unit means", "variant: clones"):
for typ in (1, 2):
print(f"{unit_label}, {'(i) sequential' if typ == 1 else '(ii) type II'}")
print(anova2[(v, unit_label, typ)].round(5).to_string())
print("\nTable S3a (R): temperature F 6.459 p .023; phage 280.752; colony 445.207; TxP 2.099 p .168; TxC 6.516; PxC 86.063; residual SS 0.043 (15)")
print("Table S3b (mu): temperature 26.621; phage 3.226; colony 29.393; TxP 4.360; TxC 0.002 p .967; PxC 2.898 p .109; residual 1.109e-4 (15)")
print("Table S3c (A): temperature 4.666; phage 19.598; colony 87.965; TxP 3.268; TxC 3.222 p .093; PxC 2.900 p .109; residual 0.032 (15)")
def term_table(v, term):
return pd.DataFrame([{"reading": f"{u}, {'(i)' if t == 1 else '(ii)'}", "F": anova2[(v, u, t)].loc[term, "F"],
"df": f"(1, {anova2[(v, u, t)].loc[term, 'df_resid']:.0f})", "p": anova2[(v, u, t)].loc[term, "p"]}
for u in ("unit means", "variant: clones") for t in (1, 2)])
--- R: unit means n = 22, clones n = 128
unit means, (i) sequential
df ss F df_resid p
C(temp) 1.0 0.01840 6.45863 15.0 0.02258
C(phage) 1.0 0.79987 280.75216 15.0 0.00000
C(ctype) 1.0 1.26840 445.20716 15.0 0.00000
C(temp):C(phage) 1.0 0.00598 2.09891 15.0 0.16799
C(temp):C(ctype) 1.0 0.01856 6.51589 15.0 0.02208
C(phage):C(ctype) 1.0 0.24519 86.06293 15.0 0.00000
Residual 15.0 0.04274 NaN 15.0 NaN
unit means, (ii) type II
df ss F df_resid p
C(temp) 1.0 0.02547 8.93854 15.0 0.00916
C(phage) 1.0 0.08688 30.49509 15.0 0.00006
C(ctype) 1.0 1.26163 442.82892 15.0 0.00000
C(temp):C(phage) 1.0 0.00143 0.50203 15.0 0.48948
C(temp):C(ctype) 1.0 0.00111 0.39014 15.0 0.54161
C(phage):C(ctype) 1.0 0.24519 86.06293 15.0 0.00000
Residual 15.0 0.04274 NaN 15.0 NaN
variant: clones, (i) sequential
df ss F df_resid p
C(temp) 1.0 0.01888 3.32179 121.0 0.07084
C(phage) 1.0 5.53346 973.57673 121.0 0.00000
C(ctype) 1.0 6.54682 1151.87172 121.0 0.00000
C(temp):C(phage) 1.0 0.00150 0.26336 121.0 0.60876
C(temp):C(ctype) 1.0 0.16898 29.73127 121.0 0.00000
C(phage):C(ctype) 1.0 1.04599 184.03573 121.0 0.00000
Residual 121.0 0.68772 NaN 121.0 NaN
variant: clones, (ii) type II
df ss F df_resid p
C(temp) 1.0 0.12787 22.49711 121.0 0.00001
C(phage) 1.0 0.29417 51.75653 121.0 0.00000
C(ctype) 1.0 6.41029 1127.84927 121.0 0.00000
C(temp):C(phage) 1.0 0.00803 1.41250 121.0 0.23697
C(temp):C(ctype) 1.0 0.00017 0.02917 121.0 0.86467
C(phage):C(ctype) 1.0 1.04599 184.03573 121.0 0.00000
Residual 121.0 0.68772 NaN 121.0 NaN
--- mu: unit means n = 22, clones n = 128
unit means, (i) sequential
df ss F df_resid p
C(temp) 1.0 0.00020 26.62071 15.0 0.00012
C(phage) 1.0 0.00002 3.22601 15.0 0.09264
C(ctype) 1.0 0.00022 29.39329 15.0 0.00007
C(temp):C(phage) 1.0 0.00003 4.36031 15.0 0.05425
C(temp):C(ctype) 1.0 0.00000 0.00180 15.0 0.96676
C(phage):C(ctype) 1.0 0.00002 2.89816 15.0 0.10931
Residual 15.0 0.00011 NaN 15.0 NaN
unit means, (ii) type II
df ss F df_resid p
C(temp) 1.0 0.00017 23.63557 15.0 0.00021
C(phage) 1.0 0.00000 0.25502 15.0 0.62090
C(ctype) 1.0 0.00024 32.43611 15.0 0.00004
C(temp):C(phage) 1.0 0.00001 1.14766 15.0 0.30097
C(temp):C(ctype) 1.0 0.00000 0.35674 15.0 0.55923
C(phage):C(ctype) 1.0 0.00002 2.89816 15.0 0.10931
Residual 15.0 0.00011 NaN 15.0 NaN
variant: clones, (i) sequential
df ss F df_resid p
C(temp) 1.0 0.00155 59.31362 121.0 0.00000
C(phage) 1.0 0.00019 7.23643 121.0 0.00815
C(ctype) 1.0 0.00100 38.27831 121.0 0.00000
C(temp):C(phage) 1.0 0.00013 5.06884 121.0 0.02616
C(temp):C(ctype) 1.0 0.00000 0.00751 121.0 0.93110
C(phage):C(ctype) 1.0 0.00009 3.46629 121.0 0.06506
Residual 121.0 0.00317 NaN 121.0 NaN
variant: clones, (ii) type II
df ss F df_resid p
C(temp) 1.0 0.00133 51.01023 121.0 0.00000
C(phage) 1.0 0.00002 0.70867 121.0 0.40155
C(ctype) 1.0 0.00109 41.83019 121.0 0.00000
C(temp):C(phage) 1.0 0.00003 1.17329 121.0 0.28088
C(temp):C(ctype) 1.0 0.00002 0.63119 121.0 0.42847
C(phage):C(ctype) 1.0 0.00009 3.46629 121.0 0.06506
Residual 121.0 0.00317 NaN 121.0 NaN
--- A: unit means n = 22, clones n = 128
unit means, (i) sequential
df ss F df_resid p
C(temp) 1.0 0.01004 4.66642 15.0 0.04735
C(phage) 1.0 0.04215 19.59854 15.0 0.00049
C(ctype) 1.0 0.18917 87.96521 15.0 0.00000
C(temp):C(phage) 1.0 0.00703 3.26783 15.0 0.09074
C(temp):C(ctype) 1.0 0.00693 3.22191 15.0 0.09283
C(phage):C(ctype) 1.0 0.00624 2.89994 15.0 0.10921
Residual 15.0 0.03226 NaN 15.0 NaN
unit means, (ii) type II
df ss F df_resid p
C(temp) 1.0 0.01080 5.02008 15.0 0.04062
C(phage) 1.0 0.00071 0.33064 15.0 0.57381
C(ctype) 1.0 0.17196 79.96093 15.0 0.00000
C(temp):C(phage) 1.0 0.00273 1.26840 15.0 0.27776
C(temp):C(ctype) 1.0 0.00279 1.29657 15.0 0.27270
C(phage):C(ctype) 1.0 0.00624 2.89994 15.0 0.10921
Residual 15.0 0.03226 NaN 15.0 NaN
variant: clones, (i) sequential
df ss F df_resid p
C(temp) 1.0 0.06669 13.60360 121.0 0.00034
C(phage) 1.0 0.37440 76.37280 121.0 0.00000
C(ctype) 1.0 0.99755 203.48882 121.0 0.00000
C(temp):C(phage) 1.0 0.05649 11.52295 121.0 0.00093
C(temp):C(ctype) 1.0 0.04595 9.37361 121.0 0.00271
C(phage):C(ctype) 1.0 0.02727 5.56300 121.0 0.01995
Residual 121.0 0.59317 NaN 121.0 NaN
variant: clones, (ii) type II
df ss F df_resid p
C(temp) 1.0 0.04486 9.15165 121.0 0.00304
C(phage) 1.0 0.00533 1.08739 121.0 0.29913
C(ctype) 1.0 0.89838 183.25894 121.0 0.00000
C(temp):C(phage) 1.0 0.01553 3.16750 121.0 0.07763
C(temp):C(ctype) 1.0 0.01811 3.69450 121.0 0.05694
C(phage):C(ctype) 1.0 0.02727 5.56300 121.0 0.01995
Residual 121.0 0.59317 NaN 121.0 NaN
Table S3a (R): temperature F 6.459 p .023; phage 280.752; colony 445.207; TxP 2.099 p .168; TxC 6.516; PxC 86.063; residual SS 0.043 (15)
Table S3b (mu): temperature 26.621; phage 3.226; colony 29.393; TxP 4.360; TxC 0.002 p .967; PxC 2.898 p .109; residual 1.109e-4 (15)
Table S3c (A): temperature 4.666; phage 19.598; colony 87.965; TxP 3.268; TxC 3.222 p .093; PxC 2.900 p .109; residual 0.032 (15)
fig2a-phage-exposed-more-resistant¶
Clones previously exposed to phage were more phage-resistant than clones from the no-phage treatments, at both temperatures. (phage F(1, 15) = 280.752)
Criterion: on unit means, reproduced when the phage effect has p < 0.05 under both sums of squares readings and the mean R of +P clones exceeds that of -P clones at each temperature.
print(term_table("R", "C(phage)").to_string(index=False))
mR = c128.groupby(["temp", "phage"]).R.agg(["mean", "size"])
print("\nMean R of clones by temperature and phage (ancestor", round(ancestor["R"], 4), "):")
print(mR.round(4).to_string())
print("Mean R by treatment and colony type:")
print(c128.groupby(["treatment", "ctype"]).R.agg(["mean", "min", "max", "size"]).round(4).to_string())
fig, axs = plt.subplots(1, 3, figsize=(11, 3.2))
for ax, v in zip(axs, clone_sheets):
groups = [(f"{t}\n{c}", g[v].values) for (t, c), g in c128.groupby(["treatment", "ctype"])]
ax.boxplot([x for _, x in groups], tick_labels=[l for l, _ in groups])
ax.axhline(ancestor[v], ls="--", color="grey")
ax.set_title({"R": "2A resistance", "mu": "2B max growth rate", "A": "2C carrying capacity"}[v]); ax.tick_params(labelsize=7)
plt.tight_layout(); plt.show()
print("Plotted ranges:", {v: (round(c128[v].min(), 4), round(c128[v].max(), 4)) for v in clone_sheets}, "; all points shown")
reading F df p
unit means, (i) 280.752161 (1, 15) 4.030262e-11
unit means, (ii) 30.495092 (1, 15) 5.856379e-05
variant: clones, (i) 973.576732 (1, 121) 1.044166e-59
variant: clones, (ii) 51.756532 (1, 121) 5.717589e-11
Mean R of clones by temperature and phage (ancestor 0.3208 ):
mean size
temp phage
25 0 0.5980 32
1 0.9481 32
35 0 0.5566 32
1 1.0381 32
Mean R by treatment and colony type:
mean min max size
treatment ctype
25_0 m 0.2459 0.1833 0.5311 16
scv 0.9500 0.8835 0.9976 16
25_1 m 0.9055 0.8637 0.9293 4
scv 0.9542 0.9118 0.9895 28
35_0 m 0.2950 0.2112 0.9326 20
scv 0.9926 0.9472 1.0375 12
35_1 scv 1.0381 0.9638 1.1187 32
Plotted ranges: {'R': (np.float64(0.1833), np.float64(1.1187)), 'mu': (np.float64(0.0367), np.float64(0.0717)), 'A': (np.float64(0.6097), np.float64(1.5112))} ; all points shown
Result: reproduced. On unit means the phage effect is F(1, 15) = 280.75 with sequential sums of squares, which reproduces Table S3a term by term, and F = 30.50 with type II, both p < 1e-4; the clone variant gives the same (F = 973.58 and 51.76). Mean R of +P clones exceeds that of -P clones at 25 degrees C (0.9481 against 0.5980) and at 35 degrees C (1.0381 against 0.5566).
fig2a-temperature-resistance¶
Clones from the 35 degrees C treatments had higher phage resistance than clones from 25 degrees C. (temperature F(1, 15) = 6.459, p < .05)
Criterion: the ANOVA and readings above; reproduced when the temperature effect has p < 0.05 with mean R higher at 35 degrees C under both sums of squares readings on unit means.
print(term_table("R", "C(temp)").to_string(index=False))
print("\nMean R, clones: 25 C", round(c128[c128.temp == 25].R.mean(), 4), "; 35 C", round(c128[c128.temp == 35].R.mean(), 4))
u = pt.unit_means(c128, "R")
print("Mean R, unit means: 25 C", round(u[u.temp == 25].R.mean(), 4), "; 35 C", round(u[u.temp == 35].R.mean(), 4))
print("Within colony type, clones:", c128.groupby(["ctype", "temp"]).R.mean().round(4).to_dict())
reading F df p
unit means, (i) 6.458627 (1, 15) 0.022580
unit means, (ii) 8.938543 (1, 15) 0.009162
variant: clones, (i) 3.321790 (1, 121) 0.070837
variant: clones, (ii) 22.497115 (1, 121) 0.000006
Mean R, clones: 25 C 0.773 ; 35 C 0.7973
Mean R, unit means: 25 C 0.7314 ; 35 C 0.7902
Within colony type, clones: {('m', 25): 0.3778, ('m', 35): 0.295, ('scv', 25): 0.9527, ('scv', 35): 1.0257}
Result: reproduced. On unit means the temperature effect has p = 0.0226 with sequential sums of squares (F(1, 15) = 6.459, as Table S3a) and p = 0.0092 with type II, and mean R is higher at 35 degrees C (0.7902 against 0.7314). In the clone variant mean R is also higher at 35 degrees C (0.7973 against 0.7730), with p = 6e-6 (type II) and p = 0.0708 (sequential, F(1, 121) = 3.32); the variant decides nothing. Within colony type the direction differs: mucoid clones have higher R at 25 degrees C (0.378 against 0.295), and all 35 degrees C +P clones are SCVs, so the temperature effect on all clones comes in part from the colony types each treatment yielded.
fig2a-no-temperature-phage-interaction¶
The increase in phage resistance after phage exposure did not differ between the two temperatures. (phage x temperature F(1, 15) = 2.099, p = .168)
Criterion: a claim of absence judged on the phage by temperature term of the ANOVA above at p < .05; reproduced when p >= 0.05. The test is on unit means; the clone unit and the sums of squares readings are shown as variants.
print(term_table("R", "C(temp):C(phage)").to_string(index=False))
inc = c128.groupby(["temp", "phage"]).R.mean().unstack()
print("\nIncrease in mean R with phage (clones): 25 C", round(inc.loc[25, 1] - inc.loc[25, 0], 4),
"; 35 C", round(inc.loc[35, 1] - inc.loc[35, 0], 4))
reading F df p
unit means, (i) 2.098915 (1, 15) 0.167986
unit means, (ii) 0.502029 (1, 15) 0.489476
variant: clones, (i) 0.263360 (1, 121) 0.608757
variant: clones, (ii) 1.412504 (1, 121) 0.236969
Increase in mean R with phage (clones): 25 C 0.3502 ; 35 C 0.4815
Result: reproduced. On unit means the phage by temperature term has p = 0.168 with sequential sums of squares (the paper's F of 2.099) and p = 0.489 with type II; the clone variant gives 0.609 and 0.237. The increase in mean R with phage is 0.3502 at 25 degrees C and 0.4815 at 35 degrees C.
fig2a-scv-more-resistant¶
SCVs were much more phage-resistant than mucoid clones, whether or not they had been exposed to phage during selection, with a significant colony type by phage interaction. (colony type F(1, 15) = 445.207; colony type x phage F(1, 15) = 86.063)
Criterion: reproduced when the colony type effect and the colony type by phage interaction both have p < 0.05 under both sums of squares readings on unit means and the mean R of SCVs exceeds that of mucoid clones both among -P and among +P clones.
print("colony type:"); print(term_table("R", "C(ctype)").to_string(index=False))
print("colony type x phage:"); print(term_table("R", "C(phage):C(ctype)").to_string(index=False))
print("\nMean R of clones by phage and colony type:")
print(c128.groupby(["phage", "ctype"]).R.agg(["mean", "size"]).round(4).to_string())
colony type:
reading F df p
unit means, (i) 445.207162 (1, 15) 1.450555e-12
unit means, (ii) 442.828924 (1, 15) 1.508147e-12
variant: clones, (i) 1151.871719 (1, 121) 1.122056e-63
variant: clones, (ii) 1127.849273 (1, 121) 3.556804e-63
colony type x phage:
reading F df p
unit means, (i) 86.062935 (1, 15) 1.328734e-07
unit means, (ii) 86.062935 (1, 15) 1.328734e-07
variant: clones, (i) 184.035734 (1, 121) 4.703043e-26
variant: clones, (ii) 184.035734 (1, 121) 4.703043e-26
Mean R of clones by phage and colony type:
mean size
phage ctype
0 m 0.2731 36
scv 0.9683 28
1 m 0.9055 4
scv 0.9989 60
Result: reproduced. On unit means colony type (F(1, 15) = 445.21, sequential) and colony type by phage (F = 86.06) have p < 1e-6 under both sums of squares readings, and the clone variant gives the same. Mean R of SCVs exceeds that of mucoid clones among -P clones (0.9683 against 0.2731) and among +P clones (0.9989 against 0.9055).
fig2a-mucoid-25p-replicate-resistant¶
Mucoid clones from one 25 degrees C phage replicate were the exception, showing increased phage resistance. (phage F(1, 18) = 174.5)
Criterion: one-way ANOVA of clone R; readings (a) against the mucoid clones of 25 degrees C -P, (b) against the mucoid clones of both -P treatments. Reproduced when R of these clones is higher and p < 0.05 under each reading.
exc_m = c128[(c128.treatment == "25_1") & (c128.ctype == "m")]
print("25 C +P mucoid clones:", exc_m.strain.tolist(), "tube", exc_m.name_rep.unique().tolist(), "R", exc_m.R.round(4).tolist())
for label, ref in (("(a) 25 C -P mucoid", c128[(c128.treatment == "25_0") & (c128.ctype == "m")]),
("(b) all -P mucoid", c128[(c128.phage == 0) & (c128.ctype == "m")])):
f, p, d1, d2 = pt.one_way(exc_m.R, ref.R)
print(f"{label}: n = {len(ref)}, mean R {ref.R.mean():.4f} against {exc_m.R.mean():.4f}; F({d1}, {d2}) = {f:.2f}, p = {p:.2e}")
25 C +P mucoid clones: ['6d2a1', '6d2a2', '6d2a3', '6d2a4'] tube [2] R [0.9235, 0.9293, 0.8637, 0.9054] (a) 25 C -P mucoid: n = 16, mean R 0.2459 against 0.9055; F(1, 18) = 174.50, p = 1.06e-10 (b) all -P mucoid: n = 36, mean R 0.2731 against 0.9055; F(1, 38) = 88.09, p = 1.93e-11
Result: reproduced. The four mucoid clones of 25 degrees C +P (all from tube 2) have mean R 0.9055. Against the 16 mucoid clones of 25 degrees C -P (mean 0.2459), F(1, 18) = 174.50, p = 1.1e-10, the paper's F and df; against the 36 mucoid clones of both -P treatments (0.2731), F(1, 38) = 88.09, p = 1.9e-11.
fig2b-scv-higher-growth-rate¶
SCVs had a higher maximum growth rate than mucoid clones, regardless of prior phage or temperature exposure. (colony type F(1, 15) = 29.393; phage x colony type p = .109; temperature x colony type p = .967)
Criterion: mu from a logistic grofit fit, with the code's best-of-four fit as a variant, then the ANOVA readings of Figure 2A (unit means, sequential and type II sums of squares; clones as a variant). Reproduced when the colony type effect has p < 0.05 with SCVs higher and both interactions have p >= 0.05 under both readings. The dataset holds one maximum growth rate per clone and not the growth curves, so the fits cannot be repeated and the logistic and best-of-four readings cannot be told apart; the deposited values are tested.
for term in ("C(ctype)", "C(phage):C(ctype)", "C(temp):C(ctype)"):
print(term); print(term_table("mu", term).to_string(index=False))
print("\nMean mu by colony type (clones):", c128.groupby("ctype").mu.mean().round(5).to_dict(), "; ancestor", round(ancestor["mu"], 5))
print("by treatment and colony type:"); print(c128.groupby(["treatment", "ctype"]).mu.mean().round(5).unstack().to_string())
C(ctype)
reading F df p
unit means, (i) 29.393292 (1, 15) 7.075080e-05
unit means, (ii) 32.436109 (1, 15) 4.245127e-05
variant: clones, (i) 38.278311 (1, 121) 8.644869e-09
variant: clones, (ii) 41.830185 (1, 121) 2.206192e-09
C(phage):C(ctype)
reading F df p
unit means, (i) 2.898159 (1, 15) 0.109307
unit means, (ii) 2.898159 (1, 15) 0.109307
variant: clones, (i) 3.466291 (1, 121) 0.065058
variant: clones, (ii) 3.466291 (1, 121) 0.065058
C(temp):C(ctype)
reading F df p
unit means, (i) 0.001796 (1, 15) 0.966759
unit means, (ii) 0.356744 (1, 15) 0.559229
variant: clones, (i) 0.007507 (1, 121) 0.931100
variant: clones, (ii) 0.631190 (1, 121) 0.428475
Mean mu by colony type (clones): {'m': 0.04585, 'scv': 0.05235} ; ancestor 0.05078
by treatment and colony type:
ctype m scv
treatment
25_0 0.04726 0.05679
25_1 0.05258 0.05600
35_0 0.04338 0.05084
35_1 NaN 0.04750
Result: reproduced, on the growth rates the dataset holds. On unit means colony type has p = 7e-5 with sequential sums of squares (F(1, 15) = 29.39, which reproduces Table S3b) and p = 4e-5 with type II, with SCVs higher (0.05235 against 0.04585 per clone). Phage by colony type has p = 0.109 under both readings and temperature by colony type p = 0.967 (sequential) and 0.559 (type II), so neither interaction is significant. The clone variant agrees (phage by colony type p = 0.065, temperature by colony type 0.93 and 0.43). The fits could not be repeated, as the dataset holds no growth curves.
fig2c-scv-lower-carrying-capacity¶
SCVs had a lower carrying capacity than mucoid clones, not affected by prior temperature or phage exposure. (colony type F(1, 15) = 87.965; temperature x colony type p = .093; phage x colony type p = .109)
Criterion: carrying capacity A from the fits of fig2b, same readings (unit means, both sums of squares; clones as a variant). Reproduced when the colony type effect has p < 0.05 with SCVs lower and both interactions have p >= 0.05 under both readings. As for fig2b, the deposited values are tested.
Result: reproduced, on the carrying capacities the dataset holds. On unit means colony type has p < 1e-6 under both sums of squares readings (F(1, 15) = 87.97 sequential, as Table S3c), with SCVs lower (1.1335 against 1.3568 per clone). Temperature by colony type has p = 0.093 (sequential) and 0.273 (type II), and phage by colony type p = 0.109 under both, so neither interaction is significant. In the clone variant phage by colony type has p = 0.0199 under both sums of squares and temperature by colony type p = 0.0027 (sequential) and 0.057 (type II); the claims file makes the clone unit a variant, so it decides nothing.
fig2d-scv-twitching-no-exposure-effect¶
Prior temperature or phage exposure had no effect on the twitching motility of SCVs. (Table S3d, all p > .05)
Criterion: a claim of absence judged on Table S3d, an ANOVA of twitching diameter over the 28 clones with temperature, phage, colony type and the three two-way interactions, sums of squares sequential and type II. Reproduced when the temperature, phage and temperature by phage terms all have p >= 0.05. The same ANOVA on SCVs alone is a variant. Colony types are those of the design reading of the data check; the reading by strain name is a variant.
twd = tw.assign(ctype=tw.type_design)
tw_anova = {typ: pt.clone_anova(twd, "diameter_mm", typ) for typ in (1, 2)}
for typ in (1, 2):
print("design reading,", "(i) sequential" if typ == 1 else "(ii) type II")
print(tw_anova[typ].round(4).to_string())
print("Table S3d: temperature SS 3.857 F 2.264 p .147; phage 0.210 F 0.123 p .729; colony 14.510 F 8.516; TxP 0.604 F 0.355 p .558;"
" TxC 0.081 F 0.047; PxC 27.040 F 15.870; residual SS 35.780 (21)")
ss_total = ((twd.diameter_mm - twd.diameter_mm.mean()) ** 2).sum()
s3d_total = 3.857 + 0.210 + 14.510 + 0.604 + 0.081 + 27.040 + 35.780
print(f"\nTotal sum of squares of the deposited diameters: {ss_total:.3f}; summed from Table S3d: {s3d_total:.3f}")
s3d_res_ms = 35.780 / 21
print("Table S3d F ratios recomputed as SS / (residual SS / 21): temperature", round(3.857 / s3d_res_ms, 3), "; colony",
round(14.510 / s3d_res_ms, 3), "; PxC", round(27.040 / s3d_res_ms, 3), "(each term has 1 df)")
print("If Table S3d's sums of squares are sequential, they add up to the total; the total does not depend on the colony types,"
" so then no assignment of colony types to the rows reproduces Table S3d.")
print("\nVariant, colony type by strain name (24 distinct names, 28 rows):")
print(pt.clone_anova(tw.assign(ctype=tw.type_as_named), "diameter_mm", 1).round(4).to_string())
scv = twd[twd.ctype == "scv"]
print("\nVariant, SCVs alone (n =", len(scv), "):")
print(sm.stats.anova_lm(ols("diameter_mm ~ C(temp) * C(phage)", scv).fit(), typ=1).round(4).to_string())
print("\nMean SCV diameter (mm) by treatment:", scv.groupby("treatment").diameter_mm.mean().round(3).to_dict())
design reading, (i) sequential
df ss F df_resid p
C(temp) 1.0 2.4800 1.8388 21.0 0.1895
C(phage) 1.0 0.8441 0.6259 21.0 0.4377
C(ctype) 1.0 13.4560 9.9771 21.0 0.0047
C(temp):C(phage) 1.0 0.0535 0.0397 21.0 0.8440
C(temp):C(ctype) 1.0 0.0806 0.0597 21.0 0.8093
C(phage):C(ctype) 1.0 27.0400 20.0492 21.0 0.0002
Residual 21.0 28.3224 NaN 21.0 NaN
design reading, (ii) type II
df ss F df_resid p
C(temp) 1.0 0.2467 0.1829 21.0 0.6732
C(phage) 1.0 2.3542 1.7455 21.0 0.2007
C(ctype) 1.0 12.2075 9.0514 21.0 0.0067
C(temp):C(phage) 1.0 5.4250 4.0225 21.0 0.0580
C(temp):C(ctype) 1.0 5.5421 4.1093 21.0 0.0555
C(phage):C(ctype) 1.0 27.0400 20.0492 21.0 0.0002
Residual 21.0 28.3224 NaN 21.0 NaN
Table S3d: temperature SS 3.857 F 2.264 p .147; phage 0.210 F 0.123 p .729; colony 14.510 F 8.516; TxP 0.604 F 0.355 p .558; TxC 0.081 F 0.047; PxC 27.040 F 15.870; residual SS 35.780 (21)
Total sum of squares of the deposited diameters: 72.277; summed from Table S3d: 82.082
Table S3d F ratios recomputed as SS / (residual SS / 21): temperature 2.264 ; colony 8.516 ; PxC 15.87 (each term has 1 df)
If Table S3d's sums of squares are sequential, they add up to the total; the total does not depend on the colony types, so then no assignment of colony types to the rows reproduces Table S3d.
Variant, colony type by strain name (24 distinct names, 28 rows):
df ss F df_resid p
C(temp) 1.0 2.4800 1.7769 22.0 0.1962
C(phage) 1.0 0.8441 0.6048 22.0 0.4450
C(ctype) 1.0 28.8459 20.6674 22.0 0.0002
C(temp):C(phage) 1.0 1.3020 0.9329 22.0 0.3446
C(temp):C(ctype) 1.0 8.0988 5.8026 22.0 0.0248
C(phage):C(ctype) 1.0 2.8170 2.0183 22.0 0.1694
Residual 22.0 30.7058 NaN 22.0 NaN
Variant, SCVs alone (n = 16 ):
df sum_sq mean_sq F PR(>F)
C(temp) 1.0 0.1502 0.1502 0.0927 0.7659
C(phage) 1.0 2.0664 2.0664 1.2762 0.2807
C(temp):C(phage) 1.0 5.4250 5.4250 3.3504 0.0921
Residual 12.0 19.4305 1.6192 NaN NaN
Mean SCV diameter (mm) by treatment: {'25_0': 8.467, '25_1': 10.35, '35_0': 9.825, '35_1': 9.379}
/home/runner/work/encore/encore/entries/dryad-tdz08kq5d-wang-2024/phage_temp.py:211: SingularMatrixWarning: The design matrix is rank-deficient. The model parameters are not uniquely determined.
a = sm.stats.anova_lm(ols(f"{dv} ~ {formula}", d).fit(), typ=typ)
Result: reproduced. Under the design reading the temperature, phage and temperature by phage terms have p = 0.190, 0.438 and 0.844 (sequential) and 0.673, 0.201 and 0.058 (type II); by strain name 0.196, 0.445 and 0.345 (a rank-deficient fit, as that reading leaves 25 degrees C -P without SCVs); and on SCVs alone 0.766, 0.281 and 0.092. None is below 0.05. The deposited twitching values do not reproduce Table S3d, under any assignment of colony types, if the sums of squares of Table S3d are sequential so that they add up to the total: the total sum of squares of the data is 72.28 against 82.08 summed from Table S3d, and a total does not depend on how the rows are labelled. Table S3d does not say which sums of squares it gives. Sequential ones are likely, since Tables S3a to S3c, S3e and S3f match with sequential sums of squares and the F ratios of S3d are its sums of squares over its residual mean square (2.264, 8.516 and 15.870 recomputed above), but this is an inference. Only the colony type by phage sum of squares (27.04) and the temperature by colony type one (0.081) agree, and the residual sum of squares is 28.32 against 35.78. The design reading of the colony types (the second 25_0 block as the SCVs) is this attempt's inference from the design, not the claims file's; the outcome is the same under the reading by strain name and on SCVs alone.
fig2d-mucoid-25p-reduced-twitching¶
Phage-resistant mucoid clones from the 25 degrees C phage treatment showed clearly reduced twitching motility. (F(1, 6) = 27.22, p < .01)
Criterion: one-way ANOVA of clone diameter, four against four; readings (a) against the mucoid clones of 25 degrees C -P, (b) against those of 35 degrees C -P. Reproduced when the 25 degrees C +P mucoid clones have the smaller mean and p < 0.05 under both. The ancestor's diameter is given alongside.
m25p = twd[(twd.treatment == "25_1") & (twd.ctype == "m")]
print("25 C +P mucoid:", m25p.strain.tolist(), m25p.diameter_mm.round(3).tolist(), "mean", round(m25p.diameter_mm.mean(), 3))
for label, tr in (("(a) 25 C -P mucoid", "25_0"), ("(b) 35 C -P mucoid", "35_0")):
ref = twd[(twd.treatment == tr) & (twd.ctype == "m")]
f, p, d1, d2 = pt.one_way(m25p.diameter_mm, ref.diameter_mm)
print(f"{label}: {ref.diameter_mm.round(3).tolist()} mean {ref.diameter_mm.mean():.3f}; F({d1}, {d2}) = {f:.2f}, p = {p:.4f}")
print("Ancestor:", round(anc_tw, 3), "mm (the data's unit)")
ref_named = tw[(tw.treatment == "25_0")]
f, p, d1, d2 = pt.one_way(m25p.diameter_mm, ref_named.diameter_mm)
print(f"Variant, all eight 25_0 rows as mucoid (strain-name reading): F({d1}, {d2}) = {f:.2f}, p = {p:.4f}")
25 C +P mucoid: ['6d2a1', '6d2a2', '6d2a3', '6d2a4'] [5.783, 6.633, 6.3, 6.25] mean 6.242 (a) 25 C -P mucoid: [11.333, 8.75, 8.75, 9.4] mean 9.558; F(1, 6) = 27.22, p = 0.0020 (b) 35 C -P mucoid: [7.033, 9.167, 9.7, 8.35] mean 8.562; F(1, 6) = 14.66, p = 0.0087 Ancestor: 9.383 mm (the data's unit) Variant, all eight 25_0 rows as mucoid (strain-name reading): F(1, 10) = 20.62, p = 0.0011
Result: reproduced. The four 25 degrees C +P mucoid clones have a mean diameter of 6.242 mm. Against the 25 degrees C -P mucoid clones (9.558 mm) F(1, 6) = 27.22, p = 0.0020, the paper's F; against the 35 degrees C -P mucoid clones (8.562 mm) F(1, 6) = 14.66, p = 0.0087. The ancestor's diameter is 9.383 mm.
fig2e-mucoid-biofilm-increased¶
Mucoid clones formed more biofilm than the ancestor, while SCVs kept ancestral-level biofilm formation. (colony type F(1, 21) = 42.370)
Criterion: the ANOVA of Table S3e; readings of Figure 2E: (a) the ancestor lies below the interquartile range of the mucoid clones and within that of the SCVs, (b) it lies below the median of the mucoid clones and within the range of the SCVs. Reproduced when the colony type effect has p < 0.05 with mucoid higher and the comparison holds under both readings; blocked when the readings differ.
bf_anova = {typ: pt.clone_anova(bf, "biofilm", typ) for typ in (1, 2)}
for typ in (1, 2):
print("(i) sequential" if typ == 1 else "(ii) type II"); print(bf_anova[typ].round(4).to_string())
print("Table S3e: temperature SS 0.346 F 1.806 p .193; phage 0.000 F 0.000 p .993; colony 8.555 F 42.370; TxP 0.053 F 0.264 p .613;"
" TxC 0.020; PxC 0.490 F 2.426; residual 4.240 (21)")
q = bf.groupby("ctype").biofilm.quantile([0, 0.25, 0.5, 0.75, 1]).unstack()
print("\nBiofilm OD570 quantiles by colony type (ancestor", round(anc_bf, 4), "):"); print(q.round(4).to_string())
m_q, s_q = q.loc["m"], q.loc["scv"]
print("(a) ancestor below mucoid IQR:", anc_bf < m_q[0.25], "; within SCV IQR:", s_q[0.25] <= anc_bf <= s_q[0.75])
print("(b) ancestor below mucoid median:", anc_bf < m_q[0.5], "; within SCV range:", s_q[0] <= anc_bf <= s_q[1])
print("SCV clones above the ancestor:", int((bf[bf.ctype == "scv"].biofilm > anc_bf).sum()), "of", int((bf.ctype == "scv").sum()),
"; SCV mean", round(bf[bf.ctype == "scv"].biofilm.mean(), 4), f"({100 * (bf[bf.ctype == 'scv'].biofilm.mean() / anc_bf - 1):.0f}% above the ancestor)")
print("By treatment and colony type:"); print(bf.groupby(["treatment", "ctype"]).biofilm.agg(["mean", "min", "max"]).round(4).to_string())
(i) sequential
df ss F df_resid p
C(temp) 1.0 0.3646 1.8056 21.0 0.1934
C(phage) 1.0 0.0000 0.0001 21.0 0.9930
C(ctype) 1.0 8.5547 42.3701 21.0 0.0000
C(temp):C(phage) 1.0 0.0533 0.2640 21.0 0.6128
C(temp):C(ctype) 1.0 0.0200 0.0989 21.0 0.7563
C(phage):C(ctype) 1.0 0.4899 2.4263 21.0 0.1343
Residual 21.0 4.2400 NaN 21.0 NaN
(ii) type II
df ss F df_resid p
C(temp) 1.0 0.0214 0.1062 21.0 0.7477
C(phage) 1.0 0.3784 1.8743 21.0 0.1855
C(ctype) 1.0 7.6965 38.1194 21.0 0.0000
C(temp):C(phage) 1.0 0.0409 0.2026 21.0 0.6572
C(temp):C(ctype) 1.0 0.2231 1.1050 21.0 0.3051
C(phage):C(ctype) 1.0 0.4899 2.4263 21.0 0.1343
Residual 21.0 4.2400 NaN 21.0 NaN
Table S3e: temperature SS 0.346 F 1.806 p .193; phage 0.000 F 0.000 p .993; colony 8.555 F 42.370; TxP 0.053 F 0.264 p .613; TxC 0.020; PxC 0.490 F 2.426; residual 4.240 (21)
Biofilm OD570 quantiles by colony type (ancestor 0.4887 ):
0.00 0.25 0.50 0.75 1.00
ctype
m 0.6033 1.4914 1.8157 2.1395 2.9090
scv 0.5763 0.6333 0.7263 0.7758 0.9507
(a) ancestor below mucoid IQR: True ; within SCV IQR: False
(b) ancestor below mucoid median: True ; within SCV range: False
SCV clones above the ancestor: 16 of 16 ; SCV mean 0.7204 (47% above the ancestor)
By treatment and colony type:
mean min max
treatment ctype
25_0 m 1.5190 1.1003 1.9343
scv 0.7772 0.7350 0.8110
25_1 m 2.1693 1.6970 2.9090
scv 0.7276 0.7100 0.7477
35_0 m 1.8263 0.6033 2.8800
scv 0.6121 0.5763 0.6373
35_1 scv 0.7648 0.5870 0.9507
Result: partially reproduced. The ANOVA reproduces Table S3e (colony type F(1, 21) = 42.37, p < 1e-4, mucoid higher), and the ancestor (0.4887) lies below the interquartile range (from 1.4914) and the median (1.8157) of the mucoid clones. It lies below every SCV clone as well (16 of 16 higher; SCV range 0.5763 to 0.9507, mean 0.7204, 47% above the ancestor), so it is within neither the interquartile range nor the range of the SCVs under either reading. The mucoid half of the claim and the colony type effect hold; the SCV half fails. The criterion defines no partial outcome; the claim holds in part, which is what the protocol calls partially reproduced (changed from not reproduced after verification round 1).
fig2e-biofilm-no-exposure-effect¶
Prior temperature or phage exposure had no significant effect on biofilm formation. (Table S3e, all p > .05)
Criterion: a claim of absence judged on the ANOVA of Table S3e; reproduced when the temperature, phage and temperature by phage terms all have p >= 0.05.
for typ in (1, 2):
t = bf_anova[typ]
print("(i) sequential" if typ == 1 else "(ii) type II", {k: (round(t.loc[k, "F"], 3), round(t.loc[k, "p"], 4))
for k in ("C(temp)", "C(phage)", "C(temp):C(phage)")})
print("Mean biofilm by temperature:", bf.groupby("temp").biofilm.mean().round(4).to_dict(), "; by phage:", bf.groupby("phage").biofilm.mean().round(4).to_dict())
(i) sequential {'C(temp)': (np.float64(1.806), np.float64(0.1934)), 'C(phage)': (np.float64(0.0), np.float64(0.993)), 'C(temp):C(phage)': (np.float64(0.264), np.float64(0.6128))}
(ii) type II {'C(temp)': (np.float64(0.106), np.float64(0.7477)), 'C(phage)': (np.float64(1.874), np.float64(0.1855)), 'C(temp):C(phage)': (np.float64(0.203), np.float64(0.6572))}
Mean biofilm by temperature: {25: 1.2983, 35: 1.0677} ; by phage: {0: 1.1836, 1: 1.2206}
Result: reproduced. Temperature, phage and temperature by phage have p = 0.193, 0.993 and 0.613 (sequential, the F values of Table S3e) and 0.748, 0.186 and 0.657 (type II).
fig2f-phage-exposure-reduces-virulence¶
Prior exposure to phage reduced virulence in tomato, irrespective of the temperature at which the bacteria evolved; the main effect of temperature was not significant. (phage F(1, 21) = 20.310; temperature p = .173)
Criterion: disease incidence per clone (the dataset holds one value per clone, the share of wilted plants over its three blocks); the ANOVA terms of Table S3f, sums of squares read both ways. Reproduced when the phage effect has p < 0.05 with lower incidence in +P clones and the temperature effect has p >= 0.05.
dz_anova = {typ: pt.clone_anova(dz, "incidence", typ) for typ in (1, 2)}
for typ in (1, 2):
print("(i) sequential" if typ == 1 else "(ii) type II"); print(dz_anova[typ].round(4).to_string())
print("Table S3f: temperature F 1.987 p .173; phage F 20.310; colony F 33.413; TxP 0.910; TxC 7.714; PxC 9.992 p .004; residual 0.214 (21)")
print("\nMean incidence by phage:", dz.groupby("phage").incidence.mean().round(4).to_dict(), "; by temperature:", dz.groupby("temp").incidence.mean().round(4).to_dict())
print("By treatment and colony type:"); print(dz.groupby(["treatment", "ctype"]).incidence.agg(["mean", "min", "max", "size"]).round(4).to_string())
(i) sequential
df ss F df_resid p
C(temp) 1.0 0.0203 1.9872 21.0 0.1733
C(phage) 1.0 0.2074 20.3104 21.0 0.0002
C(ctype) 1.0 0.3413 33.4128 21.0 0.0000
C(temp):C(phage) 1.0 0.0093 0.9097 21.0 0.3511
C(temp):C(ctype) 1.0 0.0788 7.7136 21.0 0.0113
C(phage):C(ctype) 1.0 0.1020 9.9918 21.0 0.0047
Residual 21.0 0.2145 NaN 21.0 NaN
(ii) type II
df ss F df_resid p
C(temp) 1.0 0.0046 0.4535 21.0 0.5080
C(phage) 1.0 0.0680 6.6612 21.0 0.0174
C(ctype) 1.0 0.3478 34.0504 21.0 0.0000
C(temp):C(phage) 1.0 0.0000 0.0000 21.0 1.0000
C(temp):C(ctype) 1.0 0.0069 0.6802 21.0 0.4188
C(phage):C(ctype) 1.0 0.1020 9.9918 21.0 0.0047
Residual 21.0 0.2145 NaN 21.0 NaN
Table S3f: temperature F 1.987 p .173; phage F 20.310; colony F 33.413; TxP 0.910; TxC 7.714; PxC 9.992 p .004; residual 0.214 (21)
Mean incidence by phage: {0: 0.1806, 1: 0.0} ; by temperature: {25: 0.0799, 35: 0.1343}
By treatment and colony type:
mean min max size
treatment ctype
25_0 m 0.3194 0.1667 0.4444 4
scv 0.0000 0.0000 0.0000 4
25_1 m 0.0000 0.0000 0.0000 4
scv 0.0000 0.0000 0.0000 4
35_0 m 0.4028 0.0556 0.5556 4
scv 0.0000 0.0000 0.0000 4
35_1 scv 0.0000 0.0000 0.0000 4
Result: reproduced. The ANOVA reproduces Table S3f: phage F(1, 21) = 20.310, p = 0.0002 (type II p = 0.0174), with mean incidence 0 in all +P clones against 0.1806 in -P clones; temperature p = 0.173 (type II 0.508).
fig2f-scv-avirulent¶
All SCVs showed complete loss of virulence regardless of prior phage exposure. (colony type F(1, 21) = 33.413)
Criterion: reproduced when every SCV clone has a disease incidence of 0 in every block and the colony type effect has p < 0.05 with SCVs lower; partially reproduced when the effect holds but some SCV clones wilt plants. The dataset holds the pooled incidence per clone; a pooled value of 0 means no wilted plant in any block.
s = dz[dz.ctype == "scv"]
print("SCV clones:", len(s), "; with incidence 0:", int((s.incidence == 0).sum()), "; values:", sorted(s.incidence.unique().tolist()))
for typ in (1, 2):
t = dz_anova[typ]
print("(i) sequential" if typ == 1 else "(ii) type II", "colony type F =", round(t.loc["C(ctype)", "F"], 3), "p =", f"{t.loc['C(ctype)', 'p']:.2e}")
print("Mean incidence: mucoid", round(dz[dz.ctype == "m"].incidence.mean(), 4), "; SCV", round(s.incidence.mean(), 4))
SCV clones: 16 ; with incidence 0: 16 ; values: [0.0] (i) sequential colony type F = 33.413 p = 9.75e-06 (ii) type II colony type F = 34.05 p = 8.60e-06 Mean incidence: mucoid 0.2408 ; SCV 0.0
Result: reproduced. All 16 SCV clones have an incidence of 0, so no plant wilted in any of their blocks; colony type F(1, 21) = 33.41, p = 9.7e-06, with SCVs lower (0 against 0.2408).
fig2f-mucoid-25p-avirulent¶
The phage-resistant mucoid clones from the 25 degrees C phage treatment lost their virulence. (colony type x phage F(1, 21) = 9.992)
Criterion: reproduced when the colony type by phage interaction has p < 0.05 and the mean incidence of the four 25 degrees C +P mucoid clones is below that of every -P mucoid group and of the ancestor; whether each of the four has 0 is reported.
for typ in (1, 2):
t = dz_anova[typ]
print("(i) sequential" if typ == 1 else "(ii) type II", "colony type x phage F =", round(t.loc["C(phage):C(ctype)", "F"], 3),
"p =", round(t.loc["C(phage):C(ctype)", "p"], 4))
mz = dz[dz.ctype == "m"].groupby("treatment").incidence.agg(["mean", "size"])
print(mz.round(4).to_string()); print("ancestor", round(anc_dz, 4))
print("25 C +P mucoid clones:", dz[(dz.treatment == "25_1") & (dz.ctype == "m")][["strain", "incidence"]].values.tolist())
(i) sequential colony type x phage F = 9.992 p = 0.0047
(ii) type II colony type x phage F = 9.992 p = 0.0047
mean size
treatment
25_0 0.3194 4
25_1 0.0000 4
35_0 0.4028 4
ancestor 0.7222
25 C +P mucoid clones: [['6d2a1', 0.0], ['6d2a2', 0.0], ['6d2a3', 0.0], ['6d2a4', 0.0]]
Result: reproduced. Colony type by phage F(1, 21) = 9.992, p = 0.0047. The four 25 degrees C +P mucoid clones each have an incidence of 0, below the -P mucoid groups (0.3194 at 25 degrees C, 0.4028 at 35 degrees C) and the ancestor (0.7222).
sequencing-mutation-count¶
Sequenced clones carried on average 8.5 nonsynonymous mutations, and the number did not depend on colony morphology.
Criterion: variants called with Snippy from the reads; counted per clone are variants that change a protein (small deletions, insertions and SNPs); the mean within 8.5 +- 0.85 and a one-way ANOVA of mucoid against SCV counts with p >= 0.05. The reads are not public (data check) and the dataset's variant table is aggregated by group, so the mean per clone is the sum over protein-changing variants of the number of clones carrying each, divided by the 28 clones. Per-clone counts, which the ANOVA needs, cannot be formed; group means are shown.
pc = var[var.protein_changing]
print("Protein-changing variant rows:", len(pc), "; effects:", pc.EFFECT.str.split(" ").str[0].value_counts().to_dict())
n_clones = var.groupby("group").size().size * 4
total = pc.num_samples_with_snp.sum()
print(f"Sum over rows of carrying clones: {total}; clones: {n_clones}; mean per clone {total / n_clones:.3f} (8.5 +- 0.85)")
per_group = pc.groupby("group").num_samples_with_snp.sum() / 4
print("\nMean protein-changing variants per clone by group:"); print(per_group.round(3).to_string())
by_type = pc.groupby("colonial_morphology").num_samples_with_snp.sum() / var.groupby("colonial_morphology").group.nunique() / 4
print("By colony type (mucoid 12 clones, SCV 16):", by_type.round(3).to_dict())
print("\nVariant, leaving out protein-changing variants present in all seven groups:",
round(pc[pc.groupby(["CHROM", "POS"]).group.transform("nunique") < 7].num_samples_with_snp.sum() / n_clones, 3), "per clone")
print("Variant, every variant in the table (intergenic included):", round(var.num_samples_with_snp.sum() / n_clones, 3), "per clone")
Protein-changing variant rows: 81 ; effects: {'frameshift_variant': 46, 'conservative_inframe_insertion': 13, 'disruptive_inframe_insertion': 8, 'frameshift_variant&start_lost': 6, 'missense_variant': 4, 'disruptive_inframe_deletion': 4}
Sum over rows of carrying clones: 282; clones: 28; mean per clone 10.071 (8.5 +- 0.85)
Mean protein-changing variants per clone by group:
group
25,+P,mucoid 10.50
25,+P,scv 11.75
25,-P,mucoid 10.25
25,-P,scv 9.00
35,+P,scv 8.25
35,-P,mucoid 9.75
35,-P,scv 11.00
By colony type (mucoid 12 clones, SCV 16): {'mucoid': 10.167, 'scv': 10.0}
Variant, leaving out protein-changing variants present in all seven groups: 2.357 per clone
Variant, every variant in the table (intergenic included): 16.964 per clone
Result: not reproduced. The dataset's variant table gives 282 protein-changing variants over 28 clones, 10.07 per clone, outside 8.5 +- 0.85. Per clone, mucoid clones carry 10.17 and SCVs 10.00 on average; the ANOVA on per-clone counts cannot be run, as the table does not say which clone carries which variant.
table1-phcs-in-scvs¶
SCVs carried mutations in the quorum-sensing histidine kinase gene phcS (13 of 16 SCV clones), while no mucoid clone did.
Criterion: reproduced when the number of SCV clones with a phcS mutation is within 13 +- 1.3 and of mucoid clones 0. The table does not say whether two phcS variants of a group are in different clones, so the count is shown both ways (different clones: the sum; the same clones: the larger).
var["gene_name"] = var.LOCUS_TAG.map(names)
def carriers(gene):
g = var[var.gene_name == gene]
return g.groupby(["group", "colonial_morphology"]).num_samples_with_snp.agg(["sum", "max", "size"])
phcs = carriers("phcS")
print(var[var.gene_name == "phcS"][["group", "POS", "TYPE", "EFFECT", "num_samples_with_snp"]].to_string(index=False))
print(phcs.to_string())
scv_hi = phcs.xs("scv", level=1)["sum"].clip(upper=4).sum()
scv_lo = phcs.xs("scv", level=1)["max"].sum()
m_n = phcs.xs("mucoid", level=1)["sum"].sum() if "mucoid" in phcs.index.get_level_values(1) else 0
print(f"\nSCV clones with a phcS mutation: {scv_hi} if the variants are in different clones, {scv_lo} if they overlap (13 +- 1.3); mucoid clones: {m_n}")
group POS TYPE EFFECT num_samples_with_snp
35,-P,scv 3415554 snp missense_variant c.218G>A p.Gly73Asp 1
35,-P,scv 3415879 del disruptive_inframe_deletion c.552_560delCTTCCAGCT p.Phe185_Leu187del 3
35,+P,scv 3415879 del disruptive_inframe_deletion c.552_560delCTTCCAGCT p.Phe185_Leu187del 2
25,-P,scv 3415879 del disruptive_inframe_deletion c.552_560delCTTCCAGCT p.Phe185_Leu187del 3
25,+P,scv 3415879 del disruptive_inframe_deletion c.552_560delCTTCCAGCT p.Phe185_Leu187del 4
sum max size
group colonial_morphology
25,+P,scv scv 4 4 1
25,-P,scv scv 3 3 1
35,+P,scv scv 2 2 1
35,-P,scv scv 4 3 2
SCV clones with a phcS mutation: 13 if the variants are in different clones, 12 if they overlap (13 +- 1.3); mucoid clones: 0
Result: reproduced. phcS (RS-N_03191) variants are in SCV groups only: 4 of 4 at 25 degrees C +P, 3 at 25 degrees C -P, 2 at 35 degrees C +P and at 35 degrees C -P two variants carried by 1 and 3 clones. That is 13 SCV clones if the two 35 degrees C -P variants are in different clones and 12 if not, both within 13 +- 1.3, and 0 mucoid clones.
table1-pild-in-25p-mucoid¶
Two of the four mucoid phage-resistant clones from 25 degrees C with phage had mutations in pilD, which no other clone carried.
Criterion: reproduced when two of the four carry a pilD mutation and no other clone does.
print(var[var.gene_name == "pilD"][["group", "POS", "TYPE", "EFFECT", "num_samples_with_snp"]].to_string(index=False))
group POS TYPE EFFECT num_samples_with_snp 25,+P,mucoid 3523802 ins disruptive_inframe_insertion c.669_677dupGCTGGGCGC p.Ala226_Trp227insLeuGlyAla 2
Result: reproduced. One pilD (RS-N_03289) variant, an in-frame insertion, is carried by 2 of the 4 mucoid clones of 25 degrees C +P and by no clone of another group.
table1-t2ss-in-35p-scv¶
Three of four SCVs from 35 degrees C with phage had mutations in the type II secretion genes gspE or gspF.
Criterion: reproduced when three of the four carry a gspE or gspF mutation (Table 1: gspE 2 of 4, gspF 1 of 4) and no clone of another treatment does.
t2 = var[var.gene_name.isin(["gspE", "gspF"])]
print(t2[["group", "gene_name", "POS", "TYPE", "EFFECT", "num_samples_with_snp"]].to_string(index=False))
g35 = t2[t2.group == "35,+P,scv"]
print(f"35 C +P SCV clones with a gspE or gspF mutation: {min(4, g35.num_samples_with_snp.sum())} if in different clones, "
f"{g35.num_samples_with_snp.max()} if they overlap; other groups with one: {sorted(set(t2.group) - {'35,+P,scv'})}")
group gene_name POS TYPE EFFECT num_samples_with_snp 35,+P,scv gspF 405903 del frameshift_variant c.291delC p.Tyr98fs 1 35,+P,scv gspE 406950 snp missense_variant c.848G>A p.Gly283Asp 2 35 C +P SCV clones with a gspE or gspF mutation: 3 if in different clones, 2 if they overlap; other groups with one: []
Result: blocked, metadata insufficient. In the 35 degrees C +P SCV group, 2 clones carry a gspE missense variant and 1 a gspF frameshift; no other group carries either. The table does not say whether the gspF variant is in one of the gspE clones, so the number of clones with a T2SS mutation is 3 or 2, and the claim's three cannot be confirmed.
tables6-pil-mutants-resistant¶
The pilM and pilB mutants were fully resistant to phage NNP42, unlike their wild types.
Blocked: the dataset holds no measurements of the pilB and pilM mutants or of their wild types OE1-1 and QL-Rs1115 (no sheet or column names them; data check). The paper gives only the means of Table S6.
hits = {name: [c for c in sheet.columns if any(k in str(c).lower() for k in ("pilb", "pilm", "oe1", "ql-rs"))] +
[v for v in map(str, sheet.values.ravel()) if any(k in v.lower() for k in ("pilb", "pilm", "oe1-1", "ql-rs"))][:5]
for name, sheet in wb.items()}
print("Sheets with a pilB, pilM, OE1-1 or QL-Rs1115 column or value:", {k: v for k, v in hits.items() if v} or "none")
Sheets with a pilB, pilM, OE1-1 or QL-Rs1115 column or value: none
figs3-adsorption¶
pilD mutants showed reduced phage adsorption, while the phcS and the T2SS (gspE, gspF) mutants showed no reduction.
Criterion: one-way ANOVA across strains with Tukey comparisons; reproduced when every pilD mutant adsorbs less than the ancestor at p < 0.05 and no phcS, gspE or gspF mutant adsorbs significantly less.
f, p, d1, d2 = pt.one_way(*[g.rate for _, g in ads.groupby("group")])
print(f"One-way ANOVA across {ads.group.nunique()} strains: F({d1}, {d2}) = {f:.3f}, p = {p:.2e}")
tk = pairwise_tukeyhsd(ads.rate, ads.group)
tkt = pd.DataFrame(tk.summary().data[1:], columns=tk.summary().data[0])
vs_anc = tkt[(tkt.group1 == "anc") | (tkt.group2 == "anc")].copy()
vs_anc["mutant"] = np.where(vs_anc.group1 == "anc", vs_anc.group2, vs_anc.group1)
vs_anc["mutant_minus_anc"] = np.where(vs_anc.group1 == "anc", vs_anc.meandiff, -vs_anc.meandiff)
print(vs_anc[["mutant", "mutant_minus_anc", "p-adj", "reject"]].to_string(index=False))
print("\nMean adsorption by strain:", ads.groupby("group").rate.mean().round(4).to_dict())
One-way ANOVA across 13 strains: F(12, 26) = 19.746, p = 4.57e-10
mutant mutant_minus_anc p-adj reject
gspE_B+S4 -0.0183 1.0000 False
gspF_B+S1 -0.1100 0.7357 False
phcS_B+S2 0.0367 1.0000 False
phcS_B+S3 -0.0300 1.0000 False
phcS_B-S1 0.0133 1.0000 False
phcS_B-S2 0.0500 0.9992 False
phcS_D+S2 0.0650 0.9912 False
phcS_D+S3 0.0750 0.9730 False
phcS_D-S2 0.0633 0.9930 False
phcS_D-S4 0.0150 1.0000 False
pilD_D+a1 -0.4900 0.0000 True
pilD_D+a2 -0.3617 0.0000 True
Mean adsorption by strain: {'anc': 0.7733, 'gspE_B+S4': 0.755, 'gspF_B+S1': 0.6633, 'phcS_B+S2': 0.81, 'phcS_B+S3': 0.7433, 'phcS_B-S1': 0.7867, 'phcS_B-S2': 0.8233, 'phcS_D+S2': 0.8383, 'phcS_D+S3': 0.8483, 'phcS_D-S2': 0.8367, 'phcS_D-S4': 0.7883, 'pilD_D+a1': 0.2833, 'pilD_D+a2': 0.4117}
Result: reproduced. One-way ANOVA across the 13 strains F(12, 26) = 19.75, p = 4.6e-10. By Tukey's test both pilD mutants adsorb less than the ancestor (0.283 and 0.412 against 0.773, adjusted p < 0.001), and none of the eight phcS mutants, the gspE mutant or the gspF mutant differs from it (smallest adjusted p 0.736, gspF, 0.663).
rnaseq-differential-expression-count¶
After 6 h of phage exposure, 1477 R. solanacearum genes were differentially expressed, 707 up and 770 down.
Criterion: DESeq2 between five phage and five control samples; differentially expressed means |log2 fold change| > 1 and FDR < 0.05. The dataset holds the DESeq2 table; the pipeline from the reads (PRJNA1032726, about 9.1 GB) was not re-run. The thresholds are applied to the table.
de = full[(full.padj < 0.05) & (full.log2FoldChange.abs() > 1)]
print(f"Differentially expressed: {len(de)} (1477 +- 147.7); up {(de.log2FoldChange > 0).sum()} (707 +- 70.7); "
f"down {(de.log2FoldChange < 0).sum()} (770 +- 77)")
for fdr, lfc in ((0.05, 0.5), (0.01, 1), (0.05, 2)):
x = full[(full.padj < fdr) & (full.log2FoldChange.abs() > lfc)]
print(f" threshold variant FDR < {fdr}, |log2FC| > {lfc}: {len(x)}")
print("Genes without padj (left out by the thresholds):", int(full.padj.isna().sum()))
Differentially expressed: 1477 (1477 +- 147.7); up 707 (707 +- 70.7); down 770 (770 +- 77) threshold variant FDR < 0.05, |log2FC| > 0.5: 2766 threshold variant FDR < 0.01, |log2FC| > 1: 1435 threshold variant FDR < 0.05, |log2FC| > 2: 362 Genes without padj (left out by the thresholds): 32
Result: reproduced. In the dataset's DESeq2 table, 1477 genes have padj < 0.05 and |log2FC| > 1, 707 up and 770 down, the paper's counts exactly.
rnaseq-kegg-enrichment¶
Five KEGG pathways were enriched in the phage compared to the control treatment, namely propanoate metabolism, RNA degradation, synthesis and degradation of ketone bodies, bacterial secretion system and valine, leucine and isoleucine degradation.
Criterion: KOBAS enrichment of the differentially expressed genes at p < .05, readings uncorrected and corrected p. Reproduced when the five named pathways are enriched at p < 0.05 under both readings; blocked when the readings differ. The dataset's KEGG table (KOBAS output with Pvalue and FDR) is read; KOBAS was not re-run.
named = ["Propanoate metabolism", "RNA degradation", "Synthesis and degradation of ketone bodies",
"Bacterial secretion system", "Valine, leucine and isoleucine degradation"]
print(kegg[["layer3", "ID", "Input.number", "Background.number", "Pvalue", "FDR"]].head(8).to_string(index=False))
unc = set(kegg[kegg.Pvalue < 0.05].layer3); cor = set(kegg[kegg.FDR < 0.05].layer3)
print("\nUncorrected p < 0.05:", len(unc), "pathways; the five named among them:", set(named) <= unc, "; others:", sorted(unc - set(named)))
print("Corrected (FDR) < 0.05:", len(cor), "pathways; smallest FDR", round(kegg.FDR.min(), 4))
layer3 ID Input.number Background.number Pvalue FDR
Propanoate metabolism ko00640 8 38 0.002448 0.195837
RNA degradation ko03018 5 20 0.007548 0.301925
Synthesis and degradation of ketone bodies ko00072 3 10 0.022841 0.609084
Bacterial secretion system ko03070 8 61 0.040914 0.818280
Valine, leucine and isoleucine degradation ko00280 5 31 0.046383 0.742135
Tryptophan metabolism ko00380 4 23 0.057225 0.763001
beta-Alanine metabolism ko00410 3 14 0.057238 0.654152
Melanogenesis ko04916 1 1 0.064680 0.646804
Uncorrected p < 0.05: 5 pathways; the five named among them: True ; others: []
Corrected (FDR) < 0.05: 0 pathways; smallest FDR 0.1958
Result: blocked, methods underspecified. In the dataset's KEGG table the five named pathways are exactly those with uncorrected p < 0.05 (0.0024 to 0.0464), but none has FDR < 0.05 (smallest 0.1958). The criterion's two readings give different outcomes. The match of the named five with uncorrected p suggests the paper used uncorrected p, but the criterion leaves that open and blocks the claim when the readings differ; this attempt asks the question in questions/dryad-tdz08kq5d.yaml.
rnaseq-mutated-genes-expression¶
Under phage exposure phcR and phcS were upregulated and pilD and other pilus and fimbrial genes downregulated; seven of the 13 genes mutated under phage selection were differentially expressed.
Criterion: thresholds FDR < 0.05 and |log2FC| > 1; the 13 genes are those of Figure 4B. Reproduced when phcR and phcS are up, pilD is down, and seven of the 13 meet the thresholds. The dataset's Figure 4B list holds 11 genes; the genes mutated in +P clones by Table 1 are shown as well.
fb = fig4b.copy()
fb["DE"] = (fb.padj < 0.05) & (fb.log2FoldChange.abs() > 1)
print(fb.to_string(index=False))
print(f"\nFigure 4B list: {fb.DE.sum()} of {len(fb)} meet FDR < 0.05 and |log2FC| > 1; with FDR < 0.05 alone {(fb.padj < 0.05).sum()}")
print("At most", fb.DE.sum() + 2, "of 13 could meet both thresholds if the two genes missing from the list did")
plus_tags = var[var.group.str.contains(r"\+P") & var.LOCUS_TAG.notna()].LOCUS_TAG.unique()
t1 = full.reindex([pt.rnaseq_id(t) for t in plus_tags])[["log2FoldChange", "padj"]]
t1["gene"] = [names.get(t, t) for t in plus_tags]
t1["DE"] = (t1.padj < 0.05) & (t1.log2FoldChange.abs() > 1)
print(f"\nGenes with a variant in a +P group of the variant table ({len(t1)}):")
print(t1.round(4).to_string())
print("Meeting both thresholds:", int(t1.DE.sum()), "of", len(t1))
for gname in ("phcR", "phcS", "pilD"):
r = fb[fb.label.str.endswith(gname)].iloc[0]
print(f"{gname}: log2FC {r.log2FoldChange:.3f}, padj {r.padj:.2e}, meets both thresholds: {r.DE}")
print("\nOther pilus and fimbrial genes of the Figure 4C list:")
print(fig4c[fig4c.label.str.contains("pil|fim|Flp|cpa", case=False)][["label", "log2FoldChange", "padj"]].to_string(index=False))
gene label log2FoldChange padj DE
RS-N003355 0_RSc2889 0.291628 7.636127e-02 False
RS-N002118 1_RSc1770 -0.296235 1.238303e-01 False
RS-N002189 2_RSc1771 0.317011 4.940984e-01 False
RS-N001894 3_RSc1365 1.089598 7.891401e-10 True
RS-N003191 4_phcS 0.411209 3.043844e-03 False
RS-N003289 5_pilD -1.348558 6.525955e-36 True
RS-N003094 6_RSc2653 -0.451678 5.067883e-03 False
RS-N000359 7_gspE -0.048754 8.598386e-01 False
RS-N000383 10_RSc3091 0.057861 7.856951e-01 False
RS-N000358 11_gspF 0.488300 8.604830e-03 False
RS-N003192 12_phcR 1.834769 8.035238e-28 True
Figure 4B list: 3 of 11 meet FDR < 0.05 and |log2FC| > 1; with FDR < 0.05 alone 6
At most 5 of 13 could meet both thresholds if the two genes missing from the list did
Genes with a variant in a +P group of the variant table (18):
log2FoldChange padj gene DE
gene
RS-N002118 -0.2962 0.1238 RSc1770 False
RS-N002189 0.3170 0.4941 RSc1701 False
RS-N000287 0.6904 0.1797 RSc3188 False
RS-N003191 0.4112 0.0030 phcS False
RS-N003192 1.8348 0.0000 phcR True
RS-N003355 0.2916 0.0764 RSc2889 False
RS-N003435 1.2244 0.0000 gltB True
RS-N000358 0.4883 0.0086 gspF False
RS-N000359 -0.0488 0.8598 gspE False
RS-N005118 -0.1496 0.5061 RSc0104 False
RS-N003478 1.6731 0.0000 tISRso1 True
RS-N004190 1.5832 0.1345 RSc1775 False
RS-N001695 0.4524 0.0000 RSc1174 False
RS-N001894 1.0896 0.0000 RSc1365 True
RS-N003289 -1.3486 0.0000 pilD True
RS-N003706 0.3868 0.1458 RSp0796 False
RS-N003094 -0.4517 0.0051 RSc2653 False
RS-N000383 0.0579 0.7857 RSc3091 False
Meeting both thresholds: 5 of 18
phcR: log2FC 1.835, padj 8.04e-28, meets both thresholds: True
phcS: log2FC 0.411, padj 3.04e-03, meets both thresholds: False
pilD: log2FC -1.349, padj 6.53e-36, meets both thresholds: True
Other pilus and fimbrial genes of the Figure 4C list:
label log2FoldChange padj
type II and III secretion system protein family protein _RS-N003916_cpaC, rcpA -4.405639 5.800000e-186
prepilin peptidase_RS-N003919_cpaA, tadV -5.942366 0.000000e+00
Flp family type IVb pilin_RS-N003920_flp, pilA -5.625140 8.380000e-273
pilus assembly protein_RS-N003915_cpaE, tadZ -4.738090 8.540000e-126
Flp pilus assembly protein_RS-N003917_cpaB, rcpC -4.396537 1.730000e-124
Sigma-fimbriae chaperone protein_RS-N004603_fimC -2.100763 3.900000e-18
Result: partially reproduced. phcR is up (log2FC 1.835, padj 8e-28) and pilD down (-1.349, padj 7e-36), and the pilus and fimbrial genes of Figure 4C are down (-2.10 to -5.94). phcS has padj 0.003 but log2FC 0.411, below the threshold of 1. Of the 11 genes of the dataset's Figure 4B list, 3 meet both thresholds (6 meet padj < 0.05 alone), so at most 5 of 13 could, against the paper's 7. phcR up and pilD and the other pilus and fimbrial genes down hold; phcS falls short of the fold change threshold and the count of differentially expressed mutated genes (3, at most 5) falls short of 7. The criterion defines no partial outcome; the claim holds in part, so it is partially reproduced (changed from not reproduced after verification round 1).
rnaseq-defense-systems¶
Phage exposure upregulated type I restriction-modification genes (MTases, REases) and CBASS genes (SMODS cyclase, Gros effector) and downregulated the Wadjet gene JetC.
Criterion: thresholds padj < .05 and |log2FC| > 2; the loci are those annotated as these system genes. The dataset's Figure 4C list names them with their loci; values are read from the main table.
defense = fig4c[fig4c.label.str.contains("JetC|MTase|REase|SMODS|Effector")].copy()
defense = defense.join(full[["log2FoldChange", "padj"]].add_suffix("_4A"), on="gene")
defense["meets"] = (defense.padj_4A < 0.05) & (defense.log2FoldChange_4A.abs() > 2)
print(defense[["label", "gene", "log2FoldChange_4A", "padj_4A", "regulation", "meets"]].to_string(index=False))
label gene log2FoldChange_4A padj_4A regulation meets
RS-N004717_JetC_II RS-N004717 -4.209067 6.920000e-65 DOWN True
RS-N000213_Type_I_MTases RS-N000213 2.302917 2.170000e-54 UP True
RS-N000216_Effector_1TM_DUF1043 RS-N000216 2.324671 1.500000e-35 UP True
RS-N000215_Type_I_REases RS-N000215 2.600479 2.440000e-26 UP True
RS-N000217_Cyclase_SMODS RS-N000217 2.640740 1.760000e-44 UP True
Result: reproduced. JetC is down (log2FC -4.21), and the type I MTase (2.30), REase (2.60), SMODS cyclase (2.64) and the CBASS effector locus (2.32) are up, all with padj < 1e-25. The dataset labels the effector Effector_1TM_DUF1043, not Gros; it is the locus beside the SMODS cyclase.
rnaseq-virulence-genes-down¶
Phage exposure downregulated motility genes and type II (e.g. gspG, rcpA, tadC) and type III secretion system genes (several hrp and hrc genes).
Criterion: thresholds padj < .05 and |log2FC| > 2; reproduced when gspG, rcpA and tadC and at least two hrp or hrc genes meet the thresholds downward and none of these is upregulated.
vir = fig4c[fig4c.label.str.contains("gspG|rcpA|tadC|hrp|hrc", case=False)].copy()
vir = vir.join(full[["log2FoldChange", "padj"]].add_suffix("_4A"), on="gene")
vir["meets_down"] = (vir.padj_4A < 0.05) & (vir.log2FoldChange_4A < -2)
vir["up"] = (vir.padj_4A < 0.05) & (vir.log2FoldChange_4A > 0)
print(vir[["label", "gene", "log2FoldChange_4A", "padj_4A", "regulation", "meets_down", "up"]].to_string(index=False))
print("\nhrp/hrc genes meeting the thresholds downward:", int(vir[vir.label.str.contains("hrp|hrc")].meets_down.sum()))
label gene log2FoldChange_4A padj_4A regulation meets_down up
type II and III secretion system protein family protein _RS-N003916_cpaC, rcpA RS-N003916 -4.405639 5.800000e-186 DOWN True False
type II secretion system F family protein_RS-N003912_tadC RS-N003912 -4.266422 1.790000e-130 DOWN True False
type III secretion protein_RS-N003622_hrpB2 RS-N003622 -2.490343 1.790000e-17 DOWN True False
type III secretion system apparatus protein_RS-N003623_hrpB1 RS-N003623 -2.371315 2.050000e-10 DOWN True False
type III secretion system export apparatus protein EscS/YscS/HrcS family_RS-N003629_yscS, sctS, hrcS RS-N003629 -2.065348 4.660000e-05 DOWN True False
type III secretion system export apparatus protein EscV/YscV/HrcV family_RS-N003625_yscV, sctV, hrcV RS-N003625 -2.055090 1.210000e-19 DOWN True False
type III secretion system export apparatus protein EscR/YscR/HrcR family _RS-N003628_yscR, sctR, hrcR RS-N003628 -2.048767 1.190000e-07 DOWN True False
type III secretion apparatus protein YscQ/HrcQ family_RS-N003627_yscQ, sctQ, hrcQ RS-N003627 -2.001040 9.110000e-09 DOWN True False
type II secretion system protein_RS-N004330_gspG RS-N004330 2.037027 1.527727e-03 UP False True
hrp/hrc genes meeting the thresholds downward: 6
Result: partially reproduced. rcpA (log2FC -4.41) and tadC (-4.27) are down, and six hrp or hrc genes (hrpB2, hrpB1, hrcS, hrcV, hrcR, hrcQ, log2FC -2.00 to -2.49) meet the thresholds downward. gspG, which the text names as downregulated, is up in the dataset's table (log2FC 2.04, padj 0.0015) and labelled UP in its Figure 4C list.
Summary¶
| Claim | Outcome |
|---|---|
| fig1b-temperature-density | reproduced |
| fig1b-phage-reduces-density | reproduced |
| fig1b-phage-effect-diminishes | reproduced |
| auc-phage-reduction-by-temperature | reproduced |
| fig1c-phage-density-temperature | reproduced |
| fig1d-mucoid-lost-with-phage | reproduced |
| fig1d-phage-raises-scv-frequency | partially reproduced |
| fig2a-phage-exposed-more-resistant | reproduced |
| fig2a-temperature-resistance | reproduced |
| fig2a-no-temperature-phage-interaction | reproduced |
| fig2a-scv-more-resistant | reproduced |
| fig2a-mucoid-25p-replicate-resistant | reproduced |
| fig2b-scv-higher-growth-rate | reproduced |
| fig2c-scv-lower-carrying-capacity | reproduced |
| fig2d-scv-twitching-no-exposure-effect | reproduced |
| fig2d-mucoid-25p-reduced-twitching | reproduced |
| fig2e-mucoid-biofilm-increased | partially reproduced |
| fig2e-biofilm-no-exposure-effect | reproduced |
| fig2f-phage-exposure-reduces-virulence | reproduced |
| fig2f-scv-avirulent | reproduced |
| fig2f-mucoid-25p-avirulent | reproduced |
| sequencing-mutation-count | not reproduced |
| table1-phcs-in-scvs | reproduced |
| table1-pild-in-25p-mucoid | reproduced |
| table1-t2ss-in-35p-scv | blocked: metadata insufficient |
| tables6-pil-mutants-resistant | blocked: data missing |
| figs3-adsorption | reproduced |
| rnaseq-differential-expression-count | reproduced |
| rnaseq-kegg-enrichment | blocked: methods underspecified |
| rnaseq-mutated-genes-expression | partially reproduced |
| rnaseq-defense-systems | reproduced |
| rnaseq-virulence-genes-down | partially reproduced |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. The growth parameters, the variant calls, the differential expression and the KEGG enrichment were not recomputed: the dataset holds the authors' fitted growth rates and carrying capacities (no growth curves), their variant calls aggregated by group, their DESeq2 table and their KOBAS table, and the claims are tested on these. The ANOVAs are computed in Python (statsmodels) where the paper used R 3.6.3. The colony types of the twitching sheet are a choice this attempt had to make; it is shown against the reading by strain name. The OD time points of Figure 1B, the readings of the 0.24% and of the area under the curve, and the unit of the Figure 2A to 2C ANOVAs were questions of the first run; the claims file now settles them (days 2 to 12; the mean of replicate frequencies or pooled counts; the readings by day; unit means, with clones as a variant), and the alternatives are kept as variants. The sequencing and RNA-seq claims are data: dataset after the same answers, tested on the deposited tables.
What was found about the dataset. Where the paper prints a test, the deposited data reproduce it: Tables S1a, S1b, S2, S3a to S3c, S3e and S3f, the F values of the area under the curve, the two one-way ANOVAs of Figure 2A and 2D, and the counts of Table S4. The matches to Tables S1a and S3a to S3c identify the paper's choices (days 2 to 12; replicate means by colony type with sequential sums of squares). The exceptions are the twitching sheet, whose values do not reproduce Table S3d (read as sequential sums of squares) and whose strain names do not identify the clones (four 35_1 rows carry 35_0 names, and four 25_0 names appear twice), and the paper's 20% mucoid at 25 degrees C without phage, which the data put at 15.3%. The variant table is aggregated by group, so per-clone counts cannot be formed, and the genome reads of PRJNA1028400 are not public (28 BioSamples, 14 assemblies, no runs). The RNA-seq reads are public in PRJNA1032726, though the Methods say they are shared on request. The DESeq2 table reproduces the paper's counts exactly, while phcS (log2FC 0.41) and gspG (up, not down) differ from the text. The SCVs form more biofilm than the ancestor in every clone. The pilB and pilM data of Table S6 are not deposited. Minor: the carrying capacity sheet heads its ancestor value maxgrowthrate(OD600), OD600 values are blank-corrected without saying so, and the PROKKA gene names of the variant table differ from the paper's (a mapping sheet is given).
print(f"Notebook ran in {time.time() - started:.0f} s")
Notebook ran in 4 s