| 
492
 | 
     1 # Python API Tutorial
 | 
| 
 | 
     2 
 | 
| 
 | 
     3 Learn how to use COBRAxy tools programmatically in Python scripts and analysis pipelines.
 | 
| 
 | 
     4 
 | 
| 
 | 
     5 ## Overview
 | 
| 
 | 
     6 
 | 
| 
 | 
     7 This tutorial teaches you to integrate COBRAxy into Python workflows by calling tool main functions directly with parsed arguments.
 | 
| 
 | 
     8 
 | 
| 
 | 
     9 **Time required**: ~45 minutes  
 | 
| 
 | 
    10 **Difficulty**: Intermediate  
 | 
| 
 | 
    11 **Prerequisites**: Basic Python knowledge, COBRAxy installation
 | 
| 
 | 
    12 
 | 
| 
 | 
    13 ## Understanding COBRAxy Architecture
 | 
| 
 | 
    14 
 | 
| 
 | 
    15 ### Tool Structure
 | 
| 
 | 
    16 
 | 
| 
 | 
    17 Each COBRAxy tool is a Python module with:
 | 
| 
 | 
    18 - `main(args)` function that accepts argument list
 | 
| 
 | 
    19 - Command-line argument parsing
 | 
| 
 | 
    20 - Self-contained execution logic
 | 
| 
 | 
    21 
 | 
| 
 | 
    22 ```python
 | 
| 
 | 
    23 # General pattern for all tools
 | 
| 
 | 
    24 import tool_module
 | 
| 
 | 
    25 tool_module.main(['-arg1', 'value1', '-arg2', 'value2'])
 | 
| 
 | 
    26 ```
 | 
| 
 | 
    27 
 | 
| 
 | 
    28 ### Available Tools
 | 
| 
 | 
    29 
 | 
| 
 | 
    30 | Python Module | Purpose | Key Arguments |
 | 
| 
 | 
    31 |---------------|---------|---------------|
 | 
| 
 | 
    32 | `ras_generator` | Compute reaction activity scores | `-in`, `-ra`, `-rs` |
 | 
| 
 | 
    33 | `rps_generator` | Compute reaction propensity scores | `-id`, `-rp` |
 | 
| 
 | 
    34 | `marea` | Statistical pathway analysis | `-input_data`, `-choice_map` |
 | 
| 
 | 
    35 | `ras_to_bounds` | Apply RAS constraints to model | `-ir`, `-ms`, `-idop` |
 | 
| 
 | 
    36 | `flux_simulation` | Sample metabolic fluxes | `-ms`, `-in`, `-a`, `-ns` |
 | 
| 
 | 
    37 | `flux_to_map` | Add flux data to maps | `-if`, `-mp`, `-idop` |
 | 
| 
 | 
    38 
 | 
| 
 | 
    39 ## Setup Your Environment
 | 
| 
 | 
    40 
 | 
| 
 | 
    41 ### Import Required Modules
 | 
| 
 | 
    42 
 | 
| 
 | 
    43 ```python
 | 
| 
 | 
    44 import sys
 | 
| 
 | 
    45 import os
 | 
| 
 | 
    46 from pathlib import Path
 | 
| 
 | 
    47 
 | 
| 
 | 
    48 # Add COBRAxy to Python path
 | 
| 
 | 
    49 cobraxy_path = "/path/to/COBRAxy"
 | 
| 
 | 
    50 sys.path.insert(0, cobraxy_path)
 | 
| 
 | 
    51 
 | 
| 
 | 
    52 # Import COBRAxy tools
 | 
| 
 | 
    53 import ras_generator
 | 
| 
 | 
    54 import rps_generator
 | 
| 
 | 
    55 import marea
 | 
| 
 | 
    56 import ras_to_bounds
 | 
| 
 | 
    57 import flux_simulation  
 | 
| 
 | 
    58 import flux_to_map
 | 
| 
 | 
    59 import metabolicModel2Tabular as model_setting
 | 
| 
 | 
    60 ```
 | 
| 
 | 
    61 
 | 
| 
 | 
    62 ### Set Working Directory
 | 
| 
 | 
    63 
 | 
