对于两个序列来说，如果要比较两个波形的相似程度，可以使用DWT（动态时间规整）的方法。对于dwt方法，可以解决处理两个长度不一样的序列。
 
DTW是一种衡量两个时间序列之间的相似度的方法，主要应用在语音识别领域来识别两段语音是否表示同一个单词，度量的特征量是：两个序列之间的最短距离。
 
原理：
 
DTW通过把时间序列进行延伸和缩短，来计算两个时间序列性之间的相似性。
 
1，两个要进行匹配的数据A=[A1，A2，...An]和B=[B1，B2，...Bm]
 
归整路径的形式为W=w1,w2,...,wK，其中Max(|A|,|B|)<=K<=|A|+|B|。
 
wk的形式为(i,j)，其中i表示的是A中的i坐标，j表示的是B中的j坐标.
 
归整路径W必须从w1=(1,1)开始，到wK=(|X|,|Y|)结尾，以保证X和Y中的每个坐标都在W中出现.
 
(i,j)必须是单调增加的，从w(a1,a2)沿着某条路径到达w(am,bn)。
 
找路径满足：假如当前节点是w(ai,bj)，那么下一个节点必须是在   w(i+1,j)，w(i,j+1)，w(i+1,j+1)之间选择，并且路径必须是最短的。
 
 
计算的时候是按照动态规划的思想计算，也就是说在计算到达第(i,j)个节点的最短路径时候，考虑的是左下角也即第(i-1,j)、(i-1,j-1)、(i,j-1)这三个点到(i,j)的最短距离。
 
最后的目标是要找到一个在两个序列之间的最短距离以及实现这个最短距离的路径。
 
距离选用任意经典的距离计算方法：欧几里得距离
 
2，代码实现-matlab
 
 
function [Dist,D,k,w,rw,tw]=dtw(r,t,pflag)
%
% [Dist,D,k,w,rw,tw]=dtw(r,t,pflag)
%
% Dynamic Time Warping Algorithm
% Dist is unnormalized distance between t and r
% D is the accumulated distance matrix
% k is the normalizing factor
% w is the optimal path
% t is the vector you are testing against
% r is the vector you are testing
% rw is the warped r vector
% tw is the warped t vector
% pflag  plot flag: 1 (yes), 0(no)
%
% Version comments:
% rw, tw and pflag added by Pau Mic

