plot_file10_fits.py

#!/usr/bin/env python3
"""
Plot the Lorentzian fits for File 10 mixture
"""

import sys
sys.path.insert(0, '.')
import os
import numpy as np
import matplotlib.pyplot as plt
from quantify_all_metabolites_v3_ref20 import (
    read_and_process, find_tsp_peak, integrate_peak,
    detect_peaks_in_region, fit_region
)

def plot_file10_fits():
    folder_path = "raw_data/Model_Mixtures-M1-M6_JCAMP-DX"
    filename = "10.dx"
    
    # Load spectrum
    ppm, spec = read_and_process(os.path.join(folder_path, filename))
    tsp = find_tsp_peak(ppm, spec)
    ppm_corr = ppm - tsp
    
    # Normalize
    tsp_area = integrate_peak(ppm_corr, spec, (-0.2, 0.2))
    spec_norm = spec / tsp_area
    
    # Key metabolites to visualize
    metabolites_to_plot = [
        ('Lactate CH3', (1.25, 1.40), [1.33]),
        ('Alanine CH3', (1.40, 1.55), [1.48]),
        ('Glutamate', (2.05, 2.50), [2.14, 2.42]),
        ('Glucose anomeric', (5.15, 5.35), [5.24]),
    ]
    
    fig, axes = plt.subplots(2, 2, figsize=(16, 12))
    axes = axes.flatten()
    
    for idx, (label, region, peaks) in enumerate(metabolites_to_plot):
        ax = axes[idx]
        
        # Extract data
        mask = (ppm_corr >= region[0]) & (ppm_corr <= region[1])
        x = ppm_corr[mask]
        y = spec_norm[mask]
        
        # Plot data
        ax.plot(x, y, 'b.', markersize=3, label='Data', alpha=0.6)
        
        # Try to fit
        try:
            detected = detect_peaks_in_region(ppm_corr, spec_norm, region, peaks, len(peaks))
            result, _ = fit_region(x, y, detected, len(peaks))
            
            # Plot fit
            ax.plot(x, result.best_fit, 'r-', linewidth=2, label='Fit')
            
            # Mark detected peaks
            for p in detected:
                ax.axvline(x=p, color='g', linestyle='--', alpha=0.5)
                ax.text(p, np.max(y)*0.9, f'{p:.2f}', ha='center', fontsize=9, color='green')
            
            # Get amplitude
            if len(peaks) == 1:
                amp = result.params['amplitude'].value
            else:
                amp = sum([result.params[f'p{i}_amplitude'].value for i in range(1, len(peaks)+1)])
            
            ax.set_title(f'{label}\nRegion: {region}, Amplitude: {amp:.3f}', fontsize=11, fontweight='bold')
            
        except Exception as e:
            ax.set_title(f'{label}\nError: {e}', fontsize=11)
        
        ax.set_xlabel('ppm', fontsize=10)
        ax.set_ylabel('Normalized Intensity', fontsize=10)
        ax.legend(fontsize=9)
        ax.grid(True, alpha=0.3)
    
    plt.suptitle(f'Model Mixture File 10 - Lorentzian Fits\nTSP Correction: +{-tsp:.4f} ppm', 
                 fontsize=14, fontweight='bold')
    plt.tight_layout()
    plt.savefig('file10_fits.png', dpi=150, bbox_inches='tight')
    print("Saved: file10_fits.png")
    plt.close()
    
    # Also plot full spectrum
    fig, ax = plt.subplots(figsize=(16, 8))
    mask = (ppm_corr >= 0.5) & (ppm_corr <= 5.5)
    ax.plot(ppm_corr[mask], spec_norm[mask], linewidth=0.8, color='blue')
    ax.set_xlim(5.5, 0.5)
    ax.set_xlabel('Chemical Shift (ppm)', fontsize=12)
    ax.set_ylabel('Normalized Intensity (/TSP)', fontsize=12)
    ax.set_title('Model Mixture File 10 - Full Spectrum (Normalized by TSP)', fontsize=14, fontweight='bold')
    ax.grid(True, alpha=0.3)
    
    # Annotate major peaks
    annotations = [
        (1.33, 'Lactate CH3'),
        (1.48, 'Alanine CH3'),
        (2.14, 'Glu β-CH2'),
        (2.42, 'Glu γ-CH2'),
        (2.75, 'Asn β-CH2'),
        (3.98, 'Asn α-CH'),
        (4.12, 'Lactate CH'),
        (5.24, 'Glc H-1'),
    ]
    
    for ppm_val, label in annotations:
        idx_closest = np.argmin(np.abs(ppm_corr[mask] - ppm_val))
        y_val = spec_norm[mask][idx_closest]
        if y_val > 0.5:  # Only annotate if peak is significant
            ax.annotate(label, xy=(ppm_val, y_val), xytext=(ppm_val+0.1, y_val+0.5),
                       arrowprops=dict(arrowstyle='->', color='red', alpha=0.5),
                       fontsize=9, color='red')
    
    plt.tight_layout()
    plt.savefig('file10_full_spectrum.png', dpi=150, bbox_inches='tight')
    print("Saved: file10_full_spectrum.png")
    plt.close()

if __name__ == "__main__":
    plot_file10_fits()