| 
 | 
    64 ```python
 | 
| 
 | 
    65 # Set up working directory
 | 
| 
 | 
    66 work_dir = Path("/path/to/analysis")
 | 
| 
 | 
    67 work_dir.mkdir(exist_ok=True)
 | 
| 
 | 
    68 os.chdir(work_dir)
 | 
| 
 | 
    69 
 | 
| 
 | 
    70 # COBRAxy tools expect this parameter
 | 
| 
 | 
    71 tool_dir = str(Path(cobraxy_path).absolute())
 | 
| 
 | 
    72 ```
 | 
| 
 | 
    73 
 | 
| 
 | 
    74 ## Basic Workflow Example
 | 
| 
 | 
    75 
 | 
| 
 | 
    76 ### Step 1: Prepare Sample Data
 | 
| 
 | 
    77 
 | 
| 
 | 
    78 ```python
 | 
| 
 | 
    79 import pandas as pd
 | 
| 
 | 
    80 import numpy as np
 | 
| 
 | 
    81 
 | 
| 
 | 
    82 # Create sample gene expression data
 | 
| 
 | 
    83 genes = ['HGNC:5', 'HGNC:10', 'HGNC:15', 'HGNC:25', 'HGNC:30']
 | 
| 
 | 
    84 samples = ['Control_1', 'Control_2', 'Treatment_1', 'Treatment_2']
 | 
| 
 | 
    85 
 | 
| 
 | 
    86 # Generate random expression values
 | 
| 
 | 
    87 np.random.seed(42)
 | 
| 
 | 
    88 data = np.random.lognormal(mean=2, sigma=1, size=(len(genes), len(samples)))
 | 
| 
 | 
    89 
 | 
| 
 | 
    90 # Create DataFrame
 | 
| 
 | 
    91 expression_df = pd.DataFrame(data, index=genes, columns=samples)
 | 
| 
 | 
    92 expression_df.index.name = 'Gene_ID'
 | 
| 
 | 
    93 
 | 
| 
 | 
    94 # Save to file
 | 
| 
 | 
    95 expression_file = work_dir / "expression_data.tsv"
 | 
| 
 | 
    96 expression_df.to_csv(expression_file, sep='\t')
 | 
| 
 | 
    97 print(f"Created sample data: {expression_file}")
 | 
| 
 | 
    98 ```
 | 
| 
 | 
    99 
 | 
| 
 | 
   100 ### Step 2: Extract Model Information
 | 
| 
 | 
   101 
 | 
| 
 | 
   102 ```python
 | 
| 
 | 
   103 # Extract model components (optional, for understanding model structure)
 | 
| 
 | 
   104 model_args = [
 | 
| 
 | 
   105     '-td', tool_dir,
 | 
| 
 | 
   106     '-ms', 'ENGRO2',  # Use built-in ENGRO2 model
 | 
| 
 | 
   107     '-idop', str(work_dir / 'model_info')
 | 
| 
 | 
   108 ]
 | 
| 
 | 
   109 
 | 
| 
 | 
   110 try:
 | 
| 
 | 
   111     model_setting.main(model_args)
 | 
| 
 | 
   112     print("✓ Model information extracted")
 | 
| 
 | 
   113 except Exception as e:
 | 
| 
 | 
   114     print(f"Model extraction failed: {e}")
 | 
| 
 | 
   115 ```
 | 
| 
 | 
   116 
 | 
| 
 | 
   117 ### Step 3: Generate RAS Scores
 | 
| 
 | 
   118 
 | 
| 
 | 
   119 ```python
 | 
| 
 | 
   120 # Generate Reaction Activity Scores
 | 
| 
 | 
   121 ras_output = work_dir / "ras_scores.tsv"
 | 
| 
 | 
   122 
 | 
| 
 | 
   123 ras_args = [
 | 
| 
 | 
   124     '-td', tool_dir,
 | 
| 
 | 
   125     '-in', str(expression_file),
 | 
| 
 | 
   126     '-ra', str(ras_output), 
 | 
| 
 | 
   127     '-rs', 'ENGRO2',  # Built-in model
 | 
| 
 | 
   128     '-n', 'true'  # Handle missing genes
 | 
| 
 | 
   129 ]
 | 
| 
 | 
   130 
 | 
| 
 | 
   131 try:
 | 
| 
 | 
   132     ras_generator.main(ras_args)
 | 
| 
 | 
   133     print(f"✓ RAS scores generated: {ras_output}")
 | 
| 
 | 
   134 except Exception as e:
 | 
| 
 | 
   135     print(f"RAS generation failed: {e}")
 | 
| 
 | 
   136     raise
 | 
| 
 | 
   137 ```
 | 
