run_integration_all_folders.py

#!/usr/bin/env python3
"""
Batch NMR Integration Script
============================

Runs peak integration on all amino acid folders in Reference_Raw_Date_JCAMP-DX.

Author: Assistant
Date: 2026-03-22
"""

import sys
import importlib.util
from pathlib import Path
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np

# Load the main integration script
spec = importlib.util.spec_from_file_location("nmr_integration", 
    "/home/tjiang/elise/pure_shift_nmr/new/nmr_integration_script.py")
nmr_module = importlib.util.module_from_spec(spec)
spec.loader.exec_module(nmr_module)

# Import functions from the module
get_even_numbered_files = nmr_module.get_even_numbered_files
read_nmr_data = nmr_module.read_nmr_data
integrate_peak = nmr_module.integrate_peak
process_all_files = nmr_module.process_all_files
plot_integration_results = nmr_module.plot_integration_results
plot_spectra_overlay = nmr_module.plot_spectra_overlay
print_summary_table = nmr_module.print_summary_table


def process_single_folder(folder_path, output_dir, ch3_region=(1.1, 1.8), ch_region=(3.5, 3.9)):
    """
    Process a single amino acid folder.
    
    Parameters:
    -----------
    folder_path : Path
        Path to the amino acid folder
    output_dir : Path
        Directory to save output files
    ch3_region : tuple
        CH3 integration region (ppm)
    ch_region : tuple
        CH integration region (ppm)
    
    Returns:
    --------
    dict : Summary results for this folder
    """
    folder_name = folder_path.name
    print(f"\n{'='*70}")
    print(f"PROCESSING: {folder_name}")
    print(f"{'='*70}")
    
    # Check if there are any .dx files
    dx_files = list(folder_path.glob("*.dx"))
    if not dx_files:
        print(f"  ⚠ No .dx files found in {folder_name}")
        return None
    
    # Process files
    results = process_all_files(folder_path, ch3_region, ch_region)
    
    if not results:
        print(f"  ⚠ No results for {folder_name}")
        return None
    
    # Generate output filenames
    base_name = folder_name.replace('-Reference', '').lower()
    integration_plot = output_dir / f"{base_name}_integration.png"
    overlay_plot = output_dir / f"{base_name}_spectra.png"
    
    # Generate plots
    plot_integration_results(results, output_file=str(integration_plot))
    plot_spectra_overlay(folder_path, output_file=str(overlay_plot), 
                        region=(0.5, 5.0))
    
    # Calculate average ratio
    avg_ratio_orig = np.mean([r['ratio_orig'] for r in results])
    avg_ratio_smooth = np.mean([r['ratio_smooth'] for r in results])
    
    # Calculate average areas
    avg_ch3 = np.mean([r['ch3_smoothed'] for r in results])
    avg_ch = np.mean([r['ch_smoothed'] for r in results])
    
    print(f"\n  ✅ {folder_name} complete!")
    print(f"     Files processed: {len(results)}")
    print(f"     Avg CH3/CH ratio: {avg_ratio_smooth:.2f}")
    
    return {
        'folder': folder_name,
        'num_files': len(results),
        'avg_ratio': avg_ratio_smooth,
        'avg_ch3': avg_ch3,
        'avg_ch': avg_ch,
        'results': results
    }


