alanine_quantification_analysis.py

#!/usr/bin/env python3
"""
Alanine Quantification Analysis
===============================

Analyze the linear relationship between TSP-scaled peak area and concentration
for Alanine reference samples.

The procedure follows the paper's method:
1. Read pure shift spectra
2. Find TSP peak and integrate TSP
3. Integrate Alanine CH3 peak (1.48 ppm region)
4. Calculate TSP-scaled peak area
5. Compare with known concentrations from Excel file
"""

import nmrglue as ng
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d
from pathlib import Path
from scipy import stats


def read_nmr_data(file_path):
    """Read JCAMP-DX file and return ppm axis and magnitude spectrum."""
    dic, data_list = ng.jcampdx.read(str(file_path))
    
    # Get real and imaginary parts
    data_real = data_list[0]
    data_imag = data_list[1] if len(data_list) > 1 else np.zeros_like(data_real)
    
    # Calculate magnitude spectrum
    data_magnitude = np.sqrt(data_real**2 + data_imag**2)
    
    # Calculate ppm axis
    sfo1 = float(dic['$SFO1'][0])
    o1_hz = float(dic['$O1'][0])
    sw_hz = float(dic['$SWH'][0])
    si = data_real.size
    
    o1_ppm = o1_hz / sfo1
    sw_ppm = sw_hz / sfo1
    ppm = np.linspace(o1_ppm + sw_ppm/2, o1_ppm - sw_ppm/2, si)
    
    return ppm, data_magnitude, sfo1


def find_tsp_peak(ppm, data, tsp_region=(-0.5, 0.5)):
    """Find the TSP peak position."""
    mask = (ppm >= tsp_region[0]) & (ppm <= tsp_region[1])
    if not np.any(mask):
        return None
    
    ppm_region = ppm[mask]
    data_region = data[mask]
    
    peak_idx = np.argmax(data_region)
    peak_ppm = ppm_region[peak_idx]
    peak_intensity = data_region[peak_idx]
    
    # S/N check
    noise = np.std(data_region[data_region < np.percentile(data_region, 50)])
    sn_ratio = peak_intensity / noise if noise > 0 else 0
    
    if sn_ratio < 10:
        return None
    
    return peak_ppm


def integrate_peak(ppm, intensity, region):
    """Integrate peak area within specified ppm region."""
    mask = (ppm >= region[0]) & (ppm <= region[1])
    if not np.any(mask):
        return 0
    
    ppm_region = ppm[mask]
    intensity_region = intensity[mask]
    
    # Use absolute value because ppm axis goes high->low
    area = abs(np.trapz(intensity_region, ppm_region))
    return area


def get_alanine_concentrations_from_excel():
    """Extract Alanine concentration data from Excel file."""
    df = pd.read_excel('2_Analytical Chemistry Data.xlsx', 
                       sheet_name='Reference Library Data Analysis', header=None)
    
    # Extract Alanine data (rows 27-33)
    ala_data = []
    for idx in range(28, 34):  # Rows with actual data
        row = df.iloc[idx]
        file_name = row[1]  # e.g., 'Xi-Reference-Ala-20'
        concentration = row[3]  # Concentration in mM
        
        # Extract the file number (e.g., '20' from 'Xi-Reference-Ala-20')
        if isinstance(file_name, str) and '-' in file_name:
            parts = file_name.split('-')
            if len(parts) >= 3 and parts[-1].isdigit():
                file_num = int(parts[-1])
                
                # Skip the 40mM reference (will use as reference, not sample)
                if file_num != 10:
                    ala_data.append({
                        'file_num': file_num,
                        'file_name': file_name,
                        'concentration_mM': float(concentration)
                    })
    
    return ala_data


