AMR co-occurrence networks are robust-yet-fragile: resilient to random drug pressure but fragile to targeted hub disruption

Antimicrobial resistance (AMR) gene co-occurrence networks, constructed from metagenomic or clinical isolate datasets, exhibit a "robust-yet-fragile" (RYF) topology: when nodes (AMR genes or resistance classes) are removed randomly (simulating stochastic or broad-spectrum antibiotic pressure), network connectivity and functional redundancy are preserved (>50% of the largest connected component retained after removing 30% of nodes randomly); but when high-degree hub nodes are removed in targeted fashion (simulating precision disruption of co-selected resistance determinants), network fragmentation occurs disproportionately (>60% reduction in largest connected component after removing only the top 10% of nodes by degree/betweenness centrality). This RYF property is statistically distinguishable from Erdős–Rényi random networks and scale-free null models at p < 0.05.

PHASE 1 — Network Construction (Days 1–15): Collect ≥4 independent AMR metagenomic datasets. Detect AMR genes using AMRFinderPlus or CARD RGI. Compute pairwise co-occurrence using SparCC or SPIEC-EASI. Apply FDR correction (Benjamini-Hochberg, q < 0.05). Construct weighted undirected networks. Characterize topology (degree distribution, clustering coefficient, modularity).

PHASE 2 — Robustness Simulation (Days 16–30): Implement sequential node removal: (a) random order, (b) descending degree, (c) descending betweenness centrality, (d) descending eigenvector centrality. At each removal step (1% increments), compute: largest connected component (LCC) size, average path length, network diameter, algebraic connectivity. Run 1,000 bootstrap resamples per strategy. Compare attack curves using area-under-curve (AUC) analysis.

PHASE 3 — Null Model Comparison (Days 31–40): Generate 1,000 Erdős–Rényi random networks matched for N nodes and M edges. Generate 1,000 Barabási–Albert scale-free networks matched for degree sequence. Repeat robustness simulations on all null models. Compute empirical p-values for RYF differential.

PHASE 4 — Biological Validation (Days 41–60): Cross-reference hub genes with CARD database for MGE association. Query PATRIC/NCBI for co-occurrence of hub genes on same plasmid/integron. Validate 3–5 top hub gene pairs in conjugation transfer assays (E. coli K-12 donor/recipient system). Measure transfer frequency with and without competing AMR gene presence.

# AMR Co-occurrence Network RYF Validation Pipeline
# Architecture: Modular Python pipeline with R integration for SparCC

import numpy as np
import pandas as pd
import networkx as nx
from scipy import stats
import powerlaw
from itertools import combinations

# ============================================================
# MODULE 1: DATA LOADING AND AMR PROFILING
# ============================================================
class AMRProfiler:
    def __init__(self, dataset_paths, card_db_path):
        self.datasets = dataset_paths  # dict: {name: path}
        self.card_db = card_db_path
    
    def run_rgi(self, sample_fasta):
        # Shell call: rgi main --input_sequence sample.fasta
        #             --output_file output --input_type contig
        #             --alignment_tool DIAMOND --include_loose
        pass
    
    def build_abundance_matrix(self, rgi_outputs):
        # Returns DataFrame: samples × AMR_gene_families (RPKM)
        # Filter: retain genes in ≥5% of samples
        matrix = pd.DataFrame()  # shape: (n_samples, n_genes)
        prevalence_filter = matrix.apply(lambda x: (x > 0).mean()) >= 0.05
        return matrix.loc[:, prevalence_filter]

