equation7_lineshape_fitting_explained.py

#!/usr/bin/env python3
"""
Equation 7 Implementation using Lineshape Fitting (Figure S8 Method)
====================================================================

Paper's Equation 7:
[M]₁ = ([M]ref × [TSP]₁ × scale(M)) / ([TSP]ref × scale(TSP))

For LINESHape FITTING method (Figure S8):
- scale(M) = optimal scaling factor from fitting reference to sample
- But spectra must be TSP-normalized FIRST to account for λSS

Procedure:
1. TSP-normalize both sample and reference (divide by TSP area)
2. Fit TSP-normalized reference to TSP-normalized sample
3. The fitted scale factor = scale(M) / scale(TSP) effectively
   (because TSP normalization already accounts for scale(TSP))
4. [M]₁ = [M]ref × fitted_scale

Alternative interpretation (if not pre-normalizing):
- scale(M) = fitted_scale (from lineshape fitting)
- scale(TSP) = ∫TSP_sample / ∫TSP_ref (from integration)
- [M]₁ = [M]ref × fitted_scale / scale(TSP)
"""

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 spectrum to sample by minimizing residuals.
    Both spectra should already be on the same ppm grid or interpolated.
    """
    # 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)
    
    # Remove baseline offset for better fitting
    s_centered = s_data - np.mean(s_data[:50])  # Use first 50 points as baseline
    r_centered = r_interp - np.mean(r_interp[:50])
    
    # Find scale factor that minimizes sum of squared residuals
    def residuals(scale):
        return np.sum((s_centered - scale * r_centered)**2)
    
    result = minimize_scalar(residuals, bounds=(0.001, 10), method='bounded')
    scale_factor = result.x
    
    # Calculate R²
    fitted = scale_factor * r_centered
    ss_res = np.sum((s_centered - fitted)**2)
    ss_tot = np.sum((s_centered - np.mean(s_centered))**2)
    r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0
    
    return scale_factor, r_squared, s_ppm, s_centered, fitted


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("EQUATION 7 IMPLEMENTATION USING LINESHape FITTING (Figure S8)")
    print("=" * 100)
    print()
    print("Paper's Equation 7:")
    print("  [M]₁ = ([M]ref × [TSP]₁ × scale(M)) / ([TSP]ref × scale(TSP))")
    print()
    print("For LINESHape FITTING method:")
    print()
    print("  OPTION 1: TSP-normalize first, then fit")
    print("    - TSP-normalize: divide spectrum by ∫TSP")
    print("    - Fit normalized reference to normalized sample")
    print("    - fitted_scale = scale(M) / scale(TSP)  [TSP terms already normalized out]")
    print("    - [M]₁ = [M]ref × fitted_scale")
    print()
    print("  OPTION 2: Fit raw spectra, then TSP-correct")
    print("    - Fit raw reference to raw sample → fitted_scale (this is scale(M))")
    print("    - Calculate scale(TSP) = ∫TSP_sample / ∫TSP_ref (from integration)")
    print("    - [M]₁ = [M]ref × fitted_scale / scale(TSP)")
    print()
    print("  [TSP]₁ = [TSP]ref = 0.95 mM (from Sample Preparation sheet)")
    print("  Since [TSP]₁ = [TSP]ref, they cancel out!")
    print()
    print("=" * 100)
    print()
    
    # Get data
    ala_data = get_alanine_data()
    base_dir = Path("raw_data/Reference_Raw_Date_JCAMP-DX/Alanine-Reference")
    
    # ========== REFERENCE SAMPLE (40 mM) ==========
    print("STEP 1: REFERENCE SAMPLE (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
    
    # Integrate TSP for normalization
    ref_tsp_area = integrate_peak(ref_ppm_corr, ref_smooth, (-0.2, 0.2))
    
    # TSP-normalize reference
    ref_normalized = ref_smooth / ref_tsp_area
    
    print(f"  [M]ref (from Excel): 40.00 mM")
    print(f"  [TSP]ref (from Sample Prep): 0.95 mM")
    print(f"  TSP area in ref: {ref_tsp_area:.4e}")
    print(f"  TSP-normalized reference created")
    print()
    
    # Fitting region
    fit_region = (1.38, 1.58)  # Alanine CH3
    
    # ========== OPTION 1: TSP-normalize then fit ==========
    print("STEP 2: OPTION 1 - TSP-Normalize First, Then Fit")
    print("-" * 100)
    print(f"{'File':<8} {'[M]known':<12} {'TSP Area':<14} {'Scale(fit)':<14} {'[M]calc':<12} {'Recovery':<10}")
    print("-" * 100)
    
    results_opt1 = []
    
    for item in ala_data:
        file_num = item['file_num']
        conc_known = item['concentration_mM']
        
        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
        
        # Integrate TSP
        sample_tsp_area = integrate_peak(sample_ppm_corr, sample_smooth, (-0.2, 0.2))
        
        # TSP-normalize sample
        sample_normalized = sample_smooth / sample_tsp_area
        
        # Fit TSP-normalized reference to TSP-normalized sample
        fitted_scale, r2, _, _, _ = fit_lineshape(
            sample_ppm_corr, sample_normalized,
            ref_ppm_corr, ref_normalized,
            fit_region
        )
        
        # Calculate concentration
        # fitted_scale already incorporates scale(M)/scale(TSP) because of TSP normalization
        conc_calc = 40.0 * fitted_scale
        recovery = (conc_calc / conc_known) * 100
        
        results_opt1.append({
            'file_num': file_num,
            'conc_known': conc_known,
            'tsp_area': sample_tsp_area,
            'fitted_scale': fitted_scale,
            'r2': r2,
            'conc_calc': conc_calc,
            'recovery': recovery
        })
        
        print(f"{file_num:<8} {conc_known:<12.3f} {sample_tsp_area:<14.4e} {fitted_scale:<14.4f} {conc_calc:<12.2f} {recovery:<10.1f}%")
    
    print()
    
    # ========== OPTION 2: Fit raw, then TSP-correct ==========
    print("STEP 3: OPTION 2 - Fit Raw Spectra, Then TSP-Correct")
    print("-" * 100)
    print(f"{'File':<8} {'[M]known':<12} {'Scale(fit)':<14} {'Scale(TSP)':<14} {'Ratio':<12} {'[M]calc':<12} {'Recovery':<10}")
    print("-" * 100)
    
    results_opt2 = []
    
    for item in ala_data:
        file_num = item['file_num']
        conc_known = item['concentration_mM']
        
        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
        
        # Fit RAW reference to RAW sample (get scale(M))
        scale_M, r2, _, _, _ = fit_lineshape(
            sample_ppm_corr, sample_smooth,
            ref_ppm_corr, ref_smooth,
            fit_region
        )
        
        # Get scale(TSP) from integration
        sample_tsp_area = integrate_peak(sample_ppm_corr, sample_smooth, (-0.2, 0.2))
        scale_TSP = sample_tsp_area / ref_tsp_area
        
        # Calculate concentration using full Equation 7
        conc_calc = 40.0 * scale_M / scale_TSP
        recovery = (conc_calc / conc_known) * 100
        
        results_opt2.append({
            'file_num': file_num,
            'conc_known': conc_known,
            'scale_M': scale_M,
            'scale_TSP': scale_TSP,
            'conc_calc': conc_calc,
            'recovery': recovery
        })
        
        print(f"{file_num:<8} {conc_known:<12.3f} {scale_M:<14.4f} {scale_TSP:<14.4f} {scale_M/scale_TSP:<12.4f} {conc_calc:<12.2f} {recovery:<10.1f}%")
    
    print()
    
    # ========== COMPARISON ==========
    print("=" * 100)
    print("COMPARISON OF METHODS")
    print("=" * 100)
    
    # Calculate statistics
    concs = [r['conc_known'] for r in results_opt1]
    
    # Option 1
    calc_opt1 = [r['conc_calc'] for r in results_opt1]
    rec_opt1 = [r['recovery'] for r in results_opt1]
    _, _, r2_opt1, _, _ = stats.linregress(concs, calc_opt1)
    
    # Option 2
    calc_opt2 = [r['conc_calc'] for r in results_opt2]
    rec_opt2 = [r['recovery'] for r in results_opt2]
    _, _, r2_opt2, _, _ = stats.linregress(concs, calc_opt2)
    
    print(f"\nOPTION 1 (TSP-normalize then fit):")
    print(f"  R² = {r2_opt1:.4f}")
    print(f"  Mean recovery = {np.mean(rec_opt1):.1f}%")
    
    print(f"\nOPTION 2 (Fit raw, then TSP-correct):")
    print(f"  R² = {r2_opt2:.4f}")
    print(f"  Mean recovery = {np.mean(rec_opt2):.1f}%")
    
    print()
    print("=" * 100)
    print("SUMMARY: LINESHape FITTING TERMS")
    print("=" * 100)
    print()
    print("For Equation 7 using lineshape fitting:")
    print()
    print("  [M]ref = 40.00 mM (from Excel)")
    print("  [TSP]₁ = 0.95 mM (from Sample Preparation)")
    print("  [TSP]ref = 0.95 mM (from Sample Preparation)")
    print()
    print("  OPTION 1:")
    print("    scale(M)/scale(TSP) = fitted_scale (from TSP-normalized fitting)")
    print("    [M]₁ = 40.00 × fitted_scale")
    print()
    print("  OPTION 2:")
    print("    scale(M) = fitted_scale (from raw spectrum fitting)")
    print("    scale(TSP) = ∫TSP_sample / ∫TSP_ref (from integration)")
    print("    [M]₁ = 40.00 × fitted_scale / scale_TSP")
    print()
    print("=" * 100)


if __name__ == "__main__":
    main()