2008年9月24日星期三

SourceForge建立项目

前几天,在Sourceforge 上建了一个项目。平时网络的“花边”新闻倒是看了不少,一操作起来可就立马撞墙了。SF的帮助在它的Wiki 上,help里面反倒没什么有用的信息。简单总结一下,供外行们参考。

Sourceforge支持OpenID登录,但必须要申请一个SF的ID才可以建立项目。
项目的名称是以后可以改的。但unix name不能改,必须是唯一的。
项目必须经过站方的审批。所以关于项目的描述应该是英文的(这是我猜的,SF的审核人员应该不懂中文吧)。审核需要一个工作日,通过后需要为自己的项目建立分类信息等等,这些也都可以在以后重新修改。

审核通过后就可以发布源代码等文件了。在Admin中选择File Release,建立一个新的Package,新的Release。仔细看一下里面的文字描述,一般有不少链接,指向更详细的说明。
可以通过很多方式将本地文件上传到SF服务器,Web直接上传或者使用一些工具。只有 rsync over ssh是支持断点续传的,推荐国内用户使用。Windows用户可以安装Cygwin ,它相当于Windows上的Linux虚拟机。Cygwin是在线安装的,选择一个镜像地址,在安装包里选择以rsync和openSSH开头的几个包。完成后运行,会在安装目录里生成用户文件夹,例如当前Windows用户是USER,安装目录中会有homeUSER文件夹。将需要上传的文件拷贝到这个目录,然后按照SourceForge.net Wiki上的描述,在Cygwin中执行相应的命令就可以了。
通过Rsync over SSH上传文件
[jsmith@linux ~]# rsync -avP -e ssh FILE AcountName@frs.sourceforge.net:uploads/
AcountName@frs.sourceforge.net's password:
building file list ...
1 file to consider
FILE
15000000 100% 34.13kB/s 0:07:08 (xfer#1, to-check=0/1)

sent 15001925 bytes received 42 bytes 29560.53 bytes/sec
total size is 15000000 speedup is 1.00

等到最后一行出现,才表示上传完成,有时候需要点耐心。但总比网页上传一次次报错好多了。然后到文件发布页面,选择修改Release,里面应该已经有了上传的文件,附加上去就可以了。

还可以在SourceForge上为项目建立网站,SourceForge提供了100MB的空间。对静态网页,在本地将网站建好,上传上去就OK了。我只建了一个index.html,只需要上传到服务器上的htdocs目录下就可以了。

用Cygwin中的sftp上传

 sftp AcountName,ProjectName@web.sourceforge.net
cd htdocs
put index.html
bye

然后在浏览器中输入http://ProjectName.sourceforge.net就可以访问项目主页了。

2008年9月4日星期四

SPKtool使用手册

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对读入的数据进行低通滤波,使用6Butterworth滤波器,截至频率4500Hz。可以打开MLoadRaw.m,找到butter函数对应的位置,修改滤波参数。

2.2锋电位波形数据

从原始信号中检测出的锋电位波形数据,至少应该包含以下三个变量:

posind 锋电位波形在连续数据中的采样位置标记。行向量。第n个波形发生的时刻为posind(n)/fs

waveforms 锋电位波形数据,二维数组。每行对应一个锋电位的采样数据。例如,1000个波形,每个波形含40个数据点, 则waveforms尺度为1000×40

fs 采样频率。

数据文件中还可以包含如下一些变量:

featurespace 锋电位特征空间,二维数组,每行包含一个波形的特征。如1000个波形的PC1PC2特征,则featurespace1000×2。特征的数目不限。

featnames:特征名称,包含字符串的单元数组,每个字符串对应一个特征的名称。例如: 两个特征的名称分别为PC1PC2, 则 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。如果采样频率为40KHz1ms长度的锋电位波形包含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 Parameters0.5,则Slice特征为锋电位中第32×0.5=16个采样点的数值。Slice Parameters可以包含多个值,提取多个采样点,例如 0.50.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.pcaPCA作为特征时,由于共用PCA参数,可以通过该工具重新选择。

4.3特征分布散点图

Features->Plot Features 绘制特征分布散点图,如果定义了两个特征,散点图是二维的,如果有三个或更多的特征,则只使用前三个,以三维形式显示,并同时画出特征向三个平面上的投影分布,如下图。

4.4特征密度图

Features-> Features Density:特征分布密度图。如下图,该图只能绘出前两个特征的分布密度。

4.5特征对比图

Features-> Features Vs Features,绘制特征分布的对比图,注意该功能中使用的特征与聚类分析的特征是独立的。并且,每次绘图都需要重新计算特征。

如果已经划分类别,不同类别的特征将以不同颜色区别显示。

4.6主成份

Features-> Principal Components:给出特征提取中主成份曲线。例如特征定义中使用第一个和第二个主成份,则给出PC1PC2的形状曲线。

5 聚类分析

聚类功能集中于菜单Sort中。在选择具体的聚类方法之前,应该先设置聚类的参数。

5.1参数设置

菜单Sort->Options…, 打开参数设置对话框:

其中共有四个参数:

Unit Number:划分类别的数目,K-Means,期望最大化(EM)方法,以及模板匹配方法使用该参数;

K-Means RadiusK-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中,然后按照上面的方法,先移除没有归类的边缘特征,再将单元10(未分类)合并。

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. IntervalMax 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