# ============================================================
# MODULE 2: CO-OCCURRENCE NETWORK CONSTRUCTION
# ============================================================
class CooccurrenceNetworkBuilder:
    def __init__(self, corr_threshold=0.3, fdr_alpha=0.05, n_bootstraps=1000):
        self.corr_threshold = corr_threshold
        self.fdr_alpha = fdr_alpha
        self.n_bootstraps = n_bootstraps
    
    def compute_sparcc(self, abundance_matrix):
        # R call via rpy2: SparCC(abundance_matrix, iter=20, 
        #                         bootstraps=self.n_bootstraps)
        # Returns: correlation_matrix, pvalue_matrix
        correlations = np.zeros((n_genes, n_genes))
        pvalues = np.ones((n_genes, n_genes))
        return correlations, pvalues
    
    def build_network(self, correlations, pvalues):
        # Apply BH FDR correction
        from statsmodels.stats.multitest import multipletests
        flat_pvals = pvalues[np.triu_indices_from(pvalues, k=1)]
        _, corrected_pvals, _, _ = multipletests(flat_pvals, method='fdr_bh')
        
        G = nx.Graph()
        idx = 0
        for i, j in combinations(range(n_genes), 2):
            if (abs(correlations[i,j]) > self.corr_threshold and 
                corrected_pvals[idx] < self.fdr_alpha):
                G.add_edge(gene_names[i], gene_names[j], 
                          weight=correlations[i,j])
            idx += 1
        return G
    
    def characterize_topology(self, G):
        degrees = [d for n, d in G.degree()]
        fit = powerlaw.Fit(degrees, discrete=True)
        # Test power law vs exponential
        R, p = fit.distribution_compare('power_law', 'exponential')
        return {
            'n_nodes': G.number_of_nodes(),
            'n_edges': G.number_of_edges(),
            'alpha': fit.power_law.alpha,
            'xmin': fit.power_law.xmin,
            'powerlaw_vs_exp_R': R,
            'powerlaw_vs_exp_p': p,
            'clustering': nx.average_clustering(G),
            'modularity': compute_modularity_louvain(G)
        }

# ============================================================
# MODULE 3: HUB IDENTIFICATION
# ============================================================
class HubIdentifier:
    def __init__(self, top_fraction=0.10):
        self.top_fraction = top_fraction
    
    def compute_centralities(self, G):
        n_hubs = max(1, int(G.number_of_nodes() * self.top_fraction))
        centralities = {
            'degree': nx.degree_centrality(G),
            'betweenness': nx.betweenness_centrality(G, k=500, normalized=True),
            'eigenvector': nx.eigenvector_centrality(G, max_iter=1000),
            'pagerank': nx.pagerank(G, alpha=0.85)
        }
        # Identify consensus hubs: top 10% in ≥2 metrics
        hub_votes = {node: 0 for node in G.nodes()}
        for metric, scores in centralities.items():
            threshold = sorted(scores.values())[-n_hubs]
            for node, score in scores.items():
                if score >= threshold:
                    hub_votes[node] += 1
        consensus_hubs = [n for n, v in hub_votes.items() if v >= 2]
        return centralities, consensus_hubs

