PCWA

基于CWT变换的高度并行且快速的事件检测器。 PCWA 是一种多尺度方法，可以找到任何形状的事件，其母小波与事件形状相匹配（详情请参阅目前正在审稿的 Nature Communications 手稿）。与以前的基于CWT的峰值检测器不同，PCWA 可以与任何用户定义的母小波函数相匹配，Macro- 和 u- 聚类步骤，将大数据分成更小的 M-clusters。聚类步骤利用了 x轴、尺度轴 和 系数值 来提高定位事件的精度。

要求

• Python >= 3.8.5
• numpy >= 1.19.2
• scipy >= 1.6.2
• matplotlib >= 3.3.4
• h5py >= 2.10.0
• pandas >= 1.2.1

可能适用于较旧版本（Python > 3），但撰写此文档时未进行测试。

如何使用PCWA

PCWA 被设计为Python类，需要初始化。导入 pcwa 并初始化一个新实例

import pcwa as pcwa

pcwa_analyzer = pcwa.PCWA(parallel=False)
# pcwa_analyzer.show_wavelets = True
pcwa_analyzer.w, pcwa_analyzer.h = 1.5, 1.5
pcwa_analyzer.selectivity = 0.7
pcwa_analyzer.use_scratchfile = False

可以在初始化期间或之后设置属性。以下是一些属性列表

属性

dt = 1e5                               # sampling period of the signal in s
parallel = True                        # enable/disable multiprocessing 
mcluster = True                        # enable/disable macro-clustering
logscale = True                        # enable/disable logarithmic scale for scale-axis
wavelet = ['ricker']                   # list of wavelet function names
scales = [0.01e-3,0.1e-3,30]           # scale range and count in in s
selectivity = 0.5                      # minimum number of candidates in a valid micro-cluster
w = 2                                  # spreading factor in x-axis
h = 6                                  # spreading factor in y-axis (scale-axis)
extent = 1                             # global extent along x and y axis, used in macro-clustering
trace = None                           # trace (data) variable. 1D numpy vector
events = []                            # list of detected events (valid after calling detect_events() function)
cwt = {}                               # dictionary of cwt coefficients
wavelets = {}                          # dictionary of generated scaled&normalized 1D wavelet arrays
show_wavelets = False                  # plot wavelet functions
update_cwt = True                      # if False, will use the current cwt coefficients to detect events to save time tuning threshold parameters
keep_cwt = False                       # if False, will use less memory by running conv() and local_maxima() at the same time. Otherwise will generate entire CWT coefficient before looking for local maxima (conventional method)
use_scratchfile = False                 # stores cwt coefficients in the scarach file (hdf5 formatted) file

事件检测

初始化后，可以通过调用 detect_events() 方法来检测事件。

events = pcwa_analyzer.detect_events(trace=data,wavelet=['ricker'],scales=[0.1e-3,1.0e-3,50],threshold=3)
tpr, fdr = pcwa.tprfdr(truth,events['loc'],e=7e-3/1e-5,MS=True) # e is the tolerance of error for event location, here 7ms/0.01ms (in data points), 0.01ms is the bin size

在调用 detect_events() 函数时，可以通过传递以下参数来覆盖一些 PCWA 参数。

• trace：覆盖跟踪参数。
• wavelet：覆盖小波函数。
• scales：覆盖尺度。

threshold 是每次调用所必需的唯一参数。

示例代码

以下示例展示了如何使用 PCWA 检测模拟质谱数据中的峰。

import numpy as np
import pandas as pd
import pcwa as pcwa
import matplotlib.pyplot as plt

# read the raw mass scpectroscopy data and truth values
df_raw = pd.read_csv('n100sig66_dataset_1_25/Dataset_14/RawSpectra/noisy22.txt',sep=' ')
df_true = pd.read_csv('n100sig66_dataset_1_25/Dataset_14/truePeaks/truth22.txt',sep=' ')

# create pcwa_analyzer object and set the desired parameters
pcwa_analyzer = pcwa.PCWA()
pcwa_analyzer.trace = df_raw['Intensity']
pcwa_analyzer.dt = 1
pcwa_analyzer.scales = [10,100,100]
pcwa_analyzer.wavelet = ['ricker']
pcwa_analyzer.keep_cwt = False
pcwa_analyzer.w, pcwa_analyzer.h = 0.2, 1
pcwa_analyzer.show_wavelets = False
pcwa_analyzer.use_scratchfile = False

# detect events (peaks)
events = pcwa_analyzer.detect_events(threshold=200)

# fine tune the location of detected peaks
loc = [int(e-events['scale'][n]+np.argmax(df_raw['Intensity'][int(e-events['scale'][n]):int(e+events['scale'][n])])) for n,e in enumerate(events['loc'])]

fig, ax = plt.subplots(3,1,figsize=(16,4),dpi=96,sharex=True,gridspec_kw={'height_ratios': [12,1,1]})
plt.subplots_adjust(hspace=0,wspace=0)
l0, = ax[1].plot(df_true['Mass'],df_true['Particles']*0, '|',markersize=10,color='gray',label='Truth')
ax[0].plot(df_raw['Mass'],df_raw['Intensity'],color='blue')
l1, = ax[2].plot(df_raw['Mass'].iloc[loc],[0]*len(loc),'|',markersize=10,color='red',label='PCWA')
ax[1].set_yticks([])
ax[1].set_ylim(0,0)
ax[2].set_yticks([])
ax[2].set_ylim(0,0)
ax[0].set_ylabel('Intensity')
ax[-1].set_xlabel('m/z')
plt.legend(handles=[l0,l1], bbox_to_anchor=(1.0, 4), loc='upper left')
plt.show()

分析完成后，绘图窗口应显示以下结果

如果您想计算 TPR 和 FDR 值，这对于 ROC 图非常有用，PCWA 类提供了一个函数来完成这项操作。

true_peaks = np.sort(df_true['Mass'].to_numpy())
detected_peaks = np.sort(df_raw['Mass'].iloc[loc].to_numpy())
tpr, fdr = pcwa.tprfdr(true_peaks, detected_peaks, e=0.01, MS=True)
print(f"TPR={tpr:.3f}, FDR={fdr:.3f}")

TPR=0.864, FDR=0.014

命令窗口应显示基于真实值和可接受误差范围（此处为1%）的 TPR 和 FDR 值。 MS 参数确定应用可接受误差的方式，对于质谱数据，错误被认为是相对于质量值的（见下图）。如果 MS=False，则考虑绝对误差值。完整的示例文件提供在本存储库中（ms_example.py）。

参考文献

提供的数据集是从模拟质谱数据集（DOI：10.1093/bioinformatics/bti254）中选取的一个子集。