CN104458785B - 一种核磁共振波谱谱峰对齐及谱峰提取方法 - Google Patents

Abstract

本发明公开了一种核磁共振波谱谱峰对齐及谱峰提取方法，读取各个核磁共振谱图；利用具有内标化合物的谱峰或者预定化合物的特征谱峰对各个核磁共振谱图进行初步校准；将预定选取的核磁共振谱图作为参考谱图并划分为若干个固定区间，在各个固定区间内对待校准的核磁共振谱图进行校准；将每个核磁共振谱图在各个固定区间内校正后谱图进行整合。本发明适用性广几乎适应于所有的核磁共振谱图数据处理；数据处理简单；算法更快。

Description

一种核磁共振波谱谱峰对齐及谱峰提取方法

技术领域

本发明涉及一种图像处理方法，尤其是涉及到一种自动谱峰对齐方法，几乎适用于所有液体核磁共振波谱谱图位移校准，同时可以实现谱峰变量的快速提取。该方法适应于对一组平行的核磁共振波谱数据进行数据挖掘和信息提取等。

背景技术

理想情况下，核磁共振波谱仪采集到的波谱数据应该具有一致的化学位移以及稳定的基线。在实际分析过程中，由于受到待测样品的pH值、离子浓度，以及在采样过程中磁场的均一性和实验温度等客观因素的影响，谱峰经常会产生不可预测的漂移，而且很多因素是不可克服的。因此，谱峰位移校准是核磁共振波谱数据处理的关键步骤。现已有一些谱峰校准的方法，例如：相关优化解缠法(correlation optimized warping，COW)、多尺度谱峰对齐法(multi-scale peak alignment，MSPA)、高斯平滑谱峰对齐算法(Gaussiansmoothing，GPA)等，但这些算法的普适性有限，且操作复杂，不适用于一般研究人员。因此，我们研究出一种简单高效的谱峰对齐算法，可以实现核磁共振波谱谱峰基线自动批量对齐以及谱峰变量的快速提取等。

另外，高通量数据分析已成为当代科技发展的必然趋势，因此对平行的核磁共振波谱数据的变量提取也成为影响科研效率的关键性因素。在高通量核磁数据分析过程中，常用的两种变量提取方法：分段积分法和目标性分析法。其中分段积分法就是把核磁谱图分成若干等份，把每等分的积分面积作为统计分析变量，降低了大批样本分析对分析人员的要求，简化了数据。但是该方法，牺牲了数据的分辨率，且不一定存在生物学意义。目标性分析方法有效的解决了这个问题，该方法基于数据库比对，选择谱图内所有 的目标化合物作为分析变量，但是该方法需要分析工作者具有较强的核磁背景知识，且谱图分析过于耗时！因此，在谱峰对齐的基础上挑选谱峰将会作为一个不错的选择。因此，该发明同时提供一种快速谱峰识别及提取的方法。

该方法在一定程度上为核磁数据批量快速分析，减少核磁结果分析的工作量及对谱峰信号的深度挖掘和分析提供了有力的工具。

发明内容

本发明针对现有技术存在的上述问题，提供了一种简单易懂，快速高效的谱峰对齐方法，同时将核磁共振谱图所有谱峰进行集中提取，为后续统计学分析提供变量。该方法虽然是基于Matlab软件进行开发的，但是对于使用者的编程能力没有要求。当然，若使用者具有一定的编程能力，也可以对该方法进行深层次的挖掘。该发明几乎可以适用于所有的批量核磁共振谱图，且校正结果优于已发表的方法，结果准确，且具有一定的拓展性。

本发明主要是通过下述技术方案得以解决的：

一种核磁共振波谱谱峰对齐及谱峰提取方法，包括以下步骤：

步骤1、读取各个核磁共振谱图；

步骤2、利用具有内标化合物的谱峰或者预定化合物的特征谱峰对各个核磁共振谱图进行初步校准；

步骤3、将预定选取的核磁共振谱图作为参考谱图并划分为若干个固定区间，在各个固定区间内对待校准的核磁共振谱图进行校准；

步骤4、将每个核磁共振谱图在各个固定区间内校正后谱图进行整合。

如上所述的步骤3包括以下步骤：