| 
 | 
   138 
 | 
| 
 | 
   139 ### Step 4: Generate RPS Scores (Optional)
 | 
| 
 | 
   140 
 | 
| 
 | 
   141 ```python
 | 
| 
 | 
   142 # Create sample metabolite data
 | 
| 
 | 
   143 metabolites = ['glucose', 'pyruvate', 'lactate', 'ATP', 'NADH']
 | 
| 
 | 
   144 met_data = np.random.lognormal(mean=3, sigma=0.5, size=(len(metabolites), len(samples)))
 | 
| 
 | 
   145 
 | 
| 
 | 
   146 met_df = pd.DataFrame(met_data, index=metabolites, columns=samples)
 | 
| 
 | 
   147 met_df.index.name = 'Metabolite_ID'
 | 
| 
 | 
   148 
 | 
| 
 | 
   149 metabolite_file = work_dir / "metabolite_data.tsv"
 | 
| 
 | 
   150 met_df.to_csv(metabolite_file, sep='\t')
 | 
| 
 | 
   151 
 | 
| 
 | 
   152 # Generate Reaction Propensity Scores
 | 
| 
 | 
   153 rps_output = work_dir / "rps_scores.tsv"
 | 
| 
 | 
   154 
 | 
| 
 | 
   155 rps_args = [
 | 
| 
 | 
   156     '-td', tool_dir,
 | 
| 
 | 
   157     '-id', str(metabolite_file),
 | 
| 
 | 
   158     '-rp', str(rps_output)
 | 
| 
 | 
   159 ]
 | 
| 
 | 
   160 
 | 
| 
 | 
   161 try:
 | 
| 
 | 
   162     rps_generator.main(rps_args)
 | 
| 
 | 
   163     print(f"✓ RPS scores generated: {rps_output}")
 | 
| 
 | 
   164 except Exception as e:
 | 
| 
 | 
   165     print(f"RPS generation warning: {e}")
 | 
| 
 | 
   166     # RPS generation might fail with sample data - that's OK
 | 
| 
 | 
   167 ```
 | 
| 
 | 
   168 
 | 
| 
 | 
   169 ### Step 5: Statistical Analysis with MAREA
 | 
| 
 | 
   170 
 | 
| 
 | 
   171 ```python
 | 
| 
 | 
   172 # Create enriched pathway maps
 | 
| 
 | 
   173 maps_output = work_dir / "pathway_maps"
 | 
| 
 | 
   174 
 | 
| 
 | 
   175 marea_args = [
 | 
| 
 | 
   176     '-td', tool_dir,
 | 
| 
 | 
   177     '-using_RAS', 'true',
 | 
| 
 | 
   178     '-input_data', str(ras_output),
 | 
| 
 | 
   179     '-choice_map', 'ENGRO2',
 | 
| 
 | 
   180     '-gs', 'true',  # Gene set analysis
 | 
| 
 | 
   181     '-idop', str(maps_output)
 | 
| 
 | 
   182 ]
 | 
| 
 | 
   183 
 | 
| 
 | 
   184 try:
 | 
| 
 | 
   185     marea.main(marea_args)
 | 
| 
 | 
   186     print(f"✓ Pathway maps created: {maps_output}")
 | 
| 
 | 
   187 except Exception as e:
 | 
| 
 | 
   188     print(f"MAREA analysis failed: {e}")
 | 
| 
 | 
   189 ```
 | 
| 
 | 
   190 
 | 
| 
 | 
   191 ### Step 6: Flux Simulation Pipeline
 | 
| 
 | 
   192 
 | 
