SPKtool
A tool for neurophysiological data processing
锋电位分析工具
版本1.0 2008-8 1 简介及安装
1.1简介
SPKtool是一个在MATLAB平台上编写的神经元锋电位分析工具软件,提供了从锋电位检测、特征提取、聚类分析到放电序列分析的各种功能:
锋电位检测:
-
显示记录信号中锋电位峰值的分布曲线;
-
通过阈值检测锋电位;
-
通过原始信号、能量、非线性能量进行检测;
-
三种锋电位波形对齐方式:阈值、波谷、波峰。
-
GUI交互式界面。
特征提取:
-
多种常用特征:PCA,峰值、能量等;
-
小波方法;
-
特征分别图,特征密度分布图等。
聚类分析:
-
K-means;
-
模板匹配;
-
基于高斯混合模型的EM算法;
-
Valley seeking;
-
通过绘制轮廓线,手动聚类;
-
对聚类结果进行合并,删除等。
动作电位串分析:
-
相邻动作电位间隔分布图(ISI);
-
Poincare map;
-
相关函数图;
-
锋电位分布光栅图;
-
放电频次直方图;
-
刺激前后放电直方图(PSTH);
-
刺激前后放电光栅图。
1.2安装及界面
本程序在MATLAB 7.1 SP3上开发完成,可能需要该版本或更新的版本的MATLAB才能运行。
将程序压缩包拷贝到MATLAB中的work目录下,解压缩,产生SPKtool目录。在MATLAB中,选择菜单File->Set Path … ,在“Set Path”对话框中选择“Add with Subfolders…”,把该目录及子目录添加到MATLAB的路径中去,选择“Save”按钮保存。
在命令窗口输入:spktool,打开GUI界面如下图。
程序界面中包含三个部分,菜单、信息表格和底部的状态栏。软件的所有功能均通过菜单命令来操作。信息表格列出了分析数据的主要信息,如采样频率,记录长度,文件名称、锋电位数目、特征名称、各个分类包含的锋电位数目等。对于多通道数据,通过信息表格右上角的Channel下拉菜单选择当前通道。状态栏给出程序运行状态的简明信息。
2 数据文件
SPKtool没有和其它软件进行数据交换的接口,不能直接读入记录系统产生的数据。用户需要首先将数据(连续记录信号、锋电位波形矩阵等)导入MATLAB,并按照特定的变量名保存到.mat文件(MATLAB中默认的数据文件格式),才能在SPKtool中打开,进行分析。
在菜单中选择File->Load…,打开文件选择对话框。加载数据后,根据文件中的变量,自动设置各个菜单的状态。
SPKtool通过变量名区分数据,因此,对数据文件(.mat)中的变量名有严格的要求。
2.1原始记录数据
原始记录数据文件中至少应包含两个变量:
dataflow:连续记录数据,二维数组,行对应通道,列对应采样点;
fs: 采样频率,数值。
默认情况下,SPKtool对读入的数据进行低通滤波,使用6阶Butterworth滤波器,截至频率4500Hz。可以打开MLoadRaw.m,找到butter函数对应的位置,修改滤波参数。
2.2锋电位波形数据
从原始信号中检测出的锋电位波形数据,至少应该包含以下三个变量:
posind: 锋电位波形在连续数据中的采样位置标记。行向量。第n个波形发生的时刻为posind(n)/fs。
waveforms: 锋电位波形数据,二维数组。每行对应一个锋电位的采样数据。例如,1000个波形,每个波形含40个数据点, 则waveforms尺度为1000×40;
fs: 采样频率。
数据文件中还可以包含如下一些变量:
featurespace: 锋电位特征空间,二维数组,每行包含一个波形的特征。如1000个波形的PC1和PC2特征,则featurespace为1000×2。特征的数目不限。
featnames:特征名称,包含字符串的单元数组,每个字符串对应一个特征的名称。例如: 两个特征的名称分别为PC1和PC2, 则 featnames 为{‘PC1’,’PC2’};
unitcell: 包含分类信息的单元数组。每个单元对应一个类别中的锋电位序号列表。例如,unitcell{1}对应没有归类的锋电位编号。 unitcell{3} 是第二个类别中的锋电位编号数组。
对于多通道数据,上面这些变量需要放在一个名为ResultSet的结构体数组内,每个结构体对应一个通道。
刺激信息
stimmark: 刺激的开始和结束时间。行向量。例如[20 30 60 70],表示两次刺激,第一次刺激从第20s开始,到30s结束,第二次刺激60s开始,70s结束。
上述变量可以共存于一个文件中。在程序中,如果找不到dataflow这个变量,将关闭锋电位检测功能(不能激活Detect菜单),找不到waveforms,将关闭特征提取功能(不能激活Features菜单)。
另外,程序中用到了两个全局变量,采样频率fs和定义各个类别色彩的矩阵colorary,所以在程序运行过程中,请勿执行clear all等能够清除全局变量的命令。
2.3数据保存和导出
选择菜单File->Save …,保存分析结果。文件格式同样是.mat文件,变量名称与读入数据相同。(注意,连续记录数据不包含在该数据文件内。)可以用SPKtool重新打开保存的文件,也可以直接用读入MATLAB,将特定的数据导出。
选择菜单File->Save Current Channel…,仅保存当前通道的分析数据。
菜单File-> Experimental Info …,设置试验信息,包含记录位置、时间、坐标、数据文件名称等信息。
菜单File-> Export to TXT file …,导出数据到文本文件。这将产生两个文本格式的文件,后缀名为.inf的文件包含上面的试验信息,以及数据信息(通道数、单元数等)。后缀名为.txt的文件是一个数据表格,第一列对应锋电位时间(即变量posind),第二列为类别编号,后面的各列为波形数据。
3 锋电位检测
如果读入的数据中包含原始连续记录信号(即dataflow变量),菜单Detect可用,其中包含以下两个命令:
3.1峰值密度分布曲线
菜单Detect->Plot Peak Histogram,绘出原始记录数据的峰值密度分布曲线。幅度为对数坐标。
3.2锋电位检测界面
菜单Detect->Load Detector…,加载锋电位检测,界面如下:
参数设置:
Detect On:该下拉菜单中包含了三种波形检测方式,即分别在原始信号(Raw Signal),原始信号的能量(Energy)和非线性能量(Nonlinear Energy)1曲线上进行检测;
Threshold:阈值。输入阈值后左图中的红色水平线也相应调整。信号中与阈值线有交点的波形将根据后面的参数进行截取;
Spike length:每个锋电位波形的长度,单位为ms。如果采样频率为40KHz,1ms长度的锋电位波形包含40个数据点;
Prethreshold ratio:对齐点以前数据长度比例。对齐点前保留的数据长度与波形数据总长度的比值。例如包含40个点的锋电位波形,对齐方式为valley时,在波谷位置前保留10个采样点,则该参数需要设置为10/40=0.25;
Dead time:不应期长度,单位为ms。例如,如果采样频率为40KHz,该参数设置为2时,程序在检测到一个波形后,将跳过该波形后紧邻的80个点,然后才开始寻找下一个波形。
Alignment:波形对齐方式,阈值位置(gate)、波谷 (valley)、和波峰 (peak) 对齐三种形式。阈值和波谷的对比如下:
a) gate b) valley
使用方法:
中间的坐标轴根据Detect On的设置,显示原始信号或其(非线性)能量曲线。底部坐标轴在(非线性)能量检测方式下显示原始信号波形。左侧坐标轴显示中间坐标轴中数据的峰值分布曲线(横坐标为对数坐标)。
可以利用顶部的工具栏,对中间坐标轴进行缩放和平移,左侧坐标轴的纵坐标和底部坐标轴的横坐标也将进行相应的调整。
参数设置完毕后,选择Detect进行波形检测,检测完成后,将打开一个单独的绘图窗口显示锋电位波形。如果不符合要求,可以修改参数,重新检测。最后,选择Done完成,并关闭该界面。
需要说明的是,由于原始记录信号中的数据量通常很大,为了加快运行速度,坐标轴中只绘出前0.4e6个采样点。当然,锋电位波形检测还是针对所有的采样数据。
右侧的参数中除Detect On外,具有自动保存功能,再次运行Detector时将自动加载前一次设置的参数,即使你前一次选择Cancel退出。
在切换不同的检测方式时,阈值将根据中坐标轴中的数据自动调整为:均值+/-3倍标准差。
4 锋电位特征提取
通过阈值检测出锋电位波形后,或者读入的数据中保护锋电位数据,这时,Feature菜单将可用。在该菜单中,可以进行锋电位特征的选择、计算、修改、绘图等功能。通常需要定义特征后,才能打开Sort菜单进行聚类操作。
4.1特征定义
菜单 Features->Feature Define…,打开如下的特征选择对话框:
左侧All Features 列表中列出了所有可用的特征。中间有两个特征列表:Features Vs Features 中的特征用于绘制特征比较图(菜单 Features-> Features Vs Features);Features For Cluster 才是用于聚类分析的特征。本说明中,除非特别指明,特征均指用于聚类的特征。在右侧可以修改一些特征相关参数。
可用特征:
peak:锋电位波峰处的幅值;
peak.pos:峰值对应的位置;
valley:锋电位波谷处的幅值;
valley.pos:波谷对应的位置;
peak-valley:峰峰值,波峰到波谷的振幅。
width:波形宽度,即波峰位置到波谷位置的距离:peak.pos-valley.pos;
energy:波形的能量,即采样点的平方和。
nonlinear.energy:所有点的非线性能量之和,幅度平方与频率平方的乘积,其计算方法为:NEO(t) = v(t)2- v(t-1)*v(t+1),(v(t)为记录信号);
Slice:锋电位波形中某个时刻的幅值,对应参数Slice Parameters。例如每个锋电位含32个采样点,Slice Parameters为0.5,则Slice特征为锋电位中第32×0.5=16个采样点的数值。Slice Parameters可以包含多个值,提取多个采样点,例如 0.5,0.75 则第16和第24个采样点的数值都作为特征。
pca:主成份分析方法,可以在右侧PCA Parameters中设置要提取的主成份序号;
ridge.pca:基于小波脊的PCA特征提取方法,首先对波形逐个进行小波变换(使用mexh小波),然后提取出小波脊,再对小波脊上的小波系数进行PCA,提取特征。需要三个参数:小波分析的频率范围Frequency limits,小波脊的数目Number of Ridges,以及需要提取的主成份序号PCA Parameters。
完成特征选取后,选择Done,开始提取用于聚类分析的特征,并给出特征分布散点图。
4.2特征修改
特征比例缩放
Features->Scale Features,打开如下图的对话框,在其中,可以调整特征的相对范围,可能有利于改善K-Means的聚类结果。
删除特征
Features->Remove Features,删除部分使用中的特征。一般特征的计算不需要花很长时间,重新选择并计算即可。在同时使用Ridge.pca和PCA作为特征时,由于共用PCA参数,可以通过该工具重新选择。
4.3特征分布散点图
Features->Plot Features 绘制特征分布散点图,如果定义了两个特征,散点图是二维的,如果有三个或更多的特征,则只使用前三个,以三维形式显示,并同时画出特征向三个平面上的投影分布,如下图。
4.4特征密度图
Features-> Features Density:特征分布密度图。如下图,该图只能绘出前两个特征的分布密度。
4.5特征对比图
Features-> Features Vs Features,绘制特征分布的对比图,注意该功能中使用的特征与聚类分析的特征是独立的。并且,每次绘图都需要重新计算特征。
如果已经划分类别,不同类别的特征将以不同颜色区别显示。
4.6主成份
Features-> Principal Components:给出特征提取中主成份曲线。例如特征定义中使用第一个和第二个主成份,则给出PC1和PC2的形状曲线。
5 聚类分析
聚类功能集中于菜单Sort中。在选择具体的聚类方法之前,应该先设置聚类的参数。
5.1参数设置
菜单Sort->Options…, 打开参数设置对话框:
其中共有四个参数:
Unit Number:划分类别的数目,K-Means,期望最大化(EM)方法,以及模板匹配方法使用该参数;
K-Means Radius:K-Means方法的聚类半径,值越小,聚类的半径越小。模板匹配方法也受该参数影响;
Parzen Mult:寻谷聚类算法的初始半径大小,初始半径越小,获得的类别数目越多。
Sort On:设置聚类的锋电位范围,选择All Spikes时对所有锋电位的特征进行聚类;选择Unsorted Spikes时仅对未分类的锋电位进行聚类。
5.2手动聚类
菜单Sort->Contour…,对应手动聚类,选择该菜单后,弹出如下图的手动聚类对话框。
对话框中间给出了未分类的特征的分布图。底部有四个控制按钮。首先单击New Contour按钮,光标将变成图中所示的十字形,在特征分布图中单击左键,将出现一个红色的十字形标记点,这样依次点击,就可以确定一个类别的大致轮廓,之后,鼠标右击,这些点将自动连接成一条封闭的曲线。然后选择底部的Update按钮,曲线内的点就被划分到一个类别。再选择New Contour绘制下一个类别的边界,依此操作,直到聚类完成。
也可以连续绘制多条边界线,然后再选择Update,显示聚类结果。最后选择Done将显示的聚类结果传回主程序,选择Cancel放弃。
手工聚类的类别数目不受前面聚类参数对话框中的类别数目的影响。
5.3自动聚类方法
SPKtool中包含了几种常见的聚类算法,简单说明如下:
K-Means
根据特征之间的距离来进行聚类。用户指定聚类的数目(聚类参数对话框中),该算法首先从样本(特征)中随机抽取确定聚类中心,然后计算每个样本到各个聚类中心的距离,并将其归到距离最近的类中,第一次聚类完成后,计算各个样本的均值作为新的聚类中心,重新计算每个样本到新的聚类中心的距离,并将其划分到聚类最近的类别中,如此循环,直到聚类中心的位置稳定[2]。
Template 模板匹配法
在SPKtool中,模板方法的与K-Means基本上是一样的,不同之处仅在于直接用锋电位的波形数据作为特征,而K-Means使用锋电位降维后的特征数据。
Gaussian EM 基于高斯混合模型的期望最大化算法
该算法首先用K-Means方法确定初始聚类,然后将锋电位特征的分布看成符合混合高斯分布,用期望最大化算法确定高斯混合模型的均值、协方差等参数,然后计算每个特征在各个高斯分布中的概率,并划分到概率最大的一个类别(高斯分布)中[2]。
Valley Seeking 寻谷聚类算法
在特征空间中寻找特征分布密度的极小位置,作为聚类的边界。算法的流程可以简述如下3:
1. 确定初始半径,可以取为样本标准差的某个倍数;初始半径的设置影响聚类的数目;
2. 在该半径范围内,确定每个特征的近邻数目,以及近邻的集合。如果特征的数目为N,需要初始化一个N×N的数组来存放每个特征的近邻编号,计算完成后,可以根据最大近邻集合中元素的数目,删去多余的元素。
3. 如果一个特征的近邻数目,超过它每一个近邻各自的近邻数目,则该特征确定了一个特征密度峰。如此,设置一个最小近邻标准,可以找出全部峰。
4. 五个峰确定5个类别,给每个峰的近邻特征设置类别标记;
5. 将所有特征按照近邻数目从大到小排序;
6. 从近邻数目最大的特征开始,看其近邻中属于哪一类别的最多,就将该特征也归入这一类;如此,完成全部特征的聚类。
5.4聚类结果调整
移除未分类的锋电位
菜单Sort->Edit units …中,选择复选框Remove unsorted spikes。未分类的锋电位将从数据中删去。
合并分类
同样在菜单Sort->Edit units …中,在Merge Units框中输入欲合并的单元编号;例如,输入“2 3”将单元2和单元3合并为单元2,输入“0 3”删除单元3(其中的波形转为未归类状态)。
移除特征离散度很大的锋电位
可以首先通过轮廓线聚类方法,将要保留的特征划分到单元1中,然后按照上面的方法,先移除没有归类的边缘特征,再将单元1与0(未分类)合并。
6 绘图功能
6.1锋电位波形
菜单:Plot->Spikes:绘出所有锋电位波形;
菜单:Plot-> Sorted spikes:分布绘出各个类别的波形。每个波形图左下的数字标明了该分类中波形的数目。
6.2特征分布图
菜单:Plot-> Features:绘出特征分布散点图;
菜单:Plot->Sorted Features:画出分类后的特征散点图;如下图。
6.3时间间隔(ISI)分布
菜单:Plot->ISI …,打开如下的ISI参数设置窗口:
Min. Interval: 最小间隔时间;
Max. Interval: 最大间隔时间;
Number of Bins: 阶梯图的分段数目,数值越大,ISI曲线越平滑;
Units: 选择单元,只画出指定单元的ISI分布图。
6.4 间隔-时间分布图
菜单:Plot->ISI versus time…,绘制相邻放电间隔的时间分布图。参数设置对话框:
X Min. :横坐标(时间)最小值;
X Max.:横坐标最大值;
Units:选择需要进行绘图的单元编号。
6.5庞加莱图
菜单:Plot-> Poincare maps:该图以与前一个放电的时间间隔为纵坐标,与后一个放电的时间间隔为横坐标,标出每一个放电的位置。坐标轴的范围(即最大时间间隔)由ISI中的Min. Interval和Max Interval 指定。如下图。
6.6放电频次阶梯图
菜单:Plot->Fire Rate…:绘制放电的频次直方图。参数设置对话框:
X Min. :横坐标(时间)最小值;
X Max.:横坐标最大值;
Time per Bin:每个阶梯表示的时间长度;
Units:选择需要进行绘图的单元编号。
Add stimulation markers:选中该项后,可以在下面的Markers中输入刺激标记(对应数据文件中的simmark变量)。例如,刺激开始的时刻为49.89885s,结束的时刻为262.4599s,在Markers中填写如下:
放电频次分布直方图如下,刺激时间将以一个黑框标示出来。
6.7放电光栅图
菜单:Plot->Rasters …:绘制放电的光栅分布图。参数设置对话框:
X Min. :横坐标(时间)最小值;
X Max.:横坐标最大值;
Units:选择需要进行绘图的单元编号。
Add stimulation markers:与频次直方图中相同。
6.8相关函数图
菜单:Plot->CrossCorrelograms…:绘制单元间的互相关和自相关函数。参数设置对话框:
Pre-ref. Time:横坐标(时间)最小值,计算相关函数时信号的偏移范围;
Post-ref. Time:横坐标最大值;
Number of Bins:相关函数分段数目;
Units:参与计算的单元编号;
Ref Units:作为参考的单元编号。
6.9刺激前后放电直方图PSTH
菜单:Plot->Peristimulus Histogram … :打开如下的参数设置对话框:
T. pre-Event:刺激事件前保留的时间长度;
T. post-Event:刺激事件后保留的时间长度;
Time/Bin:直方图的每段包含的时间长度;
Units:选择需要进行绘图的单元编号。
Stimulus Phase:选择刺激事件的相位,“on”表示刺激开始;“off”表示刺激结束;
Stimulus No.:当有多个刺激周期时,可以选择其中的一个刺激周期。设置为0时,将所有刺激周期的的PSTH叠加。
6.10刺激前后放电光栅图
菜单:Plot->Peristimulus Raster … :打开如下的参数设置对话框:
T. pre-Event:刺激事件前时间长度;
T. post-Event:刺激事件后时间长度;
Units:选择需要进行绘图的单元编号。
Stimulus Phase:选择刺激事件的相位,“on”表示刺激开始;“off”表示刺激结束;
Stimulus No.:当有多个刺激周期时,可以选择其中的一个刺激周期。设置为0时,将所有刺激前后放电都画出。
在该组信号中,共有三个刺激事件,Stimulus No.设置为0时,三个刺激前后的放电都绘出,如下图。纵坐标名称中,E*表示事件的编号,U*表示单元编号:
部分参考文献:
1Mukhopadhyay S, Ray GC. A new interpretation of nonlinear energy operator and its efficacy in spike detection. IEEE Transactions on Bio-Medical Engineering. 1998 Feb ;45(2):180-7.
2 (美)迪达等 著,李宏东等 译. 模式分类(原书第2版)[M]. 北京: 机械工业出版社,2003.
3 Zhang C, Zhang X, Zhang MQ, Li Y. Neighbor number, valley seeking and clustering. Pattern Recogn. Lett. 2007 ;28(2):173-180.
其它文献:
梁培基,陈爱华. 神经元活动的多电极同步记录及神经信息处理 [M]. 北京:北京工业出版社, 2003.
Lewicki MS. A review of methods for spike sorting: the detection and classification of neural action potential [J]. Net Comput Neural Syst, 1998, 9(4): 53-78.
Answer.com: Expectation-maximization algorithm