#!/usr/bin/env python3 """ Offline network enrichment and interaction analysis using local STRING-DB files. Replicates the NETWORK_ENRICHMENT Nextflow process without internet connection. """ import pandas as pd import argparse from pathlib import Path from scipy import stats import numpy as np def load_string_data(string_db_dir): """Load all required STRING-DB files.""" string_db_dir = Path(string_db_dir) print("Loading STRING-DB files...") # Load alias mapping (ENST -> ENSP) aliases = pd.read_csv( string_db_dir / "9606.protein.aliases.enst.v12.0.tsv", sep='\t', comment='#', names=['string_protein_id', 'alias', 'source'] ) enst_to_ensp = aliases[aliases['source'] == 'Ensembl_transcript'].set_index('alias')['string_protein_id'].to_dict() # Load protein info (for preferred names) info = pd.read_csv( string_db_dir / "9606.protein.info.v12.0.txt", sep='\t', comment='#', names=['string_protein_id', 'preferred_name', 'protein_size', 'annotation'] ) ensp_to_name = info.set_index('string_protein_id')['preferred_name'].to_dict() # Load protein links - USE DETAILED FILE links = pd.read_csv( string_db_dir / "9606.protein.links.detailed.v12.0.txt", sep=' ' ) # Load enrichment terms # Note: STRING files usually don't have headers, or have specific comment lines. # We ensure we capture string_protein_id, category, term, description enrichment = pd.read_csv( string_db_dir / "9606.protein.enrichment.terms.v12.0.txt", sep='\t', comment='#', names=['string_protein_id', 'category', 'term', 'description'] ) print(f"Loaded {len(enst_to_ensp)} ENST->ENSP mappings") print(f"Loaded {len(links)} protein-protein interactions") print(f"Loaded {len(enrichment)} enrichment annotations") return enst_to_ensp, ensp_to_name, links, enrichment def perform_enrichment_analysis(query_proteins, enrichment_df, ensp_to_name, ensp_to_input_map): """ Perform enrichment analysis similar to STRING-DB API. Uses Fisher's exact test for each term. """ results = [] query_set = set(query_proteins) n_query = len(query_set) print(f"Analyzing enrichment for {n_query} proteins...") # Group by category to calculate background size per category for category, cat_group in enrichment_df.groupby('category'): # Background: Total unique proteins annotated in this category category_universe = set(cat_group['string_protein_id']) background_size = len(category_universe) # Now iterate over terms within this category for term, group in cat_group.groupby('term'): term_proteins = set(group['string_protein_id']) # Intersection: proteins in both query and term intersection = query_set & term_proteins n_intersection = len(intersection) if n_intersection == 0: continue # Background term count n_term_background = len(term_proteins) # Fisher's exact test # Contingency table: # In term Not in term # In query a b # Not in query c d a = n_intersection b = n_query - a c = n_term_background - a d = background_size - n_term_background - b # Ensure no negative values (safety check) if c < 0 or d < 0: continue # One-sided test for enrichment oddsratio, pvalue = stats.fisher_exact([[a, b], [c, d]], alternative='greater') # Get description (take first) description = group['description'].iloc[0] # Get input genes (map ENSP back to ENST input) intersection_list = sorted(list(intersection)) # Reconstruct the 'inputGenes' list using the original ENST IDs # If an ENSP maps to multiple ENSTs in the input, we list them all input_genes_list = [] for ensp in intersection_list: if ensp in ensp_to_input_map: # Append the ENST that mapped to this ENSP input_genes_list.append(ensp_to_input_map[ensp]) else: input_genes_list.append(ensp) input_genes = ','.join(input_genes_list) preferred_names = ','.join([ensp_to_name.get(p, p.split('.')[-1]) for p in intersection_list]) results.append({ 'category': category, 'term': term, 'number_of_genes': n_intersection, 'number_of_genes_in_background': n_term_background, 'ncbiTaxonId': 9606, 'inputGenes': input_genes, 'preferredNames': preferred_names, 'p_value': pvalue, 'description': description }) if not results: return pd.DataFrame() # Create DataFrame and calculate FDR results_df = pd.DataFrame(results) # Benjamini-Hochberg FDR correction results_df = results_df.sort_values('p_value') n_tests = len(results_df) results_df['fdr'] = results_df['p_value'] * n_tests / (np.arange(1, n_tests + 1)) results_df['fdr'] = results_df['fdr'].clip(upper=1.0) # Sort by p-value results_df = results_df.sort_values('p_value').reset_index(drop=True) # Reorder columns to match API output # API Order: category, term, number_of_genes, number_of_genes_in_background, ncbiTaxonId, inputGenes, preferredNames, p_value, fdr, description cols = [ 'category', 'term', 'number_of_genes', 'number_of_genes_in_background', 'ncbiTaxonId', 'inputGenes', 'preferredNames', 'p_value', 'fdr', 'description' ] results_df = results_df[cols] return results_df def get_network_interactions(query_proteins, links_df, ensp_to_name): """ Extract network interactions for query proteins. Converts combined_score to normalized scores similar to STRING API. """ query_set = set(query_proteins) # Filter links where both proteins are in query set interactions = links_df[ links_df['protein1'].isin(query_set) & links_df['protein2'].isin(query_set) ].copy() if len(interactions) == 0: return pd.DataFrame() # Normalize scores (STRING files are integers 0-1000, API is float 0-1) score_cols = ['combined_score', 'neighborhood', 'fusion', 'cooccurence', 'coexpression', 'experimental', 'database', 'textmining'] target_cols = ['score', 'nscore', 'fscore', 'pscore', 'ascore', 'escore', 'dscore', 'tscore'] for src, tgt in zip(score_cols, target_cols): if src in interactions.columns: interactions[tgt] = interactions[src] / 1000.0 else: interactions[tgt] = 0.0 # Add preferred names interactions['preferredName_A'] = interactions['protein1'].map( lambda x: ensp_to_name.get(x, x.split('.')[-1]) ) interactions['preferredName_B'] = interactions['protein2'].map( lambda x: ensp_to_name.get(x, x.split('.')[-1]) ) # Rename columns to match STRING API output interactions = interactions.rename(columns={ 'protein1': 'stringId_A', 'protein2': 'stringId_B' }) interactions['ncbiTaxonId'] = 9606 # Select and order columns to match original output column_order = [ 'stringId_A', 'stringId_B', 'preferredName_A', 'preferredName_B', 'ncbiTaxonId', 'score', 'nscore', 'fscore', 'pscore', 'ascore', 'escore', 'dscore', 'tscore' ] interactions = interactions[column_order].reset_index(drop=True) return interactions def main(): parser = argparse.ArgumentParser( description='Offline network enrichment analysis using local STRING-DB files' ) parser.add_argument( 'interactions_file', type=str, help='Input TSV file with significant interactions (conplex output)' ) parser.add_argument( '--string-db-dir', type=str, default='/data/bugra/digital_trials/app_network', help='Directory containing STRING-DB files' ) parser.add_argument( '--threshold', type=float, default=0.65, help='ConPlex score threshold for protein network' ) parser.add_argument( '--max-proteins', type=int, default=450, help='Maximum number of proteins to analyze' ) parser.add_argument( '--output-dir', type=str, default='.', help='Output directory for results' ) args = parser.parse_args() # Load STRING-DB data enst_to_ensp, ensp_to_name, links, enrichment = load_string_data(args.string_db_dir) # Load interaction data print(f"\nLoading interaction data from {args.interactions_file}...") interaction_data = pd.read_csv(args.interactions_file, sep='\t') interaction_data = interaction_data.sort_values('conplex_score', ascending=False) # Remove duplicates and limit # Replicates Nextflow logic: drops duplicates on transcript BEFORE limit interaction_data = interaction_data.drop_duplicates('transcipt') if interaction_data.shape[0] > args.max_proteins: print(f"Limiting to top {args.max_proteins} proteins") interaction_data = interaction_data.iloc[:args.max_proteins] # Filter by threshold filtered_data = interaction_data[interaction_data['conplex_score'] > args.threshold] print(f"Found {len(filtered_data)} proteins above threshold {args.threshold}") # Map ENST to ENSP and keep track of mapping for output enst_list = filtered_data['transcipt'].tolist() ensp_list = [] ensp_to_input_map = {} # Map ENSP back to ENST for output generation for enst in enst_list: ensp = enst_to_ensp.get(enst) if ensp: ensp_list.append(ensp) # We map ENSP back to the ENST. If multiple ENSTs map to one ENSP, # this simple dict keeps the last one. # However, since input is deduped on transcript, this is mostly 1-to-1 # for the query set. ensp_to_input_map[ensp] = enst print(f"Mapped {len(ensp_list)} ENST IDs to ENSP IDs") if len(ensp_list) == 0: print("ERROR: No valid ENSP mappings found!") return # Perform enrichment analysis print("\nPerforming enrichment analysis...") enrichment_results = perform_enrichment_analysis(ensp_list, enrichment, ensp_to_name, ensp_to_input_map) print(f"Found {len(enrichment_results)} enriched terms") # Get network interactions print("\nExtracting network interactions...") network_interactions = get_network_interactions(ensp_list, links, ensp_to_name) print(f"Found {len(network_interactions)} protein-protein interactions") # Generate output filename input_path = Path(args.interactions_file) output_name = input_path.stem.replace('_significant_interactions', '') output_dir = Path(args.output_dir) output_dir.mkdir(parents=True, exist_ok=True) # Save results enrichment_file = output_dir / f"{output_name}_network_enrichment.tsv" interactions_file = output_dir / f"{output_name}_network_interactions.tsv" enrichment_results.to_csv(enrichment_file, sep='\t', index=False) network_interactions.to_csv(interactions_file, sep='\t', index=False) print(f"\nResults saved:") print(f" Enrichment: {enrichment_file}") print(f" Interactions: {interactions_file}") if __name__ == '__main__': main()