# ============================================================
# MODULE 4: ROBUSTNESS SIMULATION
# ============================================================
class RobustnessSimulator:
    def __init__(self, n_bootstraps=1000, removal_steps=100):
        self.n_bootstraps = n_bootstraps
        self.removal_steps = removal_steps
    
    def random_attack(self, G):
        """Remove nodes uniformly at random, track LCC size"""
        results = []
        for _ in range(self.n_bootstraps):
            G_copy = G.copy()
            nodes = list(G_copy.nodes())
            np.random.shuffle(nodes)
            lcc_trajectory = [1.0]  # normalized LCC size
            n_total = len(nodes)
            step_size = max(1, n_total // self.removal_steps)
            
            for i in range(0, n_total, step_size):
                remove_batch = nodes[i:i+step_size]
                G_copy.remove_nodes_from(remove_batch)
                if G_copy.number_of_nodes() > 0:
                    lcc = max(nx.connected_components(G_copy), key=len)
                    lcc_trajectory.append(len(lcc) / n_total)
                else:
                    lcc_trajectory.append(0.0)
            results.append(lcc_trajectory)
        return np.array(results)  # shape: (n_bootstraps, removal_steps+1)
    
    def targeted_attack(self, G, centrality_metric='degree', adaptive=True):
        """Remove nodes in descending centrality order"""
        results = []
        n_total = G.number_of_nodes()
        step_size = max(1, n_total // self.removal_steps)
        
        for _ in range(self.n_bootstraps):
            G_copy = G.copy()
            lcc_trajectory = [1.0]
            
            removed = 0
            while G_copy.number_of_nodes() > 0:
                if adaptive or removed == 0:
                    # Recompute centrality after each removal
                    if centrality_metric == 'degree':
                        scores = dict(G_copy.degree())
                    elif centrality_metric == 'betweenness':
                        scores = nx.betweenness_centrality(G_copy, k=min(100, G_copy.number_of_nodes()))
                    
                    sorted_nodes = sorted(scores, key=scores.get, reverse=True)
                
                remove_batch = sorted_nodes[:step_size]
                sorted_nodes = sorted_nodes[step_size:]
                G_copy.remove_nodes_from(remove_batch)
                removed += len(remove_batch)
                
                if G_copy.number_of_nodes() > 0:
                    lcc = max(nx.connected_components(G_copy), key=len)
                    lcc_trajectory.append(len(lcc) / n_total)
                else:
                    lcc_trajectory.append(0.0)
            
            results.append(lcc_trajectory)
        return np.array(results)
    
    def compute_auc(self, lcc_trajectory):
        """Compute area under LCC-vs-removal curve (trapezoidal)"""
        x = np.linspace(0, 1, len(lcc_trajectory))
        return np.trapz(lcc_trajectory, x)
    
    def compute_delta_auc(self, random_results, targeted_results):
        """Primary test statistic: ΔAUC = AUC(random) - AUC(targeted)"""
        auc_random = np.array([self.compute_auc(r) for r in random_results])
        auc_targeted = np.array([self.compute_auc(t) for t in targeted_results])
        delta_auc = auc_random.mean() - auc_targeted.mean()
        # Bootstrap CI
        delta_samples = auc_random - auc_targeted[:len(auc_random)]
        ci_low, ci_high = np.percentile(delta_samples, [2.5, 97.5])
        return delta_auc, ci_low, ci_high

# ============================================================
# MODULE 5: NULL MODEL COMPARISON
# ============================================================
class NullModelComparator:
    def __init__(self, n_null_networks=1000):
        self.n_null = n_null_networks
    
    def generate_erdos_renyi(self, n, m):
        p = 2 * m / (n * (n - 1))
        return nx.erdos_renyi_graph(n, p)
    
    def generate_configuration_model(self, degree_sequence):
        return nx.configuration_model(degree_sequence)
    
    def generate_barabasi_albert(self, n, m_param):
        return nx.barabasi_albert_graph(n, m_param)
    
    def compute_null_delta_aucs(self, G, simulator, model_type='configuration'):
        null_delta_aucs = []
        degrees = [d for _, d in G.degree()]
        
        for _ in range(self.n_null):
            if model_type == 'configuration':
                G_null = self.generate_configuration_model(degrees)
                G_null = nx.Graph(G_null)  # remove multi-edges
            elif model_type == 'erdos_renyi':
                G_null = self.generate_erdos_renyi(G.number_of_nodes(), G.number_of_edges())
            
            # Run abbreviated simulation (100 bootstraps for null)
            sim_null = RobustnessSimulator(n_bootstraps=100)
            rand_results = sim_null.random_attack(G_null)
            targ_results = sim_null.targeted_attack(G_null)
            delta, _, _ = sim_null.compute_delta_auc(rand_results, targ_results)
            null_delta_aucs.append(delta)
        
        return np.array(null_delta_aucs)
    
    def empirical_pvalue(self, observed_delta, null_deltas):
        return np.mean(null_deltas >= observed_delta)

# ============================================================
# MODULE 6: META-ANALYSIS
# ============================================================
class MetaAnalyzer:
    def random_effects_meta(self, effect_sizes, variances):
        """DerSimonian-Laird random effects meta-analysis"""
        # Fixed effects weights
        w_fixed = 1.0 / np.array(variances)
        theta_fixed = np.sum(w_fixed * effect_sizes) / np.sum(w_fixed)
        
        # Cochran's Q and tau^2
        Q = np