def process_alanine_files():
    """Process Alanine reference files and calculate TSP-scaled integrations."""
    
    # Get concentrations from Excel
    ala_data = get_alanine_concentrations_from_excel()
    print("Alanine concentration data from Excel:")
    for item in ala_data:
        print(f"  File {item['file_num']}: {item['concentration_mM']:.3f} mM")
    print()
    
    # Base directory
    base_dir = Path("raw_data/Reference_Raw_Date_JCAMP-DX/Alanine-Reference")
    
    # Use 40mM (file 10) as the reference sample
    ref_file = base_dir / "10.dx"
    
    print("Processing reference sample (40 mM, file 10)...")
    ref_ppm, ref_data, ref_sfo1 = read_nmr_data(ref_file)
    
    # Apply smoothing
    ref_smooth = gaussian_filter1d(ref_data, sigma=2)
    
    # Find TSP and apply correction
    ref_tsp_ppm = find_tsp_peak(ref_ppm, ref_smooth)
    ref_correction = -ref_tsp_ppm if ref_tsp_ppm else 0.0784
    ref_ppm_corrected = ref_ppm + ref_correction
    
    # Integrate TSP in reference
    ref_tsp_area = integrate_peak(ref_ppm_corrected, ref_smooth, (-0.2, 0.2))
    
    # Integrate Alanine CH3 in reference (1.4-1.55 ppm)
    ref_ala_area = integrate_peak(ref_ppm_corrected, ref_smooth, (1.4, 1.55))
    
    # TSP-normalized Alanine area for reference
    ref_ala_normalized = ref_ala_area / ref_tsp_area
    
    print(f"  Reference TSP area: {ref_tsp_area:.2e}")
    print(f"  Reference Alanine CH3 area: {ref_ala_area:.2e}")
    print(f"  Reference Alanine/TSP ratio: {ref_ala_normalized:.4f}")
    print()
    
    # Process each sample
    results = []
    
    print("Processing sample files...")
    print(f"{'File':<10} {'Conc (mM)':<12} {'TSP Area':<12} {'Ala Area':<12} {'Ala/TSP':<12} {'Expected':<12} {'Recovery %'}")
    print("-" * 90)
    
    for item in ala_data:
        file_num = item['file_num']
        concentration = item['concentration_mM']
        
        # Read sample file
        sample_file = base_dir / f"{file_num}.dx"
        if not sample_file.exists():
            print(f"  File {file_num}.dx not found, skipping...")
            continue
        
        sample_ppm, sample_data, sample_sfo1 = read_nmr_data(sample_file)
        
        # Apply smoothing
        sample_smooth = gaussian_filter1d(sample_data, sigma=2)
        
        # Find TSP and apply correction
        sample_tsp_ppm = find_tsp_peak(sample_ppm, sample_smooth)
        sample_correction = -sample_tsp_ppm if sample_tsp_ppm else 0.0784
        sample_ppm_corrected = sample_ppm + sample_correction
        
        # Integrate TSP in sample
        sample_tsp_area = integrate_peak(sample_ppm_corrected, sample_smooth, (-0.2, 0.2))
        
        # Integrate Alanine CH3 in sample
        sample_ala_area = integrate_peak(sample_ppm_corrected, sample_smooth, (1.4, 1.55))
        
        # TSP-normalized Alanine area for sample
        sample_ala_normalized = sample_ala_area / sample_tsp_area
        
        # Calculate scaling factor (ratio of normalized areas)
        # This represents the relative concentration vs reference
        scale_factor = sample_ala_normalized / ref_ala_normalized
        
        # Calculate calculated concentration
        # [M]_sample = [M]_ref * scale_factor * (TSP_ref/TSP_sample) 
        # But since we already normalized by TSP, it's just:
        calculated_conc = 40.0 * scale_factor  # Reference is 40 mM
        
        # Recovery percentage
        recovery = (calculated_conc / concentration) * 100
        
        results.append({
            'file_num': file_num,
            'concentration_mM': concentration,
            'tsp_area': sample_tsp_area,
            'ala_area': sample_ala_area,
            'ala_per_tsp': sample_ala_normalized,
            'scale_factor': scale_factor,
            'calculated_conc': calculated_conc,
            'recovery': recovery
        })
        
        print(f"{file_num:<10} {concentration:<12.3f} {sample_tsp_area:<12.2e} {sample_ala_area:<12.2e} "
              f"{sample_ala_normalized:<12.4f} {calculated_conc:<12.3f} {recovery:.1f}%")
    
    return results, ref_ala_normalized


