Overview

最近看了几种生物信息学领域的特征提取算法,在观察了已有算法的特点之后,也想尝试着设计新的算法。这里以CKSAAP为基础,介绍一下我们想设计的算法的动机和原理。

CKSAAP的核心原理

关于CKSAAP的算法介绍详细部分可以参考Young写的 蛋白质序列特征提取方法之——CKSAAP
CKSAAP试图统计刻画出每条序列中不同的氨基酸对出现的特征,它的核心原理可以用下面这个公式表述:

(NLENTotal,NEENTotal,NEYNTotal,,NDLNTotal,NLLNTotal)400

就是统计氨基酸对的出现频率。使用一个氨基酸(例如LE)对在一个序列中出现的次数,除以这个序列中包含的氨基酸对总数。

Motivation

这种方法不能准确地表示氨基酸对在一个序列中出现的关联度,如果它们真的有关联的话。
CKSAAP中,一个氨基酸对(例如LE)的在特征矩阵中取值高,并不一定是因为关联性高,也有可能是因为LE单独出现的概率很大,这样的话从概率上来说它们在一起出现的次数也会更多,但并不一定是存在某种关联;另一方面,另一个氨基酸对(例如TW),可能TW在这个序列中都只出现了两次,但是每次出现,它们都是结对出现,那这种关联性是不是更强呢。
虽然说这两种都是比较极端的情况,但是将氨基酸对的出现次数和这两个氨基酸单独出现的次数一起考虑,可能可以更好地表征氨基酸对的关联性。

设计

我们采用点互信息熵的概念来描述这种关联性,使用下面的公式:

I(x,y)=logP(x,y)P(x)P(y)

其中,

P(x)=N(x)N(total),P(y)=N(y)N(total),P(x,y)=N(x,y)N(totalpair)

P(x)P(y)P(x,y) 分别为一条序列中x出现的概率,y出现的概率,以及xy序列对出现的概率。

细节讨论

我们通过一个例子讨论这个模型可能出现的问题:
假设有一个长度为100的序列,则N(total)=100N(totalpair)=99,为了方便计算和观察,我们令N(totalpair)100
在以下两种情况下,我们使用这个模型求氨基酸AC的点互信息熵。

  • Case 1: AC各只出现1次,AC出现1次

    P(A,C)=N(A,C)/Ntotalpair=1/100
    P(A)=N(A)/Ntotal=1/100
    P(C)=N(C)/Ntotal=1/100

    P(A,C)/(P(A)P(C))=1100/(11001100)=100

    I(A,C)=logP(A,C)P(A)P(C)=log100

  • Case 2: AC各出现50次,AC出现50次

    P(A,C)=N(A,C)/Ntotalpair=50/100
    P(A)=N(A)/Ntotal=50/100
    P(C)=N(C)/Ntotal=50/100

    P(A,C)/(P(A)P(C))=50100/(5010050100)=2

    I(A,C)=logP(A,C)P(A)P(C)=log2

  • Case 3: AC各只出现10次,AC出现1次

    P(A,C)=N(A,C)/Ntotalpair=1/100
    P(A)=N(A)/Ntotal=10/100
    P(C)=N(C)/Ntotal=10/100

    P(A,C)/(P(A)P(C))=1100/(1010010100)=1

    I(A,C)=logP(A,C)P(A)P(C)=log1