plot_all_pure_shift.py

#!/usr/bin/env python3
"""
Plot all 18 pure shift spectra in one figure
"""

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

def plot_all_pure_shift():
    folder_path = "raw_data/Model_Mixtures-M1-M6_JCAMP-DX"
    
    # Get all files
    files = sorted([f for f in os.listdir(folder_path) if f.endswith('.dx')],
                   key=lambda x: int(x.replace('.dx', '')))
    
    # Get pure shift files (odd indices: 1, 3, 5, ...)
    pure_shift_files = [files[i] for i in range(1, len(files), 2)]
    
    print(f"Plotting {len(pure_shift_files)} pure shift files...")
    
    # Create figure with 6x3 subplots
    fig, axes = plt.subplots(6, 3, figsize=(18, 24))
    axes = axes.flatten()
    
    for idx, filename in enumerate(pure_shift_files):
        ax = axes[idx]
        
        try:
            ppm, spec = read_and_process(os.path.join(folder_path, filename))
            tsp = find_tsp_peak(ppm, spec)
            ppm_corr = ppm - tsp
            
            # Plot full spectrum (0.5-5.5 ppm)
            mask = (ppm_corr >= 0.5) & (ppm_corr <= 5.5)
            ax.plot(ppm_corr[mask], spec[mask], linewidth=0.8, color='blue')
            
            # Styling
            ax.set_xlim(5.5, 0.5)
            ax.set_title(f'{filename.replace(".dx", "")}', fontsize=12, fontweight='bold')
            ax.set_xlabel('ppm', fontsize=9)
            ax.set_ylabel('Intensity', fontsize=9)
            ax.tick_params(labelsize=8)
            ax.grid(True, alpha=0.2)
            
        except Exception as e:
            ax.text(0.5, 0.5, f'Error: {e}', ha='center', va='center', transform=ax.transAxes)
            ax.set_title(f'{filename}', fontsize=10)
    
    plt.suptitle('All 18 Pure Shift Mixture Spectra (TSP-corrected)', 
                 fontsize=16, fontweight='bold', y=1.00)
    plt.tight_layout()
    
    output_path = 'all_pure_shift_spectra.png'
    plt.savefig(output_path, dpi=150, bbox_inches='tight')
    print(f"Saved to: {output_path}")
    plt.close()

if __name__ == "__main__":
    plot_all_pure_shift()