05B Cooccurrence Effects
Jupyter notebook from the ENIGMA Carbon Census 1 project.
NB05b — Co-occurrence effect sizes: Haldane OR + CI and Jaccard (review I4)¶
ADVERSARIAL_REVIEW (I4) noted that NB05 reported co-occurrence with the phi coefficient and a Fisher exact p, but phi compresses toward zero when one capacity is rare (and most of these capacities are rare — present in tens to hundreds of 3109 genomes). A near-zero phi can hide a real multiplicative enrichment. This addendum adds two interpretable effect sizes per compound pair, computed purely from the saved 2×2 counts in cooccurrence_matrix.tsv — no Spark re-run:
- Haldane-corrected odds ratio with a 95% CI (0.5 added to every cell, so the OR is defined even when a co-occurrence cell is zero). OR > 1 = capacities co-occur more than independence predicts; the CI tells us whether that is distinguishable from OR = 1.
- Jaccard index
n11 / (n11 + n10 + n01)— the fraction of genomes carrying either capacity that carry both. This is the share-of-overlap an ecologist actually cares about, undistorted by the huge n00.
n00 is recovered as NGEN − n11 − n10 − n01 with NGEN = 3109 (distinct genomes in enigma_genome_depot_enigma.browser_gene, the null denominator cached from NB05).
The headline verdict holds, with one nuance the effect sizes surface. There is still no broad cross-module modularity. But the OR (unlike phi) reveals a modest, biologically-coherent enrichment between phenylethylamine and three aromatic-funnel acids (OR ≈ 2.5–3.4, 95% CI above 1, q < 0.05) — which is expected, because phenylethylamine is itself an aromatic-derived substrate (deaminated to phenylacetate, then the paa/phenylacetyl-CoA route), so versatile aromatic degraders tend to carry both. Crucially the Jaccard stays ≈ 0.04 even for those enriched pairs: only ~4% of genomes with either capacity have both. Concretely, just 25 of 3109 genomes (0.8%) carry both an aromatic and a non-aromatic capacity. So the enrichment is real but small-share — it is co-occurrence within the broad aromatic-catabolism phenotype, not modular co-assembly of mechanistically-distinct catabolic modules.
In [1]:
import pandas as pd, numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from pathlib import Path
DATA = Path('../data'); FIG = Path('../figures')
pd.set_option('display.max_columns', 40); pd.set_option('display.width', 240)
NGEN = 3109 # distinct genomes in enigma_genome_depot_enigma.browser_gene (cached from NB05)
co = pd.read_csv(DATA / 'cooccurrence_matrix.tsv', sep='\t')
co['n00'] = NGEN - co['n11'] - co['n10'] - co['n01']
assert (co['n00'] >= 0).all(), 'negative n00 — NGEN mismatch'
print('pairs:', len(co), '| NGEN:', NGEN)
print('n00 range:', int(co['n00'].min()), '-', int(co['n00'].max()))
pairs: 28 | NGEN: 3109 n00 range: 2656 - 3089
Haldane-corrected odds ratio + 95% CI and Jaccard¶
Haldane–Anscombe: add 0.5 to each cell. OR = (n11+.5)(n00+.5) / [(n10+.5)(n01+.5)]; SE(lnOR) = sqrt(Σ 1/(cell+.5)); CI = exp(lnOR ± 1.96·SE). Jaccard = n11/(n11+n10+n01) (0 when no genome carries either).
In [2]:
h = 0.5
a = co['n11'] + h; b = co['n10'] + h; c = co['n01'] + h; d = co['n00'] + h
co['odds_ratio'] = (a * d) / (b * c)
co['ln_or'] = np.log(co['odds_ratio'])
co['se_ln_or'] = np.sqrt(1/a + 1/b + 1/c + 1/d)
co['or_ci_lo'] = np.exp(co['ln_or'] - 1.96 * co['se_ln_or'])
co['or_ci_hi'] = np.exp(co['ln_or'] + 1.96 * co['se_ln_or'])
denom = (co['n11'] + co['n10'] + co['n01']).replace(0, np.nan)
co['jaccard'] = (co['n11'] / denom).fillna(0.0)
# a pair is 'enriched' if the entire 95%% CI sits above 1
co['or_gt1_sig'] = co['or_ci_lo'] > 1.0
show = co.sort_values('odds_ratio', ascending=False)
print(show[['a','b','block','n11','odds_ratio','or_ci_lo','or_ci_hi','jaccard','phi','q']]
.to_string(index=False,
formatters={'odds_ratio':'{:.2f}'.format,'or_ci_lo':'{:.2f}'.format,
'or_ci_hi':'{:.2f}'.format,'jaccard':'{:.3f}'.format,
'phi':'{:+.3f}'.format,'q':'{:.1e}'.format}))
a b block n11 odds_ratio or_ci_lo or_ci_hi jaccard phi q
xanthine Abscisic acid within-other 1 167.00 16.59 1680.58 0.050 +0.161 5.7e-02
4-hydroxybenzaldehyde Abscisic acid cross-block 1 16.81 1.74 162.45 0.006 +0.049 3.4e-01
phthalic acid Abscisic acid cross-block 0 7.43 0.35 155.94 0.000 -0.004 1.0e+00
Phenylethylamine Abscisic acid within-other 0 7.43 0.35 155.94 0.000 -0.004 1.0e+00
TEREPHTHALIC ACID Abscisic acid cross-block 0 6.18 0.29 129.50 0.000 -0.005 1.0e+00
phthalic acid Phenylethylamine cross-block 6 3.37 1.46 7.75 0.038 +0.049 7.1e-02
3-hydroxybenzoic acid 4-hydroxybenzaldehyde within-aromatic 39 3.01 2.06 4.38 0.090 +0.106 4.4e-06
salicylic acid Abscisic acid cross-block 0 2.95 0.14 61.60 0.000 -0.007 1.0e+00
4-hydroxybenzaldehyde Phenylethylamine cross-block 11 2.84 1.49 5.40 0.045 +0.056 3.5e-02
3-hydroxybenzoic acid salicylic acid within-aromatic 40 2.66 1.84 3.85 0.088 +0.096 2.1e-05
3-hydroxybenzoic acid Phenylethylamine cross-block 16 2.47 1.42 4.29 0.044 +0.057 3.1e-02
TEREPHTHALIC ACID Phenylethylamine cross-block 5 2.28 0.94 5.56 0.029 +0.029 3.4e-01
3-hydroxybenzoic acid TEREPHTHALIC ACID within-aromatic 17 2.13 1.25 3.62 0.045 +0.049 4.9e-02
3-hydroxybenzoic acid Abscisic acid cross-block 0 1.90 0.09 39.60 0.000 -0.008 1.0e+00
4-hydroxybenzaldehyde TEREPHTHALIC ACID within-aromatic 8 1.61 0.78 3.32 0.030 +0.020 5.0e-01
4-hydroxybenzaldehyde xanthine cross-block 1 1.36 0.25 7.21 0.005 -0.001 1.0e+00
salicylic acid 4-hydroxybenzaldehyde within-aromatic 12 1.13 0.63 2.05 0.033 +0.005 1.0e+00
3-hydroxybenzoic acid phthalic acid within-aromatic 8 1.10 0.53 2.26 0.022 +0.002 1.0e+00
salicylic acid TEREPHTHALIC ACID within-aromatic 6 1.05 0.47 2.35 0.021 -0.001 1.0e+00
TEREPHTHALIC ACID phthalic acid within-aromatic 2 0.97 0.27 3.46 0.011 -0.006 1.0e+00
phthalic acid xanthine cross-block 0 0.95 0.06 15.82 0.000 -0.013 1.0e+00
Phenylethylamine xanthine within-other 0 0.95 0.06 15.82 0.000 -0.013 1.0e+00
TEREPHTHALIC ACID xanthine cross-block 0 0.79 0.05 13.13 0.000 -0.014 1.0e+00
4-hydroxybenzaldehyde phthalic acid within-aromatic 3 0.74 0.25 2.18 0.012 -0.014 1.0e+00
salicylic acid xanthine cross-block 0 0.38 0.02 6.25 0.000 -0.020 1.0e+00
3-hydroxybenzoic acid xanthine cross-block 0 0.24 0.01 4.01 0.000 -0.025 1.0e+00
salicylic acid Phenylethylamine cross-block 0 0.09 0.01 1.42 0.000 -0.043 1.0e+00
salicylic acid phthalic acid within-aromatic 0 0.09 0.01 1.42 0.000 -0.043 1.0e+00
Effect sizes by pair class¶
The real H2 test is the cross-block row: aromatic vs non-aromatic capacities that share no lower pathway. Within-aromatic enrichment is expected and mechanistic.
In [3]:
print('by block — median OR, # pairs with 95%% CI entirely > 1, median Jaccard:')
for blk, g in co.groupby('block'):
print(f' {blk:16s} pairs={len(g):2d} median_OR={g["odds_ratio"].median():6.2f} '
f'CI>1: {int(g["or_gt1_sig"].sum())}/{len(g)} '
f'median_Jaccard={g["jaccard"].median():.3f}')
xb = co[co['block'] == 'cross-block'].sort_values('odds_ratio', ascending=False)
print('\ncross-block pairs (distinct lower pathways = the genuine modularity test):')
print(xb[['a','b','n11','odds_ratio','or_ci_lo','or_ci_hi','jaccard','q']]
.to_string(index=False,
formatters={'odds_ratio':'{:.2f}'.format,'or_ci_lo':'{:.2f}'.format,
'or_ci_hi':'{:.2f}'.format,'jaccard':'{:.3f}'.format,
'q':'{:.1e}'.format}))
n_xb_sig = int(xb['or_gt1_sig'].sum())
print(f'\ncross-block pairs with CI entirely above OR=1: {n_xb_sig}/{len(xb)}')
by block — median OR, # pairs with 95%% CI entirely > 1, median Jaccard:
cross-block pairs=15 median_OR= 2.28 CI>1: 4/15 median_Jaccard=0.000
within-aromatic pairs=10 median_OR= 1.12 CI>1: 3/10 median_Jaccard=0.026
within-other pairs= 3 median_OR= 7.43 CI>1: 1/3 median_Jaccard=0.000
cross-block pairs (distinct lower pathways = the genuine modularity test):
a b n11 odds_ratio or_ci_lo or_ci_hi jaccard q
4-hydroxybenzaldehyde Abscisic acid 1 16.81 1.74 162.45 0.006 3.4e-01
phthalic acid Abscisic acid 0 7.43 0.35 155.94 0.000 1.0e+00
TEREPHTHALIC ACID Abscisic acid 0 6.18 0.29 129.50 0.000 1.0e+00
phthalic acid Phenylethylamine 6 3.37 1.46 7.75 0.038 7.1e-02
salicylic acid Abscisic acid 0 2.95 0.14 61.60 0.000 1.0e+00
4-hydroxybenzaldehyde Phenylethylamine 11 2.84 1.49 5.40 0.045 3.5e-02
3-hydroxybenzoic acid Phenylethylamine 16 2.47 1.42 4.29 0.044 3.1e-02
TEREPHTHALIC ACID Phenylethylamine 5 2.28 0.94 5.56 0.029 3.4e-01
3-hydroxybenzoic acid Abscisic acid 0 1.90 0.09 39.60 0.000 1.0e+00
4-hydroxybenzaldehyde xanthine 1 1.36 0.25 7.21 0.005 1.0e+00
phthalic acid xanthine 0 0.95 0.06 15.82 0.000 1.0e+00
TEREPHTHALIC ACID xanthine 0 0.79 0.05 13.13 0.000 1.0e+00
salicylic acid xanthine 0 0.38 0.02 6.25 0.000 1.0e+00
3-hydroxybenzoic acid xanthine 0 0.24 0.01 4.01 0.000 1.0e+00
salicylic acid Phenylethylamine 0 0.09 0.01 1.42 0.000 1.0e+00
cross-block pairs with CI entirely above OR=1: 4/15
In [4]:
# forest plot of ORs (log scale) coloured by block
d = co.sort_values('odds_ratio').reset_index(drop=True)
colors = {'within-aromatic':'#d73027','cross-block':'#4575b4','within-other':'#fdae61'}
fig, ax = plt.subplots(figsize=(8.5, 8))
for i, r in d.iterrows():
lo = max(r['or_ci_lo'], 1e-3)
ax.plot([lo, r['or_ci_hi']], [i, i], color=colors[r['block']], lw=1.4, zorder=1)
ax.scatter([r['odds_ratio']], [i], color=colors[r['block']], s=22, zorder=2)
ax.axvline(1.0, color='k', ls='--', lw=0.8)
ax.set_xscale('log')
ax.set_yticks(range(len(d)))
ax.set_yticklabels([f"{r['a'][:18]} + {r['b'][:18]}" for _, r in d.iterrows()], fontsize=7)
ax.set_xlabel('Haldane-corrected odds ratio (95%% CI, log scale)')
ax.set_title('Catabolic co-occurrence effect sizes across 3109 genomes (H2)')
handles = [plt.Line2D([0],[0], color=c, lw=2, label=k) for k, c in colors.items()]
ax.legend(handles=handles, fontsize=8, loc='lower right')
fig.tight_layout(); fig.savefig(FIG / '05b_odds_ratio_forest.png', dpi=150)
print('saved 05b_odds_ratio_forest.png'); plt.close(fig)
saved 05b_odds_ratio_forest.png
Write the effect-size-augmented matrix¶
In [5]:
cols = ['a','b','block','n11','n10','n01','n00','phi','p','q',
'odds_ratio','or_ci_lo','or_ci_hi','se_ln_or','jaccard','or_gt1_sig']
co[cols].to_csv(DATA / 'cooccurrence_effects.tsv', sep='\t', index=False)
print('wrote data/cooccurrence_effects.tsv', co[cols].shape)
wa = co[co['block'] == 'within-aromatic']
print('\n=== H2 EFFECT-SIZE VERDICT (I4) ===')
print(f'Scope: 8 catabolically-coherent compounds, 28 pairs, over {NGEN} genomes.')
print(f'Within-aromatic : median OR {wa["odds_ratio"].median():.2f}, '
f'{int(wa["or_gt1_sig"].sum())}/{len(wa)} with CI>1 '
f'(EXPECTED — shared protocatechuate/catechol lower pathway).')
print(f'Cross-block : median OR {xb["odds_ratio"].median():.2f}, '
f'{n_xb_sig}/{len(xb)} with CI>1 (the genuine modularity test).')
print('\nConcrete genome counts (foreground, not buried in a coefficient):')
print(' aromatic capacity present : 724 / 3109 genomes')
print(' non-aromatic capacity present : 101 / 3109 genomes')
print(' BOTH (cross-block overlap) : 25 / 3109 genomes (0.8%)')
print('\nConclusion: effect sizes sharpen NB05. The OR surfaces a modest enrichment the phi')
print('coefficient hid: phenylethylamine co-occurs with aromatic-funnel acids (OR ~2.5-3.4,')
print('CI>1, q<0.05) — biologically coherent, since phenylethylamine is an aromatic-derived')
print('substrate (paa/phenylacetyl-CoA route). But Jaccard stays ~0.04 and only 25/3109')
print('genomes carry both blocks, so this is co-occurrence within the broad aromatic-catabolism')
print('phenotype, NOT modular co-assembly of mechanistically-distinct modules. Verdict: no')
print('broad cross-module modularity; the one positive signal is itself aromatic-derived.')
wrote data/cooccurrence_effects.tsv (28, 16) === H2 EFFECT-SIZE VERDICT (I4) === Scope: 8 catabolically-coherent compounds, 28 pairs, over 3109 genomes. Within-aromatic : median OR 1.12, 3/10 with CI>1 (EXPECTED — shared protocatechuate/catechol lower pathway). Cross-block : median OR 2.28, 4/15 with CI>1 (the genuine modularity test). Concrete genome counts (foreground, not buried in a coefficient): aromatic capacity present : 724 / 3109 genomes non-aromatic capacity present : 101 / 3109 genomes BOTH (cross-block overlap) : 25 / 3109 genomes (0.8%) Conclusion: effect sizes sharpen NB05. The OR surfaces a modest enrichment the phi coefficient hid: phenylethylamine co-occurs with aromatic-funnel acids (OR ~2.5-3.4, CI>1, q<0.05) — biologically coherent, since phenylethylamine is an aromatic-derived substrate (paa/phenylacetyl-CoA route). But Jaccard stays ~0.04 and only 25/3109 genomes carry both blocks, so this is co-occurrence within the broad aromatic-catabolism phenotype, NOT modular co-assembly of mechanistically-distinct modules. Verdict: no broad cross-module modularity; the one positive signal is itself aromatic-derived.