| 
 | 
   193 ```python
 | 
| 
 | 
   194 # Apply RAS constraints to model
 | 
| 
 | 
   195 bounds_output = work_dir / "model_bounds"
 | 
| 
 | 
   196 
 | 
| 
 | 
   197 bounds_args = [
 | 
| 
 | 
   198     '-td', tool_dir,
 | 
| 
 | 
   199     '-ms', 'ENGRO2',
 | 
| 
 | 
   200     '-ir', str(ras_output),
 | 
| 
 | 
   201     '-rs', 'true',  # Use RAS values
 | 
| 
 | 
   202     '-idop', str(bounds_output)
 | 
| 
 | 
   203 ]
 | 
| 
 | 
   204 
 | 
| 
 | 
   205 try:
 | 
| 
 | 
   206     ras_to_bounds.main(bounds_args)
 | 
| 
 | 
   207     print(f"✓ Model constraints applied: {bounds_output}")
 | 
| 
 | 
   208 except Exception as e:
 | 
| 
 | 
   209     print(f"Bounds generation failed: {e}")
 | 
| 
 | 
   210     raise
 | 
| 
 | 
   211 
 | 
| 
 | 
   212 # Sample metabolic fluxes
 | 
| 
 | 
   213 flux_output = work_dir / "flux_samples"
 | 
| 
 | 
   214 
 | 
| 
 | 
   215 flux_args = [
 | 
| 
 | 
   216     '-td', tool_dir,
 | 
| 
 | 
   217     '-ms', 'ENGRO2',
 | 
| 
 | 
   218     '-in', str(bounds_output / "*.tsv"),  # Will be expanded by tool
 | 
| 
 | 
   219     '-a', 'CBS',  # Sampling algorithm
 | 
| 
 | 
   220     '-ns', '1000',  # Number of samples
 | 
| 
 | 
   221     '-idop', str(flux_output)
 | 
| 
 | 
   222 ]
 | 
| 
 | 
   223 
 | 
| 
 | 
   224 try:
 | 
| 
 | 
   225     flux_simulation.main(flux_args)
 | 
| 
 | 
   226     print(f"✓ Flux samples generated: {flux_output}")
 | 
| 
 | 
   227 except Exception as e:
 | 
| 
 | 
   228     print(f"Flux simulation failed: {e}")
 | 
| 
 | 
   229 ```
 | 
| 
 | 
   230 
 | 
| 
 | 
   231 ### Step 7: Create Final Visualizations
 | 
| 
 | 
   232 
 | 
| 
 | 
   233 ```python
 | 
| 
 | 
   234 # Add flux data to enriched maps
 | 
| 
 | 
   235 final_maps = work_dir / "final_visualizations"
 | 
| 
 | 
   236 
 | 
| 
 | 
   237 # Check if we have both maps and flux data
 | 
| 
 | 
   238 maps_dir = maps_output
 | 
| 
 | 
   239 flux_dir = flux_output
 | 
| 
 | 
   240 
 | 
| 
 | 
   241 if maps_dir.exists() and flux_dir.exists():
 | 
| 
 | 
   242     flux_to_map_args = [
 | 
| 
 | 
   243         '-td', tool_dir,
 | 
| 
 | 
   244         '-if', str(flux_dir / "*.tsv"),
 | 
| 
 | 
   245         '-mp', str(maps_dir / "*.svg"),
 | 
| 
 | 
   246         '-idop', str(final_maps)
 | 
| 
 | 
   247     ]
 | 
| 
 | 
   248     
 | 
| 
 | 
   249     try:
 | 
| 
 | 
   250         flux_to_map.main(flux_to_map_args)
 | 
| 
 | 
   251         print(f"✓ Final visualizations created: {final_maps}")
 | 
| 
 | 
   252     except Exception as e:
 | 
| 
 | 
   253         print(f"Final mapping failed: {e}")
 | 
| 
 | 
   254 else:
 | 
| 
 | 
   255     print("Skipping final visualization - missing input files")
 | 
| 
 | 
   256 ```
 | 
| 
 | 
   257 
 | 
| 
 | 
   258 ## Advanced Usage Patterns
 | 
| 
 | 
   259 
 | 
| 
 | 
   260 ### Error Handling and Validation
 | 
| 
 | 
   261 
 | 