def create_summary_comparison(all_summaries, output_dir):
    """
    Create a summary comparison plot across all amino acids.
    
    Parameters:
    -----------
    all_summaries : list
        List of summary dictionaries from each folder
    output_dir : Path
        Directory to save output file
    """
    if not all_summaries:
        print("No summaries to plot!")
        return
    
    # Sort by folder name
    all_summaries = sorted(all_summaries, key=lambda x: x['folder'])
    
    names = [s['folder'].replace('-Reference', '') for s in all_summaries]
    ratios = [s['avg_ratio'] for s in all_summaries]
    
    fig, axes = plt.subplots(2, 1, figsize=(14, 10))
    
    # Plot 1: CH3/CH ratios
    ax1 = axes[0]
    colors = plt.cm.tab20(np.linspace(0, 1, len(names)))
    bars = ax1.bar(range(len(names)), ratios, color=colors, edgecolor='black')
    ax1.set_xlabel('Amino Acid', fontsize=12)
    ax1.set_ylabel('CH₃/CH Area Ratio', fontsize=12)
    ax1.set_title('Average CH₃/CH Peak Area Ratio by Amino Acid\n(Smoothed, σ=3)', 
                  fontsize=14, fontweight='bold')
    ax1.set_xticks(range(len(names)))
    ax1.set_xticklabels(names, rotation=45, ha='right')
    ax1.axhline(y=3, color='r', linestyle='--', alpha=0.5, 
                label='Theoretical (3:1)')
    ax1.legend()
    ax1.grid(True, alpha=0.3, axis='y')
    
    # Add values on bars
    for i, (bar, ratio) in enumerate(zip(bars, ratios)):
        height = bar.get_height()
        ax1.text(bar.get_x() + bar.get_width()/2., height,
                f'{ratio:.2f}', ha='center', va='bottom', fontsize=9)
    
    # Plot 2: Number of files processed
    ax2 = axes[1]
    num_files = [s['num_files'] for s in all_summaries]
    bars2 = ax2.bar(range(len(names)), num_files, color='steelblue', 
                    edgecolor='black', alpha=0.7)
    ax2.set_xlabel('Amino Acid', fontsize=12)
    ax2.set_ylabel('Number of Files', fontsize=12)
    ax2.set_title('Number of Even-numbered Files Processed', 
                  fontsize=14, fontweight='bold')
    ax2.set_xticks(range(len(names)))
    ax2.set_xticklabels(names, rotation=45, ha='right')
    ax2.grid(True, alpha=0.3, axis='y')
    
    # Add values on bars
    for bar in bars2:
        height = bar.get_height()
        ax2.text(bar.get_x() + bar.get_width()/2., height,
                f'{int(height)}', ha='center', va='bottom', fontsize=10)
    
    plt.tight_layout()
    output_file = output_dir / 'summary_all_amino_acids.png'
    plt.savefig(output_file, dpi=200, bbox_inches='tight')
    print(f"\n✅ Summary comparison saved to: {output_file}")
    plt.close()
    
    # Print summary table
    print("\n" + "=" * 70)
    print("SUMMARY TABLE - ALL AMINO ACIDS")
    print("=" * 70)
    print(f"\n{'Amino Acid':<20} {'Files':<8} {'Avg CH3/CH Ratio':<20} {'Status':<15}")
    print("-" * 70)
    
    for s in all_summaries:
        name = s['folder'].replace('-Reference', '')
        ratio = s['avg_ratio']
        num = s['num_files']
        # Check if ratio is close to 3.0 (±20%)
        if 2.4 <= ratio <= 3.6:
            status = "✅ Good"
        elif 1.5 <= ratio <= 4.5:
            status = "⚠️ Fair"
        else:
            status = "❌ Poor"
        print(f"{name:<20} {num:<8} {ratio:<20.2f} {status:<15}")


def main():
    """
    Main function to process all amino acid folders.
    """
    # Base directory
    BASE_DIR = Path("/home/tjiang/elise/pure_shift_nmr/new/raw_data/"
                    "Reference_Raw_Date_JCAMP-DX")
    
    # Output directory
    OUTPUT_DIR = Path("/home/tjiang/elise/pure_shift_nmr/new/integration_results")
    OUTPUT_DIR.mkdir(exist_ok=True)
    
    # Integration regions
    CH3_REGION = (1.1, 1.8)
    CH_REGION = (3.5, 3.9)
    
    print("=" * 70)
    print("BATCH NMR INTEGRATION - ALL AMINO ACIDS")
    print("=" * 70)
    print(f"\nBase directory: {BASE_DIR}")
    print(f"Output directory: {OUTPUT_DIR}")
    print(f"CH3 region: {CH3_REGION[0]}-{CH3_REGION[1]} ppm")
    print(f"CH region: {CH_REGION[0]}-{CH_REGION[1]} ppm")
    
    # Find all folders
    folders = sorted([d for d in BASE_DIR.iterdir() if d.is_dir()])
    
    print(f"\nFound {len(folders)} folders to process:")
    for f in folders:
        print(f"  - {f.name}")
    
    # Process each folder
    all_summaries = []
    
    for folder in folders:
        summary = process_single_folder(folder, OUTPUT_DIR, 
                                       CH3_REGION, CH_REGION)
        if summary:
            all_summaries.append(summary)
    
    # Create overall summary
    print("\n" + "=" * 70)
    print("CREATING SUMMARY COMPARISON")
    print("=" * 70)
    create_summary_comparison(all_summaries, OUTPUT_DIR)
    
    print("\n" + "=" * 70)
    print("BATCH PROCESSING COMPLETE!")
    print("=" * 70)
    print(f"\nResults saved to: {OUTPUT_DIR}")
    print(f"Total amino acids processed: {len(all_summaries)}")


if __name__ == "__main__":
    main()