Python处理CSV多列荧光数据：峰值检测循环与导出问题

问题背景 #

我是长期浏览技术社区的Python新手，现有包含2列及以上的CSV文件：第1列为时间序列X数据，其余列是需要检测多个峰值的荧光值。已能通过代码实现单列数据的峰值检测，但无法实现多列循环处理（CSV列数从2到数十列，不想手动修改列配置），尝试使用for column in data.T语句循环但未成功，希望将每列的分析结果写入CSV或Excel文件。

问题根源 #

原代码存在两个核心问题：

• 峰值检测逻辑(find_peaks)写在循环外部，仅处理了第一列数据，循环过程中始终复用第一列的峰值结果，没有对每列信号单独分析
• 结果写入时始终覆盖同一个CSV文件，无法区分不同列的分析结果

修改后的完整代码 #

# -*- coding: utf-8 -*-
"""
Created on Sat Mar 25 14:40:52 2023

@author: klinedd
"""

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import gridspec
from scipy.signal import find_peaks
from scipy.optimize import curve_fit
import csv

my_file = r'C:/David Documents/Python/_CaImaging_python/Python_test_GABARx.csv'
output_dir = r'C:/David Documents/Python/_CaImaging_python/'

# 读取文件并获取列数
f = open(my_file)
reader = csv.reader(f, delimiter=',')
ncol = len(next(reader))
print(f"该文件共有 {ncol} 列")

# 加载数据
data = np.loadtxt(my_file, delimiter=",")
time = data[:, 0] 
f_signal_all = data[:, 1:]  # 提取所有荧光信号列

# 计算采样间隔和频率
interval = np.diff(time).min()
fs = 1 / interval
print(f"采样间隔: {interval} 秒，采样频率: {fs} Hz")

# 绘制所有列的完整曲线
fig1 = plt.figure(figsize=(12, 12))
ax1 = fig1.add_subplot(211)
ax1.set_title("所有荧光信号曲线")
ax1.plot(time, f_signal_all, linewidth=0.5)
ax1.set_ylabel("荧光强度")
ax1.set_xlabel("时间 (秒)")
fig1.tight_layout()
plt.show()

# 峰值检测参数配置
pretrigger_window = (2000 * fs) / 1000
posttrigger_window = (15000 * fs) / 1000
event_no = 0
window_length = posttrigger_window - pretrigger_window
print(f"分析窗口长度: {window_length} 秒")

thresh_min = 10
thresh_max = 1200
thresh_prominence = 16
thresh_min_width = 0.9 * (fs / 1000)