| 
 | 
   262 ```python
 | 
| 
 | 
   263 def run_cobraxy_tool(tool_module, args, description):
 | 
| 
 | 
   264     """Helper function to run COBRAxy tools with error handling."""
 | 
| 
 | 
   265     try:
 | 
| 
 | 
   266         print(f"Running {description}...")
 | 
| 
 | 
   267         tool_module.main(args)
 | 
| 
 | 
   268         print(f"✓ {description} completed successfully")
 | 
| 
 | 
   269         return True
 | 
| 
 | 
   270     except Exception as e:
 | 
| 
 | 
   271         print(f"✗ {description} failed: {e}")
 | 
| 
 | 
   272         return False
 | 
| 
 | 
   273 
 | 
| 
 | 
   274 # Usage
 | 
| 
 | 
   275 success = run_cobraxy_tool(
 | 
| 
 | 
   276     ras_generator, 
 | 
| 
 | 
   277     ras_args, 
 | 
| 
 | 
   278     "RAS generation"
 | 
| 
 | 
   279 )
 | 
| 
 | 
   280 
 | 
| 
 | 
   281 if not success:
 | 
| 
 | 
   282     print("Pipeline stopped due to error")
 | 
| 
 | 
   283     exit(1)
 | 
| 
 | 
   284 ```
 | 
| 
 | 
   285 
 | 
| 
 | 
   286 ### Batch Processing Multiple Datasets
 | 
| 
 | 
   287 
 | 
| 
 | 
   288 ```python
 | 
| 
 | 
   289 def process_dataset(dataset_path, output_dir):
 | 
| 
 | 
   290     """Process a single dataset through COBRAxy pipeline."""
 | 
| 
 | 
   291     
 | 
| 
 | 
   292     dataset_name = dataset_path.stem
 | 
| 
 | 
   293     out_dir = Path(output_dir) / dataset_name
 | 
| 
 | 
   294     out_dir.mkdir(exist_ok=True)
 | 
| 
 | 
   295     
 | 
| 
 | 
   296     # Generate RAS
 | 
| 
 | 
   297     ras_file = out_dir / "ras_scores.tsv"
 | 
| 
 | 
   298     ras_args = [
 | 
| 
 | 
   299         '-td', tool_dir,
 | 
| 
 | 
   300         '-in', str(dataset_path),
 | 
| 
 | 
   301         '-ra', str(ras_file),
 | 
| 
 | 
   302         '-rs', 'ENGRO2'
 | 
| 
 | 
   303     ]
 | 
| 
 | 
   304     
 | 
| 
 | 
   305     if run_cobraxy_tool(ras_generator, ras_args, f"RAS for {dataset_name}"):
 | 
| 
 | 
   306         # Continue with MAREA analysis
 | 
| 
 | 
   307         maps_dir = out_dir / "maps"
 | 
| 
 | 
   308         marea_args = [
 | 
| 
 | 
   309             '-td', tool_dir,
 | 
| 
 | 
   310             '-using_RAS', 'true',
 | 
| 
 | 
   311             '-input_data', str(ras_file),
 | 
| 
 | 
   312             '-choice_map', 'ENGRO2',
 | 
| 
 | 
   313             '-idop', str(maps_dir)
 | 
| 
 | 
   314         ]
 | 
| 
 | 
   315         run_cobraxy_tool(marea, marea_args, f"MAREA for {dataset_name}")
 | 
| 
 | 
   316     
 | 
| 
 | 
   317     return out_dir
 | 
| 
 | 
   318 
 | 
| 
 | 
   319 # Process multiple datasets
 | 
| 
 | 
   320 datasets = [
 | 
| 
 | 
   321     "/path/to/dataset1.tsv",
 | 
| 
 | 
   322     "/path/to/dataset2.tsv", 
 | 
| 
 | 
   323     "/path/to/dataset3.tsv"
 | 
| 
 | 
   324 ]
 | 
| 
 | 
   325 
 | 
| 
 | 
   326 results = []
 | 
| 
 | 
   327 for dataset in datasets:
 | 
| 
 | 
   328     result_dir = process_dataset(Path(dataset), work_dir / "batch_results")
 | 
| 
 | 
   329     results.append(result_dir)
 | 
| 
 | 
   330     
 | 
| 
 | 
   331 print(f"Processed {len(results)} datasets")
 | 
| 
 | 
   332 ```
 | 
| 
 | 
   333 
 | 
| 
 | 
   334 ### Custom Analysis Pipelines
 | 
| 
 | 
   335 
 | 
