equation7_integration_explained.py

#!/usr/bin/env python3
"""
Equation 7 Implementation using Integration Method
=================================================

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

For INTEGRATION method (Figure S9):
- scale(M) = ∫(M in mixture) / ∫(M in reference)
- scale(TSP) = ∫(TSP in mixture) / ∫(TSP in reference)

Since [TSP]₁ = [TSP]ref (same TSP concentration in all samples):
[M]₁ = [M]ref × scale(M) / scale(TSP)
     = [M]ref × [∫(M in mixture)/∫(M in ref)] / [∫(TSP in mixture)/∫(TSP in ref)]
     = [M]ref × [∫(M in mixture)/∫(TSP in mixture)] / [∫(M in ref)/∫(TSP in ref)]
     = [M]ref × (Ala/TSP ratio in mixture) / (Ala/TSP ratio in reference)
"""

import nmrglue as ng
import numpy as np
import pandas as pd
from scipy.ndimage import gaussian_filter1d
from pathlib import Path


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 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("=" * 90)
    print("EQUATION 7 IMPLEMENTATION USING INTEGRATION METHOD")
    print("=" * 90)
    print()
    print("Paper's Equation 7:")
    print("  [M]₁ = ([M]ref × [TSP]₁ × scale(M)) / ([TSP]ref × scale(TSP))")
    print()
    print("For INTEGRATION method:")
    print("  scale(M) = ∫(M in sample) / ∫(M in reference)")
    print("  scale(TSP) = ∫(TSP in sample) / ∫(TSP in reference)")
    print()
    print("Since [TSP] is constant in all samples:")
    print("  [M]₁ = [M]ref × scale(M) / scale(TSP)")
    print("       = [M]ref × [∫(M in sample)/∫(TSP in sample)] / [∫(M in ref)/∫(TSP in ref)]")
    print("       = [M]ref × (Ala/TSP)_sample / (Ala/TSP)_reference")
    print()
    print("=" * 90)
    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("-" * 90)
    
    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
    
    # Integration regions
    ala_region = (1.4, 1.55)  # Alanine CH3
    tsp_region = (-0.2, 0.2)  # TSP
    
    # Integrate reference
    ref_ala_integral = integrate_peak(ref_ppm_corr, ref_smooth, ala_region)
    ref_tsp_integral = integrate_peak(ref_ppm_corr, ref_smooth, tsp_region)
    ref_ala_per_tsp = ref_ala_integral / ref_tsp_integral
    
    print(f"  [M]ref (known concentration): 40.00 mM")
    print(f"  ∫(Alanine CH3 in ref): {ref_ala_integral:.4e}")
    print(f"  ∫(TSP in ref): {ref_tsp_integral:.4e}")
    print(f"  (Ala/TSP)_ref = {ref_ala_per_tsp:.4f}")
    print()
    
    # ========== SAMPLE SAMPLES ==========
    print("STEP 2: SAMPLE ANALYSIS")
    print("-" * 90)
    print(f"{'File':<8} {'[M]known':<12} {'∫(Ala)':<14} {'∫(TSP)':<14} {'(Ala/TSP)':<12} {'scale(M)':<12} {'scale(TSP)':<12}")
    print("-" * 90)
    
    results = []
    
    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 sample
        sample_ala_integral = integrate_peak(sample_ppm_corr, sample_smooth, ala_region)
        sample_tsp_integral = integrate_peak(sample_ppm_corr, sample_smooth, tsp_region)
        sample_ala_per_tsp = sample_ala_integral / sample_tsp_integral
        
        # Calculate scale factors (Equation 7 terms)
        scale_M = sample_ala_integral / ref_ala_integral  # scale(M)
        scale_TSP = sample_tsp_integral / ref_tsp_integral  # scale(TSP)
        
        print(f"{file_num:<8} {conc_known:<12.3f} {sample_ala_integral:<14.4e} {sample_tsp_integral:<14.4e} "
              f"{sample_ala_per_tsp:<12.4f} {scale_M:<12.4f} {scale_TSP:<12.4f}")
        
        results.append({
            'file_num': file_num,
            'conc_known': conc_known,
            'ala_integral': sample_ala_integral,
            'tsp_integral': sample_tsp_integral,
            'ala_per_tsp': sample_ala_per_tsp,
            'scale_M': scale_M,
            'scale_TSP': scale_TSP
        })
    
    print()
    
    # ========== EQUATION 7 CALCULATION ==========
    print("STEP 3: EQUATION 7 CALCULATION")
    print("-" * 90)
    print(f"{'File':<8} {'[M]known':<12} {'scale(M)':<12} {'scale(TSP)':<12} {'[M]calc':<12} {'Recovery':<10}")
    print("-" * 90)
    
    for r in results:
        # Equation 7: [M]₁ = [M]ref × scale(M) / scale(TSP)
        # Since [TSP]₁ = [TSP]ref, they cancel out
        conc_calc = 40.0 * r['scale_M'] / r['scale_TSP']
        recovery = (conc_calc / r['conc_known']) * 100
        
        print(f"{r['file_num']:<8} {r['conc_known']:<12.3f} {r['scale_M']:<12.4f} {r['scale_TSP']:<12.4f} "
              f"{conc_calc:<12.2f} {recovery:<10.1f}%")
    
    print()
    print("=" * 90)
    print("SUMMARY")
    print("=" * 90)
    print()
    print("Equation 7 terms obtained from INTEGRATION:")
    print()
    print("  [M]ref = 40.00 mM (known from sample preparation)")
    print("  [TSP]₁ = [TSP]ref (same TSP concentration, cancels out)")
    print()
    print("  scale(M) = ∫(Ala in sample) / ∫(Ala in reference)")
    print("           = Peak area ratio (metabolite)")
    print()
    print("  scale(TSP) = ∫(TSP in sample) / ∫(TSP in reference)")
    print("             = Peak area ratio (internal standard)")
    print()
    print("Final calculation:")
    print("  [M]₁ = 40.00 × scale(M) / scale(TSP)")
    print("       = 40.00 × (Ala/TSP)_sample / (Ala/TSP)_reference")
    print()
    print("=" * 90)


if __name__ == "__main__":
    main()