[row,M]=size(r); if (row > M) M=row; r=r'; end;
[row,N]=size(t); if (row > N) N=row; t=t'; end;
d=sqrt((repmat(r',1,N)-repmat(t,M,1)).^2); %this makes clear the above instruction Thanks Pau Mic
dd=abs(repmat(r',1,N)-repmat(t,M,1));
dd
d
D=zeros(size(d));
D(1,1)=d(1,1);

for m=2:M
    D(m,1)=d(m,1)+D(m-1,1);
end
for n=2:N
    D(1,n)=d(1,n)+D(1,n-1);
end
for m=2:M
    for n=2:N
        D(m,n)=d(m,n)+min(D(m-1,n),min(D(m-1,n-1),D(m,n-1))); % this double MIn construction improves in 10-fold the Speed-up. Thanks Sven Mensing
    end
end

Dist=D(M,N);
n=N;
m=M;
k=1;
w=[M N];
w
while ((n+m)~=2)
    if (n-1)==0
        m=m-1;
    elseif (m-1)==0
        n=n-1;
    else
        [values,number]=min([D(m-1,n),D(m,n-1),D(m-1,n-1)]);
        switch number
            case 1
                m=m-1;
            case 2
                n=n-1;
            case 3
                m=m-1;
                n=n-1;
        end
    end
    k=k+1;
    w=[m n; w]; % this replace the above sentence. Thanks Pau Mic
end
[values,number]
D
m
n
w
% warped waves
rw=r(w(:,1));
tw=t(w(:,2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if pflag
    
    % --- Accumulated distance matrix and optimal path
    figure('Name','DTW - Accumulated distance matrix and optimal path', 'NumberTitle','off');
    
    main1=subplot('position',[0.19 0.19 0.67 0.79]);
    image(D);
    cmap = contrast(D);
    colormap(cmap); % 'copper' 'bone', 'gray' imagesc(D);
    hold on;
    x=w(:,1); y=w(:,2);
    ind=find(x==1); x(ind)=1+0.2;
    ind=find(x==M); x(ind)=M-0.2;
    ind=find(y==1); y(ind)=1+0.2;
    ind=find(y==N); y(ind)=N-0.2;
    plot(y,x,'-w', 'LineWidth',1);
    hold off;
    axis([1 N 1 M]);
    set(main1, 'FontSize',7, 'XTickLabel','', 'YTickLabel','');
    
    colorb1=subplot('position',[0.88 0.19 0.05 0.79]);
    nticks=8;
    ticks=floor(1:(size(cmap,1)-1)/(nticks-1):size(cmap,1));
    mx=max(max(D));
    mn=min(min(D));
    ticklabels=floor(mn:(mx-mn)/(nticks-1):mx);
    colorbar(colorb1);
    set(colorb1, 'FontSize',7, 'YTick',ticks, 'YTickLabel',ticklabels);
    set(get(colorb1,'YLabel'), 'String','Distance', 'Rotation',-90, 'FontSize',7, 'VerticalAlignment','bottom');
    
    left1=subplot('position',[0.07 0.19 0.10 0.79]);
    plot(r,M:-1:1,'-b');
    set(left1, 'YTick',mod(M,10):10:M, 'YTickLabel',10*rem(M,10):-10:0)
    axis([min(r) 1.1*max(r) 1 M]);
    set(left1, 'FontSize',7);
    set(get(left1,'YLabel'), 'String','Samples', 'FontSize',7, 'Rotation',-90, 'VerticalAlignment','cap');
    set(get(left1,'XLabel'), 'String','Amp', 'FontSize',6, 'VerticalAlignment','cap');
    
    bottom1=subplot('position',[0.19 0.07 0.67 0.10]);
    plot(t,'-r');
    axis([1 N min(t) 1.1*max(t)]);
    set(bottom1, 'FontSize',7, 'YAxisLocation','right');
    set(get(bottom1,'XLabel'), 'String','Samples', 'FontSize',7, 'VerticalAlignment','middle');
    set(get(bottom1,'YLabel'), 'String','Amp', 'Rotation',-90, 'FontSize',6, 'VerticalAlignment','bottom');
    
    % --- Warped signals
    figure('Name','DTW - warped signals', 'NumberTitle','off');
    
    subplot(1,2,1);
    set(gca, 'FontSize',7);
    hold on;
    plot(r,'-bx');
    plot(t,':r.');
    hold off;
    axis([1 max(M,N) min(min(r),min(t)) 1.1*max(max(r),max(t))]);
    grid;
    legend('signal 1','signal 2');
    title('Original signals');
    xlabel('Samples');
    ylabel('Amplitude');
    
    subplot(1,2,2);
    set(gca, 'FontSize',7);
    hold on;
    plot(rw,'-bx');
    plot(tw,':r.');
    hold off;
    axis([1 k min(min([rw; tw])) 1.1*max(max([rw; tw]))]);
    grid;
    legend('signal 1','signal 2');
    title('Warped signals');
    xlabel('Samples');
    ylabel('Amplitude');
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% test
clear  
clc  
a=[8 9 1 9 6 1 3 5]';  
b=[2 5 4 6 7 8 3 7 7 2]';  
[Dist,D,k,w,rw,tw] = DTW(a,b,1);  
fprintf('最短距离为%d\n',Dist)  
fprintf('最优路径为')  
w


 

 
 
 
3，代码实现-python code1--画图
 
from math import *
import matplotlib.pyplot as plt
import numpy

def print_matrix(mat) :
    print '[matrix] width : %d height : %d' % (len(mat[0]), len(mat))
    print '-----------------------------------'
    for i in range(len(mat)) :
        print mat[i]#[v[:2] for v in mat[i]]

def dist_for_float(p1, p2) :
    dist = 0.0
    elem_type = type(p1)
    if  elem_type == float or elem_type == int :
        dist = float(abs(p1 - p2))
    else :
        sumval = 0.0
        for i in range(len(p1)) :
            sumval += pow(p1[i] - p2[i], 2)
        dist = pow(sumval, 0.5)
    return dist

def dtw(s1, s2, dist_func) :
    w = len(s1)
    h = len(s2)
    
    mat = [([[0, 0, 0, 0] for j in range(w)]) for i in range(h)]
    
    #print_matrix(mat)
    
    for x in range(w) :
        for y in range(h) :            
            dist = dist_func(s1[x], s2[y])
            mat[y][x] = [dist, 0, 0, 0]
            
    #print_matrix(mat)

    elem_0_0 = mat[0][0]
    elem_0_0[1] = elem_0_0[0] * 2

    for x in range(1, w) :
        mat[0][x][1] = mat[0][x][0] + mat[0][x - 1][1]
        mat[0][x][2] = x - 1
        mat[0][x][3] = 0

    for y in range(1, h) :
        mat[y][0][1] = mat[y][0][0] + mat[y - 1][0][1]            
        mat[y][0][2] = 0
        mat[y][0][3] = y - 1

    for y in range(1, h) :
        for x in range(1, w) :
            distlist = [mat[y][x - 1][1], mat[y - 1][x][1], 2 * mat[y - 1][x - 1][1]]
            mindist = min(distlist)
            idx = distlist.index(mindist)
            mat[y][x][1] = mat[y][x][0] + mindist
            if idx == 0 :
                mat[y][x][2] = x - 1
                mat[y][x][3] = y
            elif idx == 1 :
                mat[y][x][2] = x
                mat[y][x][3] = y - 1
            else :
                mat[y][x][2] = x - 1
                mat[y][x][3] = y - 1

    result = mat[h - 1][w - 1]
    retval = result[1]
    path = [(w - 1, h - 1)]
    while True :
        x = result[2]
        y = result[3]
        path.append((x, y))

        result = mat[y][x]
        if x == 0 and y == 0 :
            break
        
    #print_matrix(mat)

    return retval, sorted(path)

def display(s1, s2) :

    val, path = dtw(s1, s2, dist_for_float)
    
    w = len(s1)
    h = len(s2)
    
    mat = [[1] * w for i in range(h)]
    for node in path :
        x, y = node
        mat[y][x] = 0

    mat = numpy.array(mat)
    
    plt.subplot(2, 2, 2)
    c = plt.pcolor(mat, edgecolors='k', linewidths=4)
    plt.title('Dynamic Time Warping (%f)' % val)

    plt.subplot(2, 2, 1)
    plt.plot(s2, range(len(s2)), 'g')
    
    plt.subplot(2, 2, 4)
    plt.plot(range(len(s1)), s1, 'r')
        

    plt.show()
    

s1 = [1, 2, 3, 4, 5, 5, 5, 4]
s2 = [3, 4, 5, 5, 5, 4]
s2 = [1, 2, 3, 4, 5, 5]
s2 = [2, 3, 4, 5, 5, 5]
#val, path = dtw(s1, s2, dist_for_float)
display(s1, s2)

 
 
 
4，代码实现-python  code2-只计算距离
 
 
import numpy as np
from numpy import *
from numpy.matlib import repmat


def DTW(r, t, plt=True):
    M = len(r)
    N = len(t)
    r_array = np.matrix(r)
    t_array = np.matrix(t)
    ### 距离计算公式
    distance_classical = abs(repmat(r_array.T, 1, N) - repmat(t_array, M, 1))
    D = np.zeros(shape(distance_classical))
    for m in range(1, M):
        D[m, 0] = distance_classical[m, 0] + D[m - 1, 0]
    for n in range(1, N):
        D[0, n] = distance_classical[0, n] + D[0, n - 1]
    for m in range(1, M):
        for n in range(1, N):
            D[m, n] = distance_classical[m, n] + min(D[m - 1, n], min(D[m - 1, n - 1], D[m, n - 1]))
    Dist = D[m, n]

 

 

 
 

 

序列的相似性可以是定量的数值，也可以是定性的描述。相似度是一个数值，反映两条序列的相似程度。关于两条序列之间的关系，有许多名词，如相同、相似、同源、同功、直向同源、共生同源等。在进行序列比较时经常使用“同源”（homology）和“相似”（similarity）这两个概念，这是两个经常容易被混淆的不同概念。两条序列同源是指它们具有共同的祖先。在这个意义上，无所谓同源的程度，两条序列要么同源，要么不同源。而相似则是有程度的差别，如两条序列的相似程度达到30%或60%。一般来说，相似性很高的两条序列往往具有同源关系。但也有例外，即两条序列的相似性很高，但它们可能并不是同源序列，这两条序列的相似性可能是由随机因素所产生的，这在进化上称为“趋同”（convergence），这样一对序列可称为同功序列。直向同源（orthologous）序列是来自于不同的种属同源序列，而共生同源（paralogous）序列则是来自于同一种属的序列，它是由进化过程中的序列复制而产生的。
 
 
序列比较的基本操作是比对（align）。两条序列的比对（alignment）是指这两条序列中各个字符的一种一一对应关系，或字符对比排列。序列的比对是一种关于序列相似性的定性描述，它反映在什么部位两条序列相似，在什么部位两条序列存在差别。最优比对揭示两条序列的最大相似程度，指出序列之间的根本差异。
 
字母表和序列
 
在生物分子信息处理过程中，将生物分子序列抽象为字符串，其中的字符取自特定的字母表。字母表是一组符号或字符，字母表中的元素组成序列。一些重要的字母表有：
 
（1）4字符DNA字母表 {A, C, G, T}；
 
（2）扩展的遗传学字母表或IUPAC编码；
 
（3）单字母氨基酸编码；
 
（4）上述字母表形成的子集。
 
下面所讨论的内容独立于特定的字母表。 首先规定一些特定的符号：
 
① A — 字母表；
 
② A* — 由字母表A中字符所形成的一系列有限长度序列或字符串的集合；
 
③ a、b、c — 单独的字符；
 
④ s、t、u、v、x — A*中的序列；
 
⑤ |s| — 序列s的长度。
 
为了说明序列s的子序列和s中单个字符，我们在s中各字符之间用数字标明分割边界。例如，设s=ACCACGTA，则s可表示为
 
0A1C2C3A4C5G6T7A8 。
 
i:s:j 指明第i位或第j位之间的子序列。当然，0 £ i £ j £ |s|。子序列0 : s : i 称为前缀，即prefix(s,i)，而子序列 i:s:|s| 称为后缀suffix(s, |s|-i+1)。有两种特殊的情况，即 i=j或i = j-1。
 
① i:s:i 表示空序列
 
② ( j-1):s: j 表示s 中的第j 个字符，简记为sj 。
 
一般认为，子序列与计算机算法中子串的概念相当。但是，严格地讲，子序列与子串的概念是有区别的：子串是子序列，而子序列不一定是子串。可以通过选取s中的某些字符（或删除s中的某些字符）而形成s的子序列，例如TTT是ATATAT的子序列。而s的子串则是由s中相继的字符所组成，例如TAC是AGTACA的子串，但不是TTGAC的子串。如果t是s的子串，则称s是t的超串。子串也可以称为连续子序列。
 
两条序列s和t的连接用s + + t来表示，如：
 
ACC++CTA = ACCCTA
 
字符串操作除连接操作之外，另有一个k操作，即删除一个字符串两端的字符。其定义如下：
 
prefix(s,l) = sk|s|-l ，
 
suffix(s,l) = k|s|-ls ，
 
i:s:j = ki-1sk|s|-j 。
 
序列比较可以分为四种基本情况，具体任务和应用说明如下：
 
（1）假设有两条长度相近的、来自同一个字母表的序列，它们之间非常相似，仅仅是有一些细微的差别，例如字符的插入、字符的删除和字符替换，要求找出这两条序列的差别。这种操作实际应用比较多，例如，有两个实验室同时测定某个基因的DNA序列，其结果可能不一样，需要通过序列比较来比较实验结果。
 
（2）假设有两条序列，要求判断是否有一条序列的前缀与另一条序列的后缀相似，如果是，则分别取出前缀和后缀。该操作常用于大规模DNA测序中序列片段的组装。
 
 
 
（3）假设有两条序列，要求判断其中的一条序列是否是另一条序列的子序列。这种操作常用于搜索特定的序列模式。
 
（4）假设有两条序列，要求判断这两条序列中是否有非常相似的子序列。这种操作可用于分析保守序列。
 
当然，进行序列比较时，往往还需要说明是采取全局比较，还是采取局部比较。全局比较是比较两条完整的序列，而局部比较是找出最大相似的子序列。
 
编辑距离（Edit Distance）
 
观察这样两条DNA序列：GCATGACGAATCAG和TATGACAAACAGC。一眼看上去，这两条序列并没有什么相似之处，然而如果将第二条序列错移一位，并对比排列起来以后，就可以发现它们的相似性。
 
 
如果进一步在第二条序列中加上一条短横线，就会发现原来这两条序列有更多的相似之处。
 
 
上面是两条序列相似性的一种定性表示方法，为了说明两条序列的相似程度，还需要定量计算。有两种方法可用于量化两条序列的相似程度：一为相似度，它是两条序列的函数，其值越大，表示两条序列越相似；与相似度对应的另一个概念是两条序列之间的距离，距离越大，则两条序列的相似度就越小。在大多数情况下，相似度和距离可以交互使用，并且距离越大，相似度越小，反之亦然。但一般而言，相似度使用得较多，并且灵活多变。
 
最简单的距离就是海明（Hamming）距离。对于两条长度相等的序列，海明距离等于对应位置字符不同的个数。例如，下图是3组序列海明距离的计算结果。
 
 
使用距离来计算不够灵活，这是因为序列可能具有不同的长度，两条序列中各位置上的字符并不一定是真正的对应关系。例如，在DNA复制的过程中，可能会发生像删除或插入一个碱基这样的错误，虽然两条序列的其他部分相同，但由于位置的移动导致海明距离的失真。就图3.1中例子最右边的情况，海明距离为6，简单地从海明距离来看，两条序列差别很大（整个序列的长度只有8bp），但是，如果从s中删除G，从t中删除T，则两条序列都成为ACACACA，这说明两条序列仅仅相差两个字符。实际上，在许多情况下，直接运用海明距离来衡量两条序列的相似程度是不合理的。
 
为了解决字符插入和删除问题，引入字符“编辑操作”（Edit Operation）的概念，通过编辑操作将一个序列转化为一个新序列。用一个新的字符“-”代表空位（或空缺，Space），并定义下述字符编辑操作：
 
Match（a，a） — 字符匹配；
 
Delete（a，-） — 从第一条序列删除一个字符，或在第二条序列相应的位置插入空白字符；
 
Replace（a，b） — 以第二条序列中的字符b替换第一条序列中的字符a，a¹b；
 
Insert（-，b） — 在第一条序列插入空位字符，或删除第二条序列中的对应字符b。
 
很显然，在比较两条序列s和t时，在s中的一个删除操作等价于在t中对应位置上的一个插入操作，反之亦然。需要注意的是，两个空位字符不能匹配，因为这样的操作没有意义。引入上述编辑操作后，重新计算两条序列的距离，就成为编辑距离。
 
以上的操作仅仅是关于序列的常用操作，在实际应用中还可以引入复杂的序列操作。下面是两条序列的一种比对：
 
 
上述比对不能反映两条序列的本质关系。但是，如果将第二条序列头尾倒置，可以发现两条序列惊人的相似：
 
 
再比如，下面两条序列有什么关系？如果将其中一条序列中的碱基替换为其互补碱基，就会发现其中的关系：
 
 
 
通过点矩阵分析两条序列的相似之处
 
进行序列比较的一个简单的方法是“矩阵作图法”或“对角线作图”，这种方法是由Gibb首先提出的。将两条待比较的序列分别放在矩阵的两个轴上，一条在X轴上，从左到右，一条在Y轴上，从下往上，如图3.2所示。当对应的行与列的序列字符匹配时，则在矩阵对应的位置作出“点”标记。逐个比较所有的字符对，最终形成点矩阵。
 
 
 
 
序列比较矩阵标记图
显然，如果两条序列完全相同，则在点矩阵主对角线的位置都有标记；如果两条序列存在相同的子串，则对于每一个相同的子串对，有一条与对角线平行的由标记点所组成的斜线，如图a中的斜线代表相同的子串“ATCC”；而对于两条互为反向的序列，则在反对角线方向上有标记点组成的斜线，如图b所示。
 
 
 
 
a:相同子串矩阵标记图
 
 
 
b:反向序列矩阵标记图
对于矩阵标记图中非重叠的与对角线平行斜线，可以组合起来，形成两条序列的一种比对。在两条子序列的中间可以插入符号“-”，表示插入空位字符。在这种对比之下分析两条序列的相似性，如下图所示。找两条序列的最佳比对（对应位置等同字符最多），实际上就是在矩阵标记图中找非重叠平行斜线最长的组合。
 
 
除非已经知道待比较的序列非常相似，一般先用点矩阵方法比较，因为这种方法可以通过观察矩阵的对角线迅速发现可能的序列比对。
 
实例一：
 
 
实例二：
 
 
两条序列中有很多匹配的字符对，因而在点矩阵中会形成很多点标记。当对比较长的序列进行比较时，这样的点阵图很快会变得非常复杂和模糊。使用滑动窗口代替一次一个位点的比较是解决这个问题的有效方法。假设窗口大小为10，相似度阈值为8。首先，将X轴序列的第1-10个字符与Y轴序列的第1-10个字符进行比较。如果在第一次比较中，这10个字符中有8个或者8个以上相同，那么就在点阵空间（1,1）的位置画上点标记。然后窗口沿X轴向右移动一个字符的位置，比较X轴序列的第2-11个字符与Y轴序列的第1-10个字符。不断重复这个过程，直到X轴上所有长度为10的子串都与Y轴第1-10个字符组成的子串比较过为止。然后，将Y轴的窗口向上移动一个字符的位置，重复以上过程，直到两条序列中所有长度为10的子串都被两两比较过为止。基于滑动窗口的点矩阵方法可以明显地降低点阵图的噪声，并且可以明确地指出两条序列间具有显著相似性的区域。
 
 
序列的两两比对
 
序列的两两比对（Pairwise Sequence Alignment）就是对两条序列进行编辑操作，通过字符匹配和替换，或者插入和删除字符，使得两条序列达到一样的长度，并使两条序列中相同的字符尽可能地一一对应。设两条序列分别是s和t，在s或t中插入空位符号，使s和t达到一样的长度。下图是对序列AGCACACA和ACACACTA的两种比对结果以及对应的字符编辑操作。
 
 
下面就不同类型的编辑操作定义函数w，它表示“代价（cost）”或“权重（weight）”。对字母表A中的任意字符a、b，定义：
 
 
这是一种简单的代价定义，在实际应用中还需使用更复杂的代价模型。一方面，可以改变各编辑操作的代价值，例如，在蛋白质序列比较时，用理化性质相近的氨基酸进行替换的代价应该比完全不同的氨基酸替换代价小；另一方面，也可以使用得分（score）函数来评价编辑操作。下面给出一种基本的得分函数：
 
 
在进行序列比对时，可根据实际情况选用代价函数或得分函数，即选用（3-1）式或（3-2）式。
 
下面给出在进行序列比对时常用的概念：
 
（1）两条序列s 和 t 的比对的得分（或代价）等于将s 转化为t 所用的所有编辑操作的得分（或代价）总和；
 
（2）s 和t 的最优比对是所有可能的比对中得分最高（或代价最小）的一个比对；
 
（3）s 和t 的真实距离应该是在得分函数p值（或代价函数w值）最优时的距离。
 
使用前面代价函数w的定义，可以得到下列比对的代价：
 
s: AGCACAC-A
 
t: A-CACACTA
 
cost（s，t）= 2
 
而使用得分函数p 的定义，可以得到下列比对的得分：
 
s: AGCACAC-A
 
t: A-CACACTA
 
score (s，t)= 5
 
进行序列比对的目的是寻找一个得分最高（或代价最小）的比对。
 
用于序列相似性的打分矩阵（scoring matrix）
 
无论是3-1式还是3-2式，都是简单相似性评价模型，在计算比对的代价或得分时，对字符替换操作只进行统一的处理，没有考虑“同类字符”替换与“非同类字符”替换的差别。实际上，不同类型的字符替换，其代价或得分是不一样的，特别是对于蛋白质序列。某些氨基酸可以很容易地相互取代而不用改变它们的理化性质。例如，考虑这样两条蛋白质序列，其中一条在某一位置上是丙氨酸，如果该位点被替换成另一个较小且疏水的氨基酸，比如缬氨酸，那么对蛋白质功能的影响可能较小；如果被替换成较大且带电的残基，比如赖氨酸，那么对蛋白功能的影响可能就要比前者大。直观地讲，比较保守的替换比起较随机替换更可能维持蛋白质的功能，且更不容易被淘汰。因此，在为比对打分时，我们可能更倾向对丙氨酸与缬氨酸的比对位点多些奖励，而对于丙氨酸与那些大而带电氨基酸（比如赖氨酸）的比对位点则相反。理化性质相近的氨基酸残基之间替换的代价显然应该比理化性质相差甚远的氨基酸残基替换得分高，或者代价小。同样，保守的氨基酸替换得分应该高于非保守的氨基酸替换。这样的打分方法在比对非常相近的序列以及差异极大的序列时，会得出不同的分值。这就是提出打分矩阵（或者称为取代矩阵）的原由。在打分矩阵中，详细地列出各种字符替换的得分，从而使得计算序列之间的相似度更为合理。在比较蛋白质时，我们可以用打分矩阵来增强序列比对的敏感性。打分矩阵是序列比较的基础，选择不同的打分矩阵将得到不同的比较结果，而了解打分矩阵的理论依据将有助于在实际应用中选择合适的打分矩阵。以下介绍一些常用的打分矩阵或代价矩阵。
 
1、核酸打分矩阵
 
设核酸序列所用的字母表为 A = { A，C，G，T }。
 
（1）等价矩阵
 
等价矩阵（见表3.1）是最简单的一种打分矩阵，其中，相同核苷酸匹配的得分为“1”，而不同核苷酸的替换得分为“0”（没有得分）。
 
（2）BLAST矩阵
 
BLAST是目前最流行的核酸序列比较程序，表3.2是其打分矩阵。这也是一个非常简单的矩阵，如果被比的两个核苷酸相同，则得分为“+5”，反之得分为“-4”。
 
（3）转换-颠换矩阵
 
核酸的碱基按照环结构分为两类，一类是嘌呤（腺嘌呤A，鸟嘌呤G），它们有两个环；另一类是嘧啶（胞嘧啶C，胸腺嘧啶T），它们的碱基只有一个环。如果DNA碱基的变化（碱基替换）保持环数不变，则称为转换（transition），如A->G，C->T；如果环数发生变化，则称为颠换（transversion），如A®C，A®T等。在进化过程中，转换发生的频率远比颠换高，而表3.3所示的矩阵正好反映了这种情况，其中转换的得分为“-1”，而颠换的得分为“-5”。
 
 
2、蛋白质打分矩阵
 
（1）等价矩阵
 
 
其中，Rij代表打分矩阵元素，i、j分别代表字母表第i个和第j个字符。
 
（2）遗传密码矩阵GCM
 
GCM矩阵通过计算一个氨基酸残基转变到另一个氨基酸残基所需的密码子变化数目而得到，矩阵元素的值对应于代价。如果变化一个碱基，就可以使一个氨基酸的密码子改变为另一个氨基酸的密码子，则这两个氨基酸的替换代价为1；如果需要2个碱基的改变，则替换代价为2；以此类推（见表3.4）。注意，Met到Tyr的转变是仅有的密码子三个位置都发生变化的转换。在表3.4中，Glx代表Gly、Gln或Glu，而Asx则代表Asn或Asp，X代表任意氨基酸。GCM常用于进化距离的计算，其优点是计算结果可以直接用于绘制进化树，但是它在蛋白质序列比对尤其是相似程度很低的序列比对中很少被使用。
 
 
（3）疏水矩阵
 
该矩阵（见表3.5）是根据氨基酸残基替换前后疏水性的变化而得到得分矩阵。若一次氨基酸替换疏水特性不发生太大的变化，则这种替换得分高，否则替换得分低。
 
 
（4）PAM矩阵
 
为了得到打分矩阵，更常用的方法是统计自然界中各种氨基酸残基的相互替换率。如果两种特定的氨基酸之间替换发生得比较频繁，那么这一对氨基酸在打分矩阵中的互换得分就比较高。PAM矩阵就是这样一种打分矩阵。PAM矩阵是第一个广泛使用的最优矩阵，它是基于进化原理的，建立在进化的点接受突变模型PAM（Point Accepted Mutation）基础上，通过统计相似序列比对中的各种氨基酸替换发生率而得到该矩阵。Dayhoff和她的同事们研究了71个相关蛋白质家族的1572个突变，发现蛋白质家族中氨基酸的替换并不是随机的，由此，断言一些氨基酸的替换比其他替换更容易发生，其主要原因是这些替换不会对蛋白质的结构和功能产生太大的影响。如果氨基酸的替换是随机的，那么，每一种可能的取代频率仅仅取决于不同氨基酸出现的背景频率。然而，在相关蛋白中，存在取代频率大大地倾向于那些不影响蛋白质功能的取代，换句话说，这些点突变已经被进化所接受。这意味着，在进化历程上，相关的蛋白质在某些位置上可以出现不同的氨基酸。
 
一个PAM就是一个进化的变异单位，即1%的氨基酸改变。但是，这并不意味着经过100次PAM后，每个氨基酸都发生变化，因为其中一些位置可能会经过多次改变，甚至可能变回到原先的氨基酸。因此，另外一些氨基酸可能不发生改变。PAM有一系列的替换矩阵，每个矩阵用于比较具有特定进化距离的两条序列。例如，PAM-120矩阵用于比较相距120个PAM单位的序列。一个PAM-N矩阵元素（i，j）的值反映两条相距N个PAM单位的序列中第i种氨基酸替换第j种氨基酸的概率。从理论上讲，PAM-0是一个单位矩阵，主对角线上的元素值为1，其它矩阵元素的值为0。其他PAM-N矩阵可以通过统计计算而得到。首先针对那些确信是相距一个PAM单位的序列进行统计分析，得到PAM-1矩阵。PAM-1矩阵对角线上的元素值接近于1，而其它矩阵元素值接近于0。例如，可以按下述方法构建PAM-1矩阵。首先，构建一个序列间相似度很高（通常大于85％）的比对。接着，计算每个氨基酸j的相对突变率mj。相对突变率就是某种氨基酸被其它任意氨基酸替换的次数。比如，丙氨酸的相对突变率是通过计算丙氨酸与非丙氨酸残基比对的次数来得到。然后，针对每个氨基酸对i和j，计算氨基酸j被氨基酸i替换的次数。最后，将以上替换次数除以对应的相对替换率，利用每个氨基酸出现的频度对其进行标准化，并将以上计算结果取常用对数，于是得到了PAM-1矩阵中的元素PAM-1（i，j）。这种矩阵被称作对数几率矩阵(log odds matrix)，因为其中的元素是根据每个氨基酸替换率的对数值来得到的。
 
 
 
将PAM-1自乘N次，可以得到矩阵PAM-N。虽然Dayhoff等人只发表了PAM-250，但潜在的突变数据可以外推至其他PAM值，产生一组矩阵。可以根据待比较序列的长度以及序列间的先验相似程度来选用特定的PAM矩阵，以发现最适合的序列比对。一般，在比较差异极大的序列时，通常在较高的PAM值处得到最佳结果，比如在PAM-200到PAM-250之间，而较低值的PAM矩阵一般用于高度相似的序列。实践中用得最多的且比较折衷的矩阵是PAM-250。
 
 
（5）BLOSUM矩阵
 
BLOSUM矩阵是由Henikoff首先提出的另一种氨基酸替换矩阵，它也是通过统计相似蛋白质序列的替换率而得到的。PAM矩阵是从蛋白质序列的全局比对结果推导出来的，而BLOSUM矩阵则是从蛋白质序列块（短序列）比对而推导出来的。但在评估氨基酸替换频率时，应用了不同的策略。基本数据来源于BLOCKS数据库，其中包括了局部多重比对（包含较远的相关序列，与在PAM中使用较近的相关序列相反）。虽然在这种情况下没有用进化模型，但它的优点在于可以通过直接观察而不是通过外推获得数据。同PAM模型一样，也有一系列的BLOSUM矩阵，可以根据亲缘关系的不同来选择不同的BLOSUM矩阵进行序列比较。然而，BLOSUM矩阵阶数的意义与PAM矩阵正好相反。低阶PAM矩阵适合用来比较亲缘较近的序列，而低阶BLOSUM矩阵更多是用来比较亲缘较远的序列。一般来说，BLOSUM-62矩阵适于用来比较大约具有62％相似度的序列，而BLOSUM-80矩阵更适合于相似度为80%左右的序列。
 

前言
 
前段时间一直在看时间序列相似性搜索(Time Series Similarity Search)的相关论文，现在终于放暑假了，开心度假中，也正好对那段时间读的论文做些总结。
 
首先来说明一下什么是时间序列(Time Series，以下简称时序)，时序就是按相等的时间采样的数据点构成的序列，数据点是几维的就叫几维时序。实际中一般以一维和二维时序居多。与时序类似的关键词还有轨迹(trajectory)，按我的理解他们的区别就是轨迹不一定是等时间间隔采样的，而时序一般是指按等时间间隔采样的序列。
 
其实时序数据挖掘是一个比较热门的研究领域，因为现实生活中很多东西都能看做时序，比如股票的波动数据，病人的心电图数据等，或者更复杂的采样时间间隔不等的轨迹，也能通过等间隔插值的方法近似成时序然后进行处理。 因此时序能描述的东西是很多的，对它的研究也进行了很长时间了，现在依然非常热门。
 
回到主题，那么时序搜索有什么用呢？如果能在历史的数据中找到与现在相似的时序，这可能能帮助我们预测时序未来的走势。此外，时序相似性搜索也是很多其他时序数据挖掘的基础，比如时序数据的分类和聚类等。
 
正文
 
下面开始说明怎样进行时序相似性搜索，我们的目的是给定一个时序查询Q，然后从一个时序数据库中返回与Q最相似的时序。
 
数据预处理
 
首先，由于时序是典型的高维数据，数据点可能相当多（你想啊，心电图可能几毫秒就采样一个点，这存下来几年的数据得有多少，当然，一条时序具体有多少个点还得看你怎么分割这些数据），我们需要预处理，需要将原始数据以一种合适的方式来呈现，以方便后期的处理，这种技术叫做时序呈现(time series representation),其实也就是维度缩减技术，目的是压缩数据但是保留主要信息。这类技术包括离散傅里叶变换(DFT)，离散小波变换(DWT)，主成分分析(PCA),奇异值分解(SVD),PAA,SAX等，这些东西可以在一些综述中找到，如09年的一篇综述“高效时序相似搜索技术“。
 
另外，为了减少平移和缩放对相似性的影响，我们需要对原始数据进行标准化，常用的方法是z标准化(Z-score)，就是每个数据点减去平均值，再除以偏标准差。公式为
 $$\frac{{x - \mu }}{\sigma }$$  

 其中

$\mu$表示平均值，

$\sigma$表示标准差。 

 对每一维数据经过这样处理后，每一维数据的平均值就变为0，标准差为1. 

 至此，数据预处理部分完毕。
 
相似性测量
 
为了比较两个时序，我们需要一种评价方法来测量他们的相似性。最简单的，可以用欧氏距离(Euclidean Distance，以下简称ED)来测量。公式为
 $$ED(S,Q)=\sum\limits_{i = 1}^n {{{({S_i} - {Q_i})}^2}}$$  

 但是欧式距离的缺点也很明显，一是它只能测量长度相等的时序，二是它对噪声是敏感的，因为个别很远的偏离平均值的点可能对结果造成很大的影响。 

 为了克服这些缺点，人们又发明出了其他更好的测度。比如动态时间规整(dynamic time warping，DTW)，最长公共子序列(LCSS)等。 

 他们的公式为
 

 公式中dist为点与点之间的距离计算函数，一般选用欧氏距离。 

 推荐一篇论文”Experimental comparison of representation methods and distance measures for time series data”，这篇论文通过实验比较分析了不同的数据呈现技术以及时序测量技术。有很多论文认为DTW是最好的测量，原因一是DTW不需要设定参数，二是DTW虽然是平方算法，但是它有很多能线性计算的下界，利用下界可以进行剪枝，据”Addressing Big Data Time Series: Mining Trillions of Time Series Subsequences under Dynamic Time Warping”称，利用多种下界进行剪枝，在大规模数据集中剪枝率能达到99.9999%，也就是说实际计算DTW效率基本可以达到线性。三是利用DTW可以快速的计算子序列相似性，所谓子序列相似性就是给定序列S,Q,找到S的子序列使得它与Q的相似性最大。利用SPRING算法计算子序列相似性与计算全序列相似性的代价是一样的，具体可以参见论文”Stream Monitoring under the Time Warping Distance”
 
论文”Addressing Big Data Time Series: Mining Trillions of Time Series Subsequences under Dynamic Time Warping”给出了可以用于工业进行时序相似性搜索的方法，基本思想是结合DTW的多种下界进行剪枝。关于DTW的下界可以参考论文”Boundary-based lower-bound functions for dynamic time warping 
 and their indexing”

这是一篇阅读笔记。
 原文 An Introduction to Sequence Similarity (“Homology”) Searching by William R. Pearson（ 原文地址），作者是FASTA格式的发明者之一。
 
同源 Homology
 
定义
 
In biology, homology is similarity due to shared ancestry between a pair of structures or genes in different taxa (wikipedia).
 
同源 （Homology） = 共同进化祖先（Common evolutionary ancestry）
 只要有共同祖先，无论基因序列、其编码蛋白质的结构、功能是什么，都可看作是同源的。
 
为什么要寻找同源基因
 
一旦发现同源序列，就可以通过多序列比对建立更准确的比对，为后续的表型预测和进化分析奠定基础。
 
识别同源序列的策略
 
相似性搜索 （similarity searching）
 序列相似性搜索可以通过检测过高的相似性来识别同源蛋白质或基因：当两个序列的相似性超过偶然的预期时，我们推断这两个序列存在同源性。 当观察到过高的相似性时，最简单的解释是，这两个序列不是独立出现的，它们起源于一个共同的祖先。
 所以这是统计学意义上的同源性，显著的相似性一定程度上反映了同源性。
 
需要注意的是同源性与相似性是两个不同的概念！
 两条高度相似的序列可能不存在同源性；同样的，同源序列的相似性也可能很低。例如两条同源序列的相似性比对结果不显著，但如果它们在结构上相似性上显著，或者它们都与第三条序列的相似性显著，那么它们显然是同源序列。 因此，当相似性搜索发现统计学上显着的匹配时，我们可以放心地推断出这两个序列是同源的。 但是，如果在数据库中找不到统计上显着的匹配项，则不能确定没有同源物。
 
常见的序列比对工具，如BLAST， FASTA，HMMER等在算法上尽量减少假阳性（false positives, non-homologs with significant scores; Type I errors）的发生，但对假阴性（false negatives, homologs with non-significant scores; Type II errors）没有约束。
 
如果在InterPro和Pfam等域注释库中没有找到注释的蛋白质域，那是因为查询序列与已知的同源序列的同源关系太远。
 
期望（$E$）的计算公式：
 $$E=kmn{e}^{-\lambda S}$$
 $E$：期望值，即分数为S时，期望的高分序列（HSP）出现的数量 ；
 $\lambda ,k$：常数（Karlin Altschul statistics）；
 $m$：查询序列长度；
 $n$：数据库序列的长度。
 
期望值取决于数据库的大小，通过对比拥有10,000,000个序列数据库得到的e值比只有100,000个序列的数据库中找到相同分数时的e值低100倍。但并不是说在大的数据库中找到的序列是同源的，而小的数据库中找到的序列不是同源的。
 
与蛋白质/蛋白质比对相比，DNA/DNA序列比对比可能更不容易发现同源性。蛋白质（或者翻译后的DNA）相似性搜索要比DNA/DNA搜索敏感得多。 经过200-400亿年的演化后，DNA：DNA比对比对很少能检测到同源性，而对于蛋白质/蛋白质比对能检测到25亿年前的共同祖先。
 此外，DNA/DNA比对不如蛋白质/蛋白质准确。E值<0.001的蛋白质/蛋白质比对可以可靠地推断同源性，DNA/DNA期望值<10e-6经常是偶然发生的，一般阈值设为10e-10。提高DNA序列搜索灵敏度的最有效方法是使用翻译的DNA/蛋白质比对，例如BLASTx和FASTX产生的比对，​​而不是DNA/DNA比对。