步骤3.1、选取预定的核磁共振谱图作为参考谱图，将参考谱图划分为若干个固定区间，选定其中一个固定区间作为当前校准固定区间；

步骤3.2、求取当前校准固定区间内的谱峰最大值Y_max，同时确定谱峰最大值Y_max所对应的化学位移δ_ref，选取当前待校准的核磁共振谱图；

步骤3.3、预定搜索区间[δ_ref-Δδ’，δ_ref+Δδ’]，在当前待校准的核磁共振谱图的搜索区间[δ_ref-Δδ’，δ_ref+Δδ’]内寻找谱峰的所有峰值Y_p1,…,Y_pn，确定当前待校准的核磁共振谱图中与各个峰值相对应的化学位移，即δ_p1,…,δ_pn；pn为当前待校准的核磁共振谱图在搜索区间的谱峰个数；

步骤3.4、保持当前待校准的核磁共振谱图的谱形不变，将当前待校准的核磁共振谱图在搜索区间中的各个寻找到的峰值分别平移到化学位移δ_ref处，依次分别计算平移到δ_ref后的当前待校准的核磁共振谱图与参考谱图在当前校准固定区间内的相关系数R_px，其中px＝p1～pn；

步骤3.5、寻找在当前校准固定区间内相关系数R_px的最大的相关系数R_max，与R_max相对应的谱峰的化学位移为δ_max1，然后平移当前待校准的核磁共振谱图δ_max1-δ_ref化学位移，截取平移后的当前待校准的核磁共振谱图在当前校准固定区间内的全部谱峰信号；

步骤3.6、选取下一个待校准的核磁共振谱图的作为当前待校准的核磁共振谱图并返回步骤3.3直至所有的待校准的核磁共振谱图在当前校准固定区间内校准完毕；

步骤3.7、保存当前校准固定区间的所有核磁共振谱图；

步骤3.8、选取下一个固定区间作为当前校准固定区间并返回步骤3.2，直至所有固定区间校准完毕。

如上所述的步骤3.3中寻找谱峰的判断依据为：f(δ_pn)>f(δ_pn-1)，同时f(δ_pn)>f(δ_pn+1)，化学位移δ_pn对应的即为谱峰，f(δ_pn)、f(δ_pn-1)、f(δ_pn+1)分别表示化学位移δ_pn、δ_pn-1、δ_pn+1处的核磁共振波谱的信号强度。

一种核磁共振波谱谱峰对齐及谱峰提取方法，还包括谱峰变量提取步骤：

步骤5、确定噪声Y_noise；

步骤6、寻找参考谱图的谱峰的最高值Y_max1，

若Y_max1>Y_noise，则谱峰的最高值Y_max1作为谱峰变量，δ_max1为谱峰的最高值Y_max1所 对应的化学位移，否则终止挑选谱峰变量，进入步骤9；

步骤7：寻找谱峰区间，在δ_max1前后寻找波谷，前/后波谷之间的区间即为谱峰变量所对应的谱峰的区间[δ_max1A，δ_max1B]；

步骤8、将参考谱图在[δ_max1A，δ_max1B]内的数值去除，返回步骤6；

步骤9：挑选出参考谱图所有的谱峰变量后按照化学位移对谱峰变量进行排序，若谱峰变量对应的谱峰的区间的绝对值大于预定值，则将对应的谱峰变量删除。

如上所述的前/后波谷的搜索包括以下步骤：以化学位移δ_max1为中心向谱图两侧搜索谱值变大的前拐点和后拐点，前拐点和后拐点即分别为前波谷和后波谷。

本发明与现有技术相比，具有以下优点：

(1)适用性广：本方法对核磁共振谱图中的谱峰信号没有特殊要求，适应性更广，几乎适应于所有的核磁共振谱图数据处理；

(2)原理简单易懂：本发明在不改变谱峰性质的情况下，只对谱峰进行简单漂移处理，使用者不需要具有很强的核磁背景和编程能力。

(3)数据处理简单：处理结果直接导入Excel文件中，供研究人员使用；

(4)可拓展性：该算法是基于matlab程序进行编译的，使用者若具有一定的编程基础，即可以在此基础上对核磁共振谱图进行其他图形处理或数据分析。