def plot_linearity(results):
    """Plot linearity of TSP-scaled peak area vs concentration."""
    
    concentrations = [r['concentration_mM'] for r in results]
    ala_per_tsp = [r['ala_per_tsp'] for r in results]
    calculated_concs = [r['calculated_conc'] for r in results]
    
    # Linear regression
    slope, intercept, r_value, p_value, std_err = stats.linregress(concentrations, ala_per_tsp)
    
    # Create figure with 2 subplots
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # Plot 1: TSP-normalized area vs Concentration
    ax1 = axes[0]
    ax1.scatter(concentrations, ala_per_tsp, s=100, c='blue', alpha=0.7, edgecolors='black', zorder=3)
    
    # Fit line
    x_fit = np.linspace(0, max(concentrations) * 1.1, 100)
    y_fit = slope * x_fit + intercept
    ax1.plot(x_fit, y_fit, 'r--', linewidth=2, label=f'Linear fit: y = {slope:.4f}x + {intercept:.4f}')
    
    ax1.set_xlabel('Prepared Concentration (mM)', fontsize=12, fontweight='bold')
    ax1.set_ylabel('Alanine CH₃ Area / TSP Area', fontsize=12, fontweight='bold')
    ax1.set_title('TSP-Scaled Peak Area vs Concentration\n(Alanine Reference Samples)', 
                  fontsize=13, fontweight='bold')
    ax1.legend(loc='upper left', fontsize=10)
    ax1.grid(True, alpha=0.3)
    
    # Add R² annotation
    ax1.text(0.05, 0.95, f'R² = {r_value**2:.4f}\np-value = {p_value:.2e}', 
             transform=ax1.transAxes, fontsize=11, verticalalignment='top',
             bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    # Plot 2: Calculated vs Prepared Concentration
    ax2 = axes[1]
    ax2.scatter(concentrations, calculated_concs, s=100, c='green', alpha=0.7, edgecolors='black', zorder=3)
    
    # 1:1 line
    max_val = max(max(concentrations), max(calculated_concs)) * 1.1
    ax2.plot([0, max_val], [0, max_val], 'r--', linewidth=2, label='1:1 line (ideal)')
    
    ax2.set_xlabel('Prepared Concentration (mM)', fontsize=12, fontweight='bold')
    ax2.set_ylabel('Calculated Concentration (mM)', fontsize=12, fontweight='bold')
    ax2.set_title('Calculated vs Prepared Concentration\n(Alanine Quantification)', 
                  fontsize=13, fontweight='bold')
    ax2.legend(loc='upper left', fontsize=10)
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, max_val)
    ax2.set_ylim(0, max_val)
    
    # Add recovery info
    avg_recovery = np.mean([r['recovery'] for r in results])
    std_recovery = np.std([r['recovery'] for r in results])
    ax2.text(0.05, 0.95, f'Average Recovery: {avg_recovery:.1f} ± {std_recovery:.1f}%', 
             transform=ax2.transAxes, fontsize=11, verticalalignment='top',
             bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.5))
    
    plt.tight_layout()
    plt.savefig('alanine_quantification_linearity.png', dpi=200, bbox_inches='tight')
    print(f"\n✅ Plot saved to: alanine_quantification_linearity.png")
    plt.close()
    
    return r_value**2, slope, intercept


def main():
    print("=" * 80)
    print("ALANINE QUANTIFICATION ANALYSIS")
    print("TSP-Scaled Peak Area vs Concentration Linearity")
    print("=" * 80)
    print()
    
    # Process files
    results, ref_normalized = process_alanine_files()
    
    print()
    print("=" * 80)
    print("LINEARITY ANALYSIS")
    print("=" * 80)
    
    # Plot and calculate linearity
    r_squared, slope, intercept = plot_linearity(results)
    
    print(f"\nLinear Regression Results:")
    print(f"  Slope: {slope:.6f}")
    print(f"  Intercept: {intercept:.6f}")
    print(f"  R²: {r_squared:.4f}")
    
    if r_squared > 0.99:
        print(f"  ✅ Excellent linearity (R² > 0.99)")
    elif r_squared > 0.95:
        print(f"  ✅ Good linearity (R² > 0.95)")
    else:
        print(f"  ⚠️  Moderate linearity (R² < 0.95)")
    
    # Recovery statistics
    recoveries = [r['recovery'] for r in results]
    print(f"\nRecovery Statistics:")
    print(f"  Mean: {np.mean(recoveries):.1f}%")
    print(f"  Std Dev: {np.std(recoveries):.1f}%")
    print(f"  Range: {min(recoveries):.1f}% - {max(recoveries):.1f}%")
    
    print()
    print("=" * 80)
    print("CONCLUSION")
    print("=" * 80)
    print(f"The TSP-scaled peak area shows {'excellent' if r_squared > 0.99 else 'good' if r_squared > 0.95 else 'moderate'} ")
    print(f"linearity with concentration (R² = {r_squared:.4f}).")
    print(f"Average recovery is {np.mean(recoveries):.1f}%, indicating the method is quantitative.")
    print("=" * 80)


if __name__ == "__main__":
    main()