| 
 | 
   336 ```python
 | 
| 
 | 
   337 class COBRAxyPipeline:
 | 
| 
 | 
   338     """Custom COBRAxy analysis pipeline."""
 | 
| 
 | 
   339     
 | 
| 
 | 
   340     def __init__(self, tool_dir, work_dir):
 | 
| 
 | 
   341         self.tool_dir = tool_dir
 | 
| 
 | 
   342         self.work_dir = Path(work_dir)
 | 
| 
 | 
   343         self.work_dir.mkdir(exist_ok=True)
 | 
| 
 | 
   344         
 | 
| 
 | 
   345     def run_enrichment_analysis(self, expression_file, model='ENGRO2'):
 | 
| 
 | 
   346         """Run enrichment-focused analysis."""
 | 
| 
 | 
   347         
 | 
| 
 | 
   348         # Generate RAS
 | 
| 
 | 
   349         ras_file = self.work_dir / "ras_scores.tsv"
 | 
| 
 | 
   350         ras_args = ['-td', self.tool_dir, '-in', str(expression_file), 
 | 
| 
 | 
   351                    '-ra', str(ras_file), '-rs', model]
 | 
| 
 | 
   352         
 | 
| 
 | 
   353         if not run_cobraxy_tool(ras_generator, ras_args, "RAS generation"):
 | 
| 
 | 
   354             return None
 | 
| 
 | 
   355             
 | 
| 
 | 
   356         # Run MAREA
 | 
| 
 | 
   357         maps_dir = self.work_dir / "enrichment_maps"  
 | 
| 
 | 
   358         marea_args = ['-td', self.tool_dir, '-using_RAS', 'true',
 | 
| 
 | 
   359                      '-input_data', str(ras_file), '-choice_map', model,
 | 
| 
 | 
   360                      '-gs', 'true', '-idop', str(maps_dir)]
 | 
| 
 | 
   361                      
 | 
| 
 | 
   362         if run_cobraxy_tool(marea, marea_args, "MAREA analysis"):
 | 
| 
 | 
   363             return maps_dir
 | 
| 
 | 
   364         return None
 | 
| 
 | 
   365     
 | 
| 
 | 
   366     def run_flux_analysis(self, expression_file, model='ENGRO2', n_samples=1000):
 | 
| 
 | 
   367         """Run flux sampling analysis."""
 | 
| 
 | 
   368         
 | 
| 
 | 
   369         # Generate RAS and apply bounds
 | 
| 
 | 
   370         ras_file = self.work_dir / "ras_scores.tsv" 
 | 
| 
 | 
   371         bounds_dir = self.work_dir / "bounds"
 | 
| 
 | 
   372         flux_dir = self.work_dir / "flux_samples"
 | 
| 
 | 
   373         
 | 
| 
 | 
   374         # RAS generation
 | 
| 
 | 
   375         ras_args = ['-td', self.tool_dir, '-in', str(expression_file),
 | 
| 
 | 
   376                    '-ra', str(ras_file), '-rs', model]
 | 
| 
 | 
   377         if not run_cobraxy_tool(ras_generator, ras_args, "RAS generation"):
 | 
| 
 | 
   378             return None
 | 
| 
 | 
   379             
 | 
| 
 | 
   380         # Apply bounds
 | 
| 
 | 
   381         bounds_args = ['-td', self.tool_dir, '-ms', model, '-ir', str(ras_file),
 | 
| 
 | 
   382                       '-rs', 'true', '-idop', str(bounds_dir)]
 | 
| 
 | 
   383         if not run_cobraxy_tool(ras_to_bounds, bounds_args, "Bounds application"):
 | 
| 
 | 
   384             return None
 | 
| 
 | 
   385             
 | 
| 
 | 
   386         # Flux sampling
 | 
| 
 | 
   387         flux_args = ['-td', self.tool_dir, '-ms', model, 
 | 
| 
 | 
   388                     '-in', str(bounds_dir / "*.tsv"),
 | 
| 
 | 
   389                     '-a', 'CBS', '-ns', str(n_samples), '-idop', str(flux_dir)]
 | 
| 
 | 
   390         
 | 
| 
 | 
   391         if run_cobraxy_tool(flux_simulation, flux_args, "Flux simulation"):
 | 
| 
 | 
   392             return flux_dir
 | 
| 
 | 
   393         return None
 | 
| 
 | 
   394 
 | 
| 
 | 
   395 # Usage
 | 
| 
 | 
   396 pipeline = COBRAxyPipeline(tool_dir, work_dir / "custom_analysis")
 | 
| 
 | 
   397 
 | 
| 
 | 
   398 # Run enrichment analysis
 | 
| 
 | 
   399 enrichment_results = pipeline.run_enrichment_analysis(expression_file)
 | 
| 
 | 
   400 if enrichment_results:
 | 
| 
 | 
   401     print(f"Enrichment analysis completed: {enrichment_results}")
 | 
| 
 | 
   402 
 | 
| 
 | 
   403 # Run flux analysis  
 | 
| 
 | 
   404 flux_results = pipeline.run_flux_analysis(expression_file, n_samples=500)
 | 
| 
 | 
   405 if flux_results:
 | 
| 
 | 
   406     print(f"Flux analysis completed: {flux_results}")
 | 
| 
 | 
   407 ```
 | 