(5)算法的优势更为明显：与现有算法相比，以文献报道(见图(1))的包括红、白和桃红三种葡萄酒共40个样品的¹H-NMR谱图片断(化学位移0.5-6.0)，每张谱图8712个数据点。针对于该组数据，三种文献报道方法所得到的平均相关系数为：COW(0.98)，MSPA(0.91)，GPA(0.98)，该算法则将谱峰相关系数提高至0.99。运行时间则由文献报道得最快62秒提高到1秒左右。。具体数据分析在后面举例中，进行详细说明。

附图说明

图1是本发明的流程图；

图2是红、白和桃红三种葡萄酒中乳酸的核磁共振信号校准前(A)后(B)对比图；

图3是红、白和桃红三种葡萄酒共40个样品的¹H-NMR谱图片断用目前算法对齐前后的相关谱图。(A)：对齐前；(B)：对齐后

具体实施方式

下面通过实例，并结合附图，对本发明的操作，作进一步具体的说明。

实施例1：

一种核磁共振波谱谱峰对齐及谱峰提取方法，包括以下步骤：

谱峰校准：

步骤1：读取各个核磁共振谱图，优选经过相位校正和基线调整后的核磁共振谱图；

步骤2：选取核磁共振谱图内标化合物(内标化合物为已知浓度和已知化学位移的添加物)的化学位移δ，对步骤1读取的批量谱图的化学位移进行初步校准。

选取方法一般选择具有内标化合物的谱峰的化学位移作为基准δ_s(例如TMSP-0ppm或者Formate-8.46ppm)。

在没有内标化合物的时候，需要选取预定化合物的特征谱峰的化学位移作为初次校准的基准δ_s(例如大脑组织样品提取物中的NAA-2.02ppm)。

校准方法如下所述：(A)：选择当前需校正的谱图在设定区间(δ_s-Δδ～δ_s+Δδ)的谱峰最大值Y_max＝max(f₁(δ_i)‥‥‥f₁(δ_j))，f₁(δ_i)表示在δ_i的谱峰强度。确定谱峰最大值Y_max所对应的化学位移为δ_max值，即δ_max＝f^-1(Y_max)。当前需校正的谱图的新的化学位移即为δ_new＝δ_old-(δ_max-δ_s)。δ_old为校正前谱图的化学位移。这样所有的核磁共振谱图完成了初次校准。

步骤3：固定区域内谱峰校准。测试样品内的成分较为复杂，受外界环境的影响具有不可预见性。因此，该方法中，实验者采取将预定的一个核磁共振谱图作为参考谱图，并将参考谱图划分为若干个固定区间，在各个固定区间内对待校准的核磁共振谱图进行校准；以谱峰相关系数R的大小对校准的效果进行评价。区间的划分原则以谱峰相关系 数R和基线调整结果为准。

两个核磁共振谱图的相关系数R的计算方法如下所示：

