MEMORIES/check_pure_shift_pattern.py at main · HolobiomicsLab/MEMORIES · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
#!/usr/bin/env python3
"""
Check if odd-numbered files are pure shift and even-numbered are regular
"""

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 compare_odd_even():
    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', '')))

    odd_files = [f for f in files if int(f.replace('.dx', '')) % 2 == 1]
    even_files = [f for f in files if int(f.replace('.dx', '')) % 2 == 0]

    print(f"Odd files (potential pure shift): {odd_files}")
    print(f"Even files (potential regular): {even_files}")

    # Compare a few pairs
    fig, axes = plt.subplots(3, 2, figsize=(16, 12))

    pairs_to_compare = [
        ('9.dx', '10.dx'),
        ('13.dx', '12.dx'),  # Note: order swapped to match concentration
        ('15.dx', '16.dx'),
    ]

    for idx, (odd_f, even_f) in enumerate(pairs_to_compare):
        # Odd file (pure shift?)
        ax = axes[idx, 0]
        try:
            ppm, spec = read_and_process(os.path.join(folder_path, odd_f))
            tsp = find_tsp_peak(ppm, spec)
            ppm_corr = ppm - tsp
            mask = (ppm_corr >= 0.5) & (ppm_corr <= 5.5)
            ax.plot(ppm_corr[mask], spec[mask], linewidth=0.8, color='blue')
            ax.set_xlim(5.5, 0.5)
            ax.set_title(f'{odd_f} (Odd - Pure Shift?)', fontsize=12, fontweight='bold')
            ax.set_xlabel('ppm')
            ax.set_ylabel('Intensity')
            ax.grid(True, alpha=0.3)
        except Exception as e:
            ax.text(0.5, 0.5, f'Error: {e}', ha='center', va='center', transform=ax.transAxes)

        # Even file (regular?)
        ax = axes[idx, 1]
        try:
            ppm, spec = read_and_processing(os.path.join(folder_path, even_f))
            tsp = find_tsp_peak(ppm, spec)
            ppm_corr = ppm - tsp
            mask = (ppm_corr >= 0.5) & (ppm_corr <= 5.5)
            ax.plot(ppm_corr[mask], spec[mask], linewidth=0.8, color='red')
            ax.set_xlim(5.5, 0.5)
            ax.set_title(f'{even_f} (Even - Regular?)', fontsize=12, fontweight='bold')
            ax.set_xlabel('ppm')
            ax.set_ylabel('Intensity')
            ax.grid(True, alpha=0.3)
        except Exception as e:
            ax.text(0.5, 0.5, f'Error: {e}', ha='center', va='center', transform=ax.transAxes)

    plt.suptitle('Odd (Pure Shift) vs Even (Regular) Comparison', fontsize=14, fontweight='bold')
    plt.tight_layout()
    plt.savefig('odd_vs_even_comparison.png', dpi=150, bbox_inches='tight')
    print("\nSaved comparison to: odd_vs_even_comparison.png")
    plt.close()

    # Check linewidths to confirm
    print("\n" + "="*60)
    print("Analyzing peak linewidths to verify pure shift vs regular...")
    print("="*60)

    for odd_f, even_f in [('9.dx', '10.dx'), ('13.dx', '12.dx')]:
        try:
            # Odd file
            ppm, spec = read_and_process(os.path.join(folder_path, odd_f))
            tsp = find_tsp_peak(ppm, spec)
            ppm_corr = ppm - tsp
            # Look at a specific peak region (e.g., ~1.3 ppm for lactate)
            mask = (ppm_corr >= 1.2) & (ppm_corr <= 1.4)
            if np.sum(mask) > 0:
                peak_idx = np.argmax(spec[mask])
                # Estimate FWHM roughly
                peak_region = spec[mask]
                half_max = np.max(peak_region) / 2
                above_half = peak_region > half_max
                fwhm_points = np.sum(above_half)
                print(f"{odd_f}: Peak at ~1.3 ppm, FWHM ~{fwhm_points} points")

            # Even file
            ppm, spec = read_and_process(os.path.join(folder_path, even_f))
            tsp = find_tsp_peak(ppm, spec)
            ppm_corr = ppm - tsp
            mask = (ppm_corr >= 1.2) & (ppm_corr <= 1.4)
            if np.sum(mask) > 0:
                peak_idx = np.argmax(spec[mask])
                peak_region = spec[mask]
                half_max = np.max(peak_region) / 2
                above_half = peak_region > half_max
                fwhm_points = np.sum(above_half)
                print(f"{even_f}: Peak at ~1.3 ppm, FWHM ~{fwhm_points} points")

        except Exception as e:
            print(f"Error analyzing {odd_f}/{even_f}: {e}")

if __name__ == "__main__":
    compare_odd_even()