| 
 | 
   408 
 | 
| 
 | 
   409 ## Integration with Data Analysis Libraries
 | 
| 
 | 
   410 
 | 
| 
 | 
   411 ### Pandas Integration
 | 
| 
 | 
   412 
 | 
| 
 | 
   413 ```python
 | 
| 
 | 
   414 # Read COBRAxy results back into pandas
 | 
| 
 | 
   415 ras_df = pd.read_csv(ras_output, sep='\t', index_col=0)
 | 
| 
 | 
   416 print(f"RAS data shape: {ras_df.shape}")
 | 
| 
 | 
   417 print(f"Sample statistics:\n{ras_df.describe()}")
 | 
| 
 | 
   418 
 | 
| 
 | 
   419 # Filter highly variable reactions
 | 
| 
 | 
   420 ras_std = ras_df.std(axis=1)
 | 
| 
 | 
   421 variable_reactions = ras_std.nlargest(20).index
 | 
| 
 | 
   422 print(f"Most variable reactions: {list(variable_reactions)}")
 | 
| 
 | 
   423 ```
 | 
| 
 | 
   424 
 | 
| 
 | 
   425 ### Matplotlib Visualization
 | 
| 
 | 
   426 
 | 
| 
 | 
   427 ```python
 | 
| 
 | 
   428 import matplotlib.pyplot as plt
 | 
| 
 | 
   429 import seaborn as sns
 | 
| 
 | 
   430 
 | 
| 
 | 
   431 # Visualize RAS distributions
 | 
| 
 | 
   432 plt.figure(figsize=(12, 8))
 | 
| 
 | 
   433 sns.heatmap(ras_df.iloc[:50], cmap='RdBu_r', center=0, cbar_kws={'label': 'RAS Score'})
 | 
| 
 | 
   434 plt.title('Reaction Activity Scores (Top 50 Reactions)')
 | 
| 
 | 
   435 plt.xlabel('Samples')
 | 
| 
 | 
   436 plt.ylabel('Reactions')
 | 
| 
 | 
   437 plt.tight_layout()
 | 
| 
 | 
   438 plt.savefig(work_dir / 'ras_heatmap.png', dpi=300)
 | 
| 
 | 
   439 plt.show()
 | 
| 
 | 
   440 ```
 | 
| 
 | 
   441 
 | 
| 
 | 
   442 ## Best Practices
 | 
| 
 | 
   443 
 | 
| 
 | 
   444 ### 1. Environment Management
 | 
| 
 | 
   445 ```python
 | 
| 
 | 
   446 # Use pathlib for cross-platform compatibility
 | 
| 
 | 
   447 from pathlib import Path
 | 
| 
 | 
   448 
 | 
| 
 | 
   449 # Use absolute paths
 | 
| 
 | 
   450 tool_dir = str(Path(cobraxy_path).absolute())
 | 
| 
 | 
   451 work_dir = Path("/analysis").absolute()
 | 
| 
 | 
   452 ```
 | 
| 
 | 
   453 
 | 
| 
 | 
   454 ### 2. Error Handling
 | 
| 
 | 
   455 ```python
 | 
| 
 | 
   456 # Always wrap tool calls in try-except
 | 
| 
 | 
   457 try:
 | 
| 
 | 
   458     ras_generator.main(ras_args)
 | 
| 
 | 
   459 except Exception as e:
 | 
| 
 | 
   460     print(f"RAS generation failed: {e}")
 | 
| 
 | 
   461     # Log details, cleanup, or alternative action
 | 
| 
 | 
   462 ```
 | 
| 
 | 
   463 
 | 
| 
 | 
   464 ### 3. Argument Validation
 | 
| 
 | 
   465 ```python
 | 
| 
 | 
   466 def validate_file_exists(filepath):
 | 
| 
 | 
   467     """Validate input file exists."""
 | 
| 
 | 
   468     path = Path(filepath)
 | 
| 
 | 
   469     if not path.exists():
 | 
| 
 | 
   470         raise FileNotFoundError(f"Input file not found: {filepath}")
 | 
| 
 | 
   471     return str(path.absolute())
 | 
| 
 | 
   472 
 | 
| 
 | 
   473 # Use before calling tools
 | 
| 
 | 
   474 expression_file = validate_file_exists(expression_file)
 | 
| 
 | 
   475 ```
 | 
| 
 | 
   476 
 | 
| 
 | 
   477 
 | 
| 
 | 
   478 
 | 
| 
 | 
   479 ## Troubleshooting
 | 
| 
 | 
   480 
 | 
| 
 | 
   481 ### Common Issues
 | 
| 
 | 
   482 
 | 
| 
 | 
   483 **Import errors**
 | 
| 
 | 
   484 ```python
 | 
| 
 | 
   485 # Check if COBRAxy path is correct
 | 
| 
 | 
   486 import sys
 | 
| 
 | 
   487 print("Python path includes:")
 | 
| 
 | 
   488 for p in sys.path:
 | 
| 
 | 
   489     print(f"  {p}")
 | 
| 
 | 
   490     
 | 
| 
 | 
   491 # Add COBRAxy path
 | 
| 
 | 
   492 sys.path.insert(0, "/correct/path/to/COBRAxy")
 | 
| 
 | 
   493 ```
 | 
| 
 | 
   494 
 | 
| 
 | 
   495 **Tool execution failures**
 | 
| 
 | 
   496 ```python
 | 
| 
 | 
   497 # Enable verbose output
 | 
| 
 | 
   498 import logging
 | 
| 
 | 
   499 logging.basicConfig(level=logging.DEBUG)
 | 
| 
 | 
   500 
 | 
| 
 | 
   501 # Check working directory
 | 
| 
 | 
   502 print(f"Current directory: {os.getcwd()}")
 | 
| 
 | 
   503 print(f"Directory contents: {list(Path('.').iterdir())}")
 | 
| 
 | 
   504 ```
 | 
| 
 | 
   505 
 | 
| 
 | 
   506 **File path issues**
 | 
| 
 | 
   507 ```python
 | 
| 
 | 
   508 # Use absolute paths
 | 
| 
 | 
   509 ras_args = [
 | 
| 
 | 
   510     '-td', str(Path(tool_dir).absolute()),
 | 
| 
 | 
   511     '-in', str(Path(expression_file).absolute()),
 | 
| 
 | 
   512     '-ra', str(Path(ras_output).absolute()),
 | 
| 
 | 
   513     '-rs', 'ENGRO2'
 | 
| 
 | 
   514 ]
 | 
| 
 | 
   515 ```
 | 
| 
 | 
   516 
 | 
| 
 | 
   517 ## Next Steps
 | 
| 
 | 
   518 
 | 
| 
 | 
   519 Now that you can use COBRAxy programmatically:
 | 
| 
 | 
   520 
 | 
| 
 | 
   521 1. **[Tools Reference](../tools/)** - Detailed parameter documentation  
 | 
| 
 | 
   522 2. **[Examples](../examples/)** - Real-world analysis scripts
 | 
| 
 | 
   523 3. **Build custom analysis pipelines** for your research needs
 | 
| 
 | 
   524 4. **Integrate with workflow managers** like Snakemake or Nextflow
 | 
| 
 | 
   525 
 | 
| 
 | 
   526 ## Resources
 | 
| 
 | 
   527 
 | 
| 
 | 
   528 - [COBRApy Documentation](https://cobrapy.readthedocs.io/) - Underlying metabolic modeling library
 | 
| 
 | 
   529 - [Pandas Documentation](https://pandas.pydata.org/) - Data manipulation
 | 
| 
 | 
   530 - [Matplotlib Gallery](https://matplotlib.org/gallery/) - Visualization examples
 | 
| 
 | 
   531 - [Python Pathlib](https://docs.python.org/3/library/pathlib.html) - Modern path handling |