$R_{(a,b)}=\frac{\mathrm{Cov}(a,b)}{\sqrt{V\left(a\right)V\left(b\right)}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{a_i-\stackrel{\‾}{a}}{S_a}\·\frac{b_i-\stackrel{\‾}{b}}{S_b}$

式中，a_i和b_i分别为待比较的两个核磁共振谱图中化学位移δ_i所对应的谱峰强度，i＝1～n，n为待比较的核磁共振谱图中所包含数据点的个数。和分别为待比较的两个核磁共振谱图中的谱峰强度的平均值，S_a和S_b分别为待比较的两个核磁共振谱图的标准偏差。谱峰相关系数R已被广泛应用于评价谱峰校准的效果。

校准方法：由于样品之间受到pH值，温度以及样品内部微环境的差异，化学位移的变化以及谱峰的高低，会在一定程度上存在一定的差异，因此不能简单利用步骤2的方法进行谱峰漂移校准。

在谱峰漂移对齐的基础上，该方法又新引入了临近谱峰位移谱峰校准方法。具体包括以下步骤：

步骤3.1、选取预定的核磁共振谱图作为参考谱图，将参考谱图划分为若干个固定区间，选定其中一个固定区间作为当前校准固定区间[δ_A～δ_B]；

步骤3.2、求取当前校准固定区间[δ_A～δ_B]内的谱峰最大值Y_max，即Y_max＝max(f₁(δ_A)‥‥‥f₁(δ_B))，同时确定谱峰最大值Y_max所对应的δ_ref值，即δ_ref＝f^-1(Y_max)；选取当前待校准的核磁共振谱图；

步骤3.3、预定搜索区间[δ_ref-Δδ’，δ_ref+Δδ’]，即(δ_ref±Δδ’)。在当前待校准的核磁共振谱图的搜索区间[δ_ref-Δδ’，δ_ref+Δδ’]内寻找谱峰的所有峰值(Y_p1,…,Y_pn)，确定当前待校准的核磁共振谱图中与各个峰值相对应的化学位移，即δ_p1,…,δ_pn；pn为当前待校准的核磁共振谱图在搜索区间的谱峰个数。寻找峰值的判断依据为：f(δ_pn)>f(δ_pn-1)，同时f(δ_pn)>f(δ_pn+1)，即δ_pn即对应于一个峰值，依次类推。

步骤3.4、寻找到峰值后，保持当前待校准的核磁共振谱图的谱形不变，将当前待校准的核磁共振谱图在搜索区间中的各个寻找到的峰值分别平移到化学位移δ_ref处，即δ_new1＝δ_old-(δ_px-δ_ref)，其中px＝p1，…，pn，依次分别计算平移到δ_ref后的当前待校准的核磁共振谱图与参考谱图在当前校准固定区间[δ_A～δ_B]内的相关系数R_px。

步骤3.5、寻找在当前校准固定区间[δ_A～δ_B]内相关系数R_px的最大的相关系数R_max，即R_max＝max(R_px)。与R_max相对应的谱峰的化学位移为δ_max1，然后平移当前待校准的核磁共振谱图δ_max1-δ_ref化学位移，即δ_new1＝δ_old-(δ_max1-δ_ref)，截取平移后的待校准的核磁共振谱图[δ_A-δ_B]区间内的全部谱峰信号。δ_new1为平移后的当前待校准的核磁共振谱图的化学位移；δ_old为平移前的当前待校准的核磁共振谱图的化学位移；

步骤3.6、选取下一个待校准的核磁共振谱图的作为当前待校准的核磁共振谱图并返回步骤3.3直至所有的待校准的核磁共振谱图在当前校准固定区间内校准完毕；

步骤3.7、保存当前校准固定区间[δ_A～δ_B]的所有核磁共振谱图；

步骤3.8、选取下一个固定区间作为当前校准固定区间并返回步骤3.2，直至所有固定区间校准完毕。

步骤4、将每个核磁共振谱图在所有固定区间内校正后谱图进行整合，计算校正后的核磁共振谱图与参考谱图之间的相关系数R_j。j为1～m，m为核磁共振谱图的个数。计算各个相关系数R_j的平均值，即完成谱峰校准。

谱峰变量提取：

步骤5、确定噪音值，即Y_noise。

步骤6：寻找参考谱图的谱峰的最高值Y_max1，即Y_max1＝max(f(δ_i))，i＝1，…，N。N为参考谱图的谱峰个数，若Y_max1>Y_noise，即确定Y_max1所对应的化学位移δ_max1，即δ_max1＝f^-1(Y_max1)，其中max1为谱峰化学位移所对应的数据点位置，否则终止挑选谱峰，进入步骤9。

步骤7：寻找谱峰区间，在δ_max1前后寻找波谷，前后波谷之间的区间即为δ_max1所对应的谱峰的区间。

前/后波谷的搜索方法为：

若f(δ_max1-m)>f(δ_max1-m+1)，即停止搜索，则max1-m所对应的谱峰位置为前波谷，化学位移为δ_max1A；

若f(δ_max1+k)>f(δ_max1+k-1)，即停止搜索，则max1+k所对应的谱峰位置为后波谷，化学位移为δ_max1B。

其中，m，k为逐渐增加的数据点个数。

即寻找到δ_max1A到δ_max1B之间的谱峰，并对其进行积分。

步骤8、将参考谱图在[δ_max1A，δ_max1B]内的数值去除，返回步骤6；

步骤9：剔除噪音变量的干扰：由于该方法可能会引入一些只包含两三个点的谱峰区间，该区间对于实际样品分析不具有任何生物学意义。挑选出参考谱图所有的谱峰变量后按照化学位移对谱峰变量进行排序，若谱峰变量对应的max1A-max1B的绝对值大于预定值，则将[δ_max1A，δ_max1B]中的谱峰变量删除。即最终完成所有谱峰变量的提取和积分工作。

实施例2：

利用实施例1中记载的方法进行谱峰对齐及其谱峰提取方法，以文献报道的葡萄酒样品的核磁共振谱图为例，该样本已被用来评估多种谱峰对齐算法。样本中包含红、白和桃红三种葡萄酒共40个样品的¹H-NMR谱图片断(化学位移δ＝0.5～6.0)，每张谱图包含8172个数据点。

步骤1：从互联网上下载公共数据：http://www.models.life.ku.dk/wine_nmr；读入数据，数据中ppm变量为化学位移变量，X为核磁共振信号强度。

步骤2：所有葡萄酒样品的核磁共振谱图数据中均未加入内标化合物，因此需选取预定化合物的谱峰，作为初次校准的基准。以酒精的-CH₃基团在1.169位置的化学位移作为初次校准的标准，搜索设定区间Δδ为0.2ppm，对所有核磁共振谱图进行初次校准；

步骤3：逐级确定校准固定区间(表1)，同时对该区间内的所有谱峰进行平移校准。以校准区间(1.3598～1.4400)为例：由于受到葡萄酒pH值以及内部环境的影响，该区间 谱图的化学位移以及谱峰形状在所有核磁共振谱图中存在极为显著的差异，如图2A所示。在校准时，以0.03(Δδ’)为搜索区间，校准前后对比结果如图2所示(A：校准前，R：0.39；B：校准后：R：0.71)。其中校准前后各个区间的相关系数见表1；

表1：葡萄酒核磁共振谱图谱峰校准固定区间

步骤4、将每个核磁共振谱图在所有固定区间内校正后谱图进行整合，计算校正后的核磁共振谱图相关系数R_j的平均值为0.99，即完成谱峰校准。

步骤5：确定谱图噪音值：Y_noise＝50000；

步骤6：寻找参考谱图的谱峰的最高值Y_max1(例如首次挑选谱峰的强度为1.021*10⁹)，即Y_max1＝max(f(δ_i))，i＝1，…，n。若Y_max1>Y_noise，即确定Y_max1所对应的化学位移δ_max1(例如首次挑选谱峰所对应的化学位移为1.169)，即δ_max1＝f^-1(Y_max1)，其中max1为谱峰化学位移所对应的数据点位置(例如首次挑选谱峰所对应的数据点位置为7652)，否则终止挑选谱峰，进入步骤9。

步骤7：寻找谱峰区间，在δ_max1前后寻找波谷，前后波谷之间的区间即为δ_max1所对应的谱峰的区间，并对其进行积分。例如首次挑选谱峰所对应的前后波谷的化学位移分别为δ_max1A:1.1597；δ_max1B:1.1786。

步骤8：将参考谱图在[δ_max1A，δ_max1B]内的数值去除，返回步骤6；

步骤9：剔除噪音变量的干扰：由于该方法为全谱谱峰搜索，可能会引入一些只包含几个点的谱峰区间，该区间对于实际样品分析并不具有任何生物学意义，因此需要将这些噪音点剔除。挑选出参考谱图所有的谱峰变量后按照化学位移对谱峰变量进行排序，max1A-max1B的绝对值大于预定值(例如：0.002)，则将[δ_max1A，δ_max1B]中的谱峰变量删除。即最终完成所有谱峰变量的提取和积分工作。利用上述“谱峰变量提取”方法，获得变量1065个谱峰变量，剔出噪音干扰后，谱峰变量骤减为388个变量。而每一个变量即对应于一个核磁共振谱峰，即某一个化合物的特征峰或者几个化学物的重叠峰。该变量可以直接作为统计分析变量为大数据分析提供便利。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (4)

1.一种核磁共振波谱谱峰对齐及谱峰提取方法，包括以下步骤：

步骤1、读取各个核磁共振谱图；

步骤2、利用具有内标化合物的谱峰或者预定化合物的特征谱峰对各个核磁共振谱图进行初步校准；

步骤3、将预定选取的核磁共振谱图作为参考谱图并划分为若干个固定区间，在各个固定区间内对待校准的核磁共振谱图进行校准；

步骤4、将每个核磁共振谱图在各个固定区间内校正后谱图进行整合，

步骤5、确定噪声Y_noise；

步骤6、寻找参考谱图的谱峰的最高值Y_max1，

若Y_max1>Y_noise，则谱峰的最高值Y_max1作为谱峰变量，δ_max1为谱峰的最高值Y_max1所对应的化学位移，否则终止挑选谱峰变量，进入步骤9；

步骤7：寻找谱峰区间，在δ_max1前后寻找波谷，前/后波谷之间的区间即为谱峰变量所对应的谱峰的区间[δ_max1A，δ_max1B]；

步骤8、将参考谱图在[δ_max1A，δ_max1B]内的数值去除，返回步骤6；

步骤9：挑选出参考谱图所有的谱峰变量后按照化学位移对谱峰变量进行排序，若谱峰变量对应的谱峰的区间的绝对值大于预定值，则将对应的谱峰变量删除。

2.根据权利要求1所述的一种核磁共振波谱谱峰对齐及谱峰提取方法，其特征在于，所述的步骤3包括以下步骤：

步骤3.1、选取预定的核磁共振谱图作为参考谱图，将参考谱图划分为若干个固定区间，选定其中一个固定区间作为当前校准固定区间；

步骤3.2、求取当前校准固定区间内的谱峰最大值Y_max，同时确定谱峰最大值Y_max所对应的化学位移δ_ref，选取当前待校准的核磁共振谱图；

步骤3.3、预定搜索区间[δ_ref-Δδ’，δ_ref+Δδ’]，在当前待校准的核磁共振谱图的搜索区间[δ_ref-Δδ’，δ_ref+Δδ’]内寻找谱峰的所有峰值Y_p1,…,Y_pn，确定当前待校准的核磁共振谱图中与各个峰值相对应的化学位移，即δ_p1,…,δ_pn；pn为当前待校准的核磁共振谱图在搜索区间的谱峰个数；

步骤3.4、保持当前待校准的核磁共振谱图的谱形不变，将当前待校准的核磁共振谱图在搜索区间中的各个寻找到的峰值分别平移到化学位移δ_ref处，依次分别计算平移到δ_ref后的当前待校准的核磁共振谱图与参考谱图在当前校准固定区间内的相关系数R_px，其中px＝p1～pn；

步骤3.5、寻找在当前校准固定区间内相关系数R_px的最大的相关系数R_max，与R_max相对应的谱峰的化学位移为δ_max，然后平移当前待校准的核磁共振谱图δ_max-δ_ref化学位移，截取平移后的当前待校准的核磁共振谱图在当前校准固定区间内的全部谱峰信号；

步骤3.6、选取下一个待校准的核磁共振谱图的作为当前待校准的核磁共振谱图并返回步骤3.3直至所有的待校准的核磁共振谱图在当前校准固定区间内校准完毕；

步骤3.7、保存当前校准固定区间的所有核磁共振谱图；

步骤3.8、选取下一个固定区间作为当前校准固定区间并返回步骤3.2，直至所有固定区间校准完毕。

3.根据权利要求2所述的一种核磁共振波谱谱峰对齐及谱峰提取方法，其特征在于，所述的步骤3.3中寻找谱峰的判断依据为：f(δ_pn)>f(δ_pn-1)，同时f(δ_pn)>f(δ_pn+1)，化学位移δ_pn对应的即为谱峰，f(δ_pn)、f(δ_pn-1)、f(δ_pn+1)分别表示化学位移δ_pn、δ_pn-1、δ_pn+1处的核磁共振波谱的信号强度。

4.根据权利要求1所述的一种核磁共振波谱谱峰对齐及谱峰提取方法，其特征在于，所述的前/后波谷的寻找包括以下步骤：以化学位移δ_max1为中心向谱图两侧搜索谱值变大的前拐点和后拐点，前拐点和后拐点即分别为前波谷和后波谷。