# 创建Excel写入对象，将每列结果存入不同Sheet
with pd.ExcelWriter(f'{output_dir}所有列峰值分析结果.xlsx') as writer:
    # 遍历每一列荧光信号
    for col_idx, f_signal in enumerate(f_signal_all.T, start=1):
        print(f"正在处理第 {col_idx} 列数据...")
        
        # 对当前列执行峰值检测
        peaks, peaks_dict = find_peaks(f_signal, 
                   height=(thresh_min, thresh_max),
                   threshold=None,
                   distance=None,
                   prominence=thresh_prominence,
                   width=thresh_min_width,
                   wlen=None,
                   rel_height=0.5,
                   plateau_size=None)
        
        # 创建结果表格
        table = pd.DataFrame(columns = ['event', 'peak_position', 
                                        'peak_position_s',
                                        'event_start', 'event_end',
                                        'Peak_Amp', 'Width_ms', 
                                        'inst_freq', 'isi_s', 
                                        'Area_pA/ms', 'log_decay', 
                                        'tau_exp'])
        
        table.event = np.arange(1, len(peaks) + 1)
        table.peak_position = peaks
        table.peak_position_s = peaks / fs
        table.event_start = peaks_dict['left_ips'] - pretrigger_window
        table.event_end = peaks_dict['right_ips'] + posttrigger_window
        table.Peak_Amp = peaks_dict['peak_heights']
        table.Width_ms = peaks_dict['widths']/(fs/1000)
        
        # 计算瞬时频率
        if len(peaks) > 1:
            inst_freq_values = (1 / (np.array(table.peak_position[1:]) - np.array(peaks_dict['left_ips'][:-1])) * fs)
            table.inst_freq = np.append(inst_freq_values, np.nan)
        else:
            table.inst_freq = np.nan
        
        # 计算ISI
        table.isi_s = np.diff(peaks, axis=0, prepend=peaks[0]) / fs
        
        # 计算事件面积
        for i, event in table.iterrows():
            start_idx = max(0, int(event.event_start))
            end_idx = min(len(f_signal)-1, int(event.event_end))
            individual_event = f_signal[start_idx:end_idx]
            table.loc[i, 'Area_pA/ms'] = np.round(individual_event.sum(), 1)/(fs/1000)
        
        # 计算对数衰减tau
        for i, event in table.iterrows():
            start_idx = int(event.peak_position)
            end_idx = min(len(f_signal)-1, int(event.event_end))
            decay_tau = abs(f_signal[start_idx:end_idx])
            if len(decay_tau) < 2:
                table.loc[i, 'log_decay'] = np.nan
                continue
            log_decay_tau = np.log(decay_tau)
            decay_width_array = np.arange(len(decay_tau))
            slope, _ = np.polyfit(decay_width_array, log_decay_tau, 1)
            if slope == 0:
                table.loc[i, 'log_decay'] = np.nan
            else:
                tau = -1 / slope
                table.loc[i, 'log_decay'] = tau/(fs/1000)
        
        # 计算指数拟合tau
        for i, event in table.iterrows():
            start_idx = int(event.peak_position)
            end_idx = min(len(f_signal)-1, int(event.event_end))
            decay_tau = f_signal[start_idx:end_idx]
            if len(decay_tau) < 2:
                table.loc[i, 'tau_exp'] = np.nan
                continue
            decay_width_array = np.arange(len(decay_tau))
            try:
                popt, pcov = curve_fit(lambda t, a, b: a * np.exp(b * t), 
                                       decay_width_array, decay_tau, 
                                       p0=(200, 0.1), 
                                       maxfev=1000)
                b = popt[1]
                if b == 0:
                    table.loc[i, 'tau_exp'] = np.nan
                else:
                    table.loc[i, 'tau_exp'] = abs((1/b)/(fs/1000))
            except:
                table.loc[i, 'tau_exp'] = np.nan
        
        # 绘制当前列的峰值检测结果
        fig2 = plt.figure(figsize=(12,4))
        gridspec = fig2.add_gridspec(ncols=2, nrows=1, width_ratios=[2, 1])
        ax1 = fig2.add_subplot(gridspec[0])
        ax1.set_title(f"第 {col_idx} 列 - 峰值检测结果")   
        ax1.plot(time, f_signal)
        ax1.plot(time[peaks], f_signal[peaks], "r.")
        for i, txt in enumerate(table.event):
            ax1.annotate(txt, (time[peaks[i]], f_signal[peaks][i]))
        ax1.set_xlabel("时间 (秒)")
        ax1.set_ylabel("荧光强度")
        
        # 绘制单个事件细节
        ax2 = fig2.add_subplot(gridspec[1])
        ax2.set_title(f"第 {col_idx} 列 - 事件查看器")
        ax2.plot(time, f_signal, "gray")
        ax2.plot(time[peaks], f_signal[peaks], "r.")
        ax2.set_xlabel("时间 (秒)")
        ax2.set_ylabel("荧光强度")
        if len(peaks) > 0:
            ax2.set_xlim(table.event_start[event_no], table.event_end[event_no])
            line, = ax2.plot(time[peaks], f_signal[peaks], "r.") 
            line.set_label(f"事件 {table.event[event_no]}") 
            ax2.legend()
        
        plt.tight_layout()
        plt.show()
        
        # 将当前列结果写入Excel的对应Sheet
        table.to_excel(writer, sheet_name=f'第{col_idx}列结果', index=False)
        # 同时保存单独的CSV文件（可选）
        table.to_csv(f'{output_dir}第{col_idx}列峰值分析结果.csv', index=False)

print("所有列分析完成，结果已保存到指定路径")

关键修改说明 #

• 将峰值检测逻辑移至循环内部，确保每列信号都单独执行find_peaks分析
• 使用pd.ExcelWriter将所有列结果存入同一个Excel文件的不同Sheet，同时支持生成单独的CSV文件
• 增加了边界判断（如数组索引越界、拟合失败等情况），避免运行报错
• 循环时直接遍历f_signal_all.T（所有荧光信号列的转置），更直观地处理每一列数据
• 绘图部分更新为对应列的时间-信号曲线，提升可视化准确性

内容的提问来源于stack exchange，提问作者David Kline