lineshape_fitting_quantification_principle.py

#!/usr/bin/env python3
"""
Lineshape Fitting Quantification Principle
==========================================

According to NMR colleague:
- Take one sample as reference (e.g., 40 mM)
- Perform lineshape fitting: adjust reference to match sample signal
- The fitted scale factor should grow LINEARLY with concentration
- Scale factor ∝ [M]sample / [M]reference
- This factor can be corrected with TSP factor

Expected relationship:
  fitted_scale ≈ [M]sample / [M]reference
  
  For 20 mM sample vs 40 mM reference: scale ≈ 0.5
  For 10 mM sample vs 40 mM reference: scale ≈ 0.25
  etc.
"""

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


def read_nmr_data(file_path):
    """Read JCAMP-DX file."""
    dic, data_list = ng.jcampdx.read(str(file_path))
    
    data_real = data_list[0]
    data_imag = data_list[1] if len(data_list) > 1 else np.zeros_like(data_real)
    data_magnitude = np.sqrt(data_real**2 + data_imag**2)
    
    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


def find_tsp_peak(ppm, data):
    """Find TSP peak position."""
    mask = (ppm >= -0.5) & (ppm <= 0.5)
    if not np.any(mask):
        return None
    ppm_region = ppm[mask]
    data_region = data[mask]
    peak_idx = np.argmax(data_region)
    return ppm_region[peak_idx]


def integrate_peak(ppm, intensity, region):
    """Integrate peak area."""
    mask = (ppm >= region[0]) & (ppm <= region[1])
    if not np.any(mask):
        return 0
    ppm_region = ppm[mask]
    intensity_region = intensity[mask]
    return abs(np.trapz(intensity_region, ppm_region))


def fit_lineshape(sample_ppm, sample_data, ref_ppm, ref_data, fit_region):
    """
    Fit reference to sample by finding optimal scale factor.
    Returns the scale factor that minimizes residuals.
    """
    # Extract fitting region
    s_mask = (sample_ppm >= fit_region[0]) & (sample_ppm <= fit_region[1])
    r_mask = (ref_ppm >= fit_region[0]) & (ref_ppm <= fit_region[1])
    
    s_ppm = sample_ppm[s_mask]
    s_data = sample_data[s_mask]
    r_ppm = ref_ppm[r_mask]
    r_data = ref_data[r_mask]
    
    # Interpolate reference to match sample ppm grid
    r_interp = np.interp(s_ppm, r_ppm, r_data)
    
    # Find scale factor that minimizes sum of squared residuals
    def residuals(scale):
        return np.sum((s_data - scale * r_interp)**2)
    
    result = minimize_scalar(residuals, bounds=(0.001, 10), method='bounded')
    scale_factor = result.x
    
    # Calculate R²
    fitted = scale_factor * r_interp
    ss_res = np.sum((s_data - fitted)**2)
    ss_tot = np.sum((s_data - np.mean(s_data))**2)
    r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0
    
    return scale_factor, r_squared, s_ppm, s_data, fitted, r_interp


def get_alanine_data():
    """Get Alanine concentration data."""
    df = pd.read_excel('2_Analytical Chemistry Data.xlsx', 
                       sheet_name='Reference Library Data Analysis', header=None)
    
    ala_data = []
    for idx in range(28, 34):
        row = df.iloc[idx]
        file_name = row[1]
        concentration = row[3]
        
        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])
                if file_num != 10:
                    ala_data.append({
                        'file_num': file_num,
                        'concentration_mM': float(concentration)
                    })
    return ala_data


def main():
    print("=" * 100)
    print("LINESHape FITTING QUANTIFICATION PRINCIPLE")
    print("=" * 100)
    print()
    print("PRINCIPLE (from NMR colleague):")
    print("-" * 100)
    print("1. Take one sample as REFERENCE (e.g., 40 mM)")
    print("2. Perform LINESHAPE FITTING: fit reference to sample signal")
    print("3. The fitted SCALE FACTOR should grow LINEARLY with concentration")
    print("4. Scale factor ∝ [M]sample / [M]reference")
    print()
    print("EXPECTED relationship:")
    print("  fitted_scale ≈ [M]sample / [M]reference")
    print()
    print("  Example: Reference = 40 mM")
    print("    - Sample 20 mM: fitted_scale ≈ 0.5")
    print("    - Sample 10 mM: fitted_scale ≈ 0.25")
    print("    - Sample 5 mM:  fitted_scale ≈ 0.125")
    print()
    print("5. This factor can be corrected with TSP factor for accuracy")
    print("-" * 100)
    print()
    
    # Get data
    ala_data = get_alanine_data()
    base_dir = Path("raw_data/Reference_Raw_Date_JCAMP-DX/Alanine-Reference")
    
    # REFERENCE: 40 mM (file 10)
    print("STEP 1: LOAD REFERENCE (40 mM, file 10)")
    print("-" * 100)
    ref_file = base_dir / "10.dx"
    ref_ppm, ref_data = read_nmr_data(ref_file)
    ref_smooth = gaussian_filter1d(ref_data, sigma=2)
    
    # TSP 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_corr = ref_ppm + ref_correction
    
    print(f"  Reference: 40 mM Alanine")
    print(f"  TSP correction: +{ref_correction:.4f} ppm")
    print()
    
    # Fitting region
    fit_region = (1.38, 1.58)
    
    # Process samples
    print("STEP 2: PERFORM LINESHape FITTING")
    print("-" * 100)
    print(f"{'File':<8} {'[M]known':<12} {'Expected':<12} {'Fitted':<12} {'Ratio':<12} {'Status':<15}")
    print("-" * 100)
    
    results = []
    
    for item in ala_data:
        file_num = item['file_num']
        conc_known = item['concentration_mM']
        expected_scale = conc_known / 40.0  # Expected: conc/ref_conc
        
        sample_file = base_dir / f"{file_num}.dx"
        if not sample_file.exists():
            continue
        
        # Read sample
        sample_ppm, sample_data = read_nmr_data(sample_file)
        sample_smooth = gaussian_filter1d(sample_data, sigma=2)
        
        # TSP 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_corr = sample_ppm + sample_correction
        
        # LINESHape FITTING: fit reference to sample
        fitted_scale, r2, _, _, _, _ = fit_lineshape(
            sample_ppm_corr, sample_smooth,
            ref_ppm_corr, ref_smooth,
            fit_region
        )
        
        # Ratio of fitted to expected
        ratio = fitted_scale / expected_scale if expected_scale > 0 else 0
        
        # Status
        if 0.8 <= ratio <= 1.2:
            status = "✓ Good"
        elif 0.5 <= ratio <= 1.5:
            status = "~ Fair"
        else:
            status = "✗ Poor"
        
        results.append({
            'file_num': file_num,
            'conc': conc_known,
            'expected': expected_scale,
            'fitted': fitted_scale,
            'ratio': ratio,
            'r2': r2
        })
        
        print(f"{file_num:<8} {conc_known:<12.3f} {expected_scale:<12.3f} {fitted_scale:<12.3f} {ratio:<12.2f} {status:<15}")
    
    print()
    
    # Analysis
    print("STEP 3: ANALYSIS")
    print("-" * 100)
    
    concs = [r['conc'] for r in results]
    fitted_scales = [r['fitted'] for r in results]
    expected_scales = [r['expected'] for r in results]
    
    # Linear regression: fitted_scale vs concentration
    slope_fitted, intercept_fitted, r2_fitted, _, _ = stats.linregress(concs, fitted_scales)
    
    # Linear regression: expected_scale vs concentration (should be perfect)
    slope_expected, intercept_expected, r2_expected, _, _ = stats.linregress(concs, expected_scales)
    
    print(f"\nFitted scale vs Concentration:")
    print(f"  Slope: {slope_fitted:.6f}")
    print(f"  Intercept: {intercept_fitted:.6f}")
    print(f"  R²: {r2_fitted:.4f}")
    
    print(f"\nExpected scale vs Concentration:")
    print(f"  Slope: {slope_expected:.6f} (should be 1/40 = 0.025)")
    print(f"  Intercept: {intercept_expected:.6f} (should be ~0)")
    print(f"  R²: {r2_expected:.4f}")
    
    print(f"\nConclusion:")
    if r2_fitted > 0.95:
        print(f"  ✓ Lineshape fitting scale factor shows LINEAR relationship with concentration!")
        print(f"    R² = {r2_fitted:.4f}")
    else:
        print(f"  ✗ Lineshape fitting scale factor does NOT show expected linearity")
        print(f"    R² = {r2_fitted:.4f}")
        print(f"  This suggests the data files may not match the expected concentrations.")
    
    # Plot
    fig, axes = plt.subplots(1, 2, figsize=(14, 6))
    
    # Plot 1: Scale factor vs Concentration
    ax1 = axes[0]
    ax1.scatter(concs, expected_scales, s=100, c='green', label='Expected (conc/40)', 
                alpha=0.7, edgecolors='black', marker='o', zorder=3)
    ax1.scatter(concs, fitted_scales, s=100, c='blue', label='Fitted (lineshape)', 
                alpha=0.7, edgecolors='black', marker='s', zorder=3)
    
    # Fit lines
    x_fit = np.linspace(0, max(concs), 100)
    ax1.plot(x_fit, slope_expected * x_fit + intercept_expected, 'g--', 
             linewidth=2, label=f'Expected fit: R²={r2_expected:.4f}')
    ax1.plot(x_fit, slope_fitted * x_fit + intercept_fitted, 'b--', 
             linewidth=2, label=f'Fitted fit: R²={r2_fitted:.4f}')
    
    ax1.set_xlabel('Prepared Concentration (mM)', fontsize=12, fontweight='bold')
    ax1.set_ylabel('Scale Factor (fitted_scale)', fontsize=12, fontweight='bold')
    ax1.set_title('Scale Factor vs Concentration\n(Lineshape Fitting)', 
                  fontsize=13, fontweight='bold')
    ax1.legend(loc='upper left')
    ax1.grid(True, alpha=0.3)
    
    # Plot 2: Fitted vs Expected
    ax2 = axes[1]
    ax2.scatter(expected_scales, fitted_scales, s=100, c='purple', 
                alpha=0.7, edgecolors='black', zorder=3)
    
    # 1:1 line
    max_val = max(max(expected_scales), max(fitted_scales)) * 1.1
    ax2.plot([0, max_val], [0, max_val], 'r--', linewidth=2, label='1:1 line (ideal)')
    
    ax2.set_xlabel('Expected Scale Factor (conc/40)', fontsize=12, fontweight='bold')
    ax2.set_ylabel('Fitted Scale Factor (from lineshape)', fontsize=12, fontweight='bold')
    ax2.set_title('Fitted vs Expected Scale Factor', fontsize=13, fontweight='bold')
    ax2.legend(loc='upper left')
    ax2.grid(True, alpha=0.3)
    ax2.set_xlim(0, max_val)
    ax2.set_ylim(0, max_val)
    
    plt.tight_layout()
    plt.savefig('lineshape_fitting_principle.png', dpi=200, bbox_inches='tight')
    print(f"\n✅ Plot saved to: lineshape_fitting_principle.png")
    
    print()
    print("=" * 100)
    print("INTERPRETATION")
    print("=" * 100)
    print()
    print("IDEAL scenario (what NMR colleague expects):")
    print("  - Fitted scale factor should match Expected scale factor")
    print("  - Both should be linear with concentration")
    print("  - R² should be close to 1.0")
    print()
    print("ACTUAL result:")
    if r2_fitted > 0.95:
        print(f"  ✓ Lineshape fitting works! Fitted scale ∝ concentration")
        print(f"    R² = {r2_fitted:.4f}")
    else:
        print(f"  ✗ Lineshape fitting scale does NOT match expected concentration ratio")
        print(f"    R² = {r2_fitted:.4f}")
        print()
        print("  POSSIBLE REASONS:")
        print("  1. The .dx files don't correspond to the Excel concentrations")
        print("  2. These are different samples than the reference library")
        print("  3. Data preprocessing issue in JCAMP-DX export")
    print()
    print("=" * 100)


if __name__ == "__main__":
    main()