您正在查看: 2015年10月

使用Apache反向代理转发Tomcat请求

Overview

最近在为我们的一篇论文做Online Server,我们使用了Java的Struts框架搭建了Server,并且部署到了Tomcat上,这也带来了一个问题。服务器上本来已经有了Apache,并且占用了80端口,加上Tomcat之后使用域名必须加上Tomcat的端口号访问,使用起来很不方便,Online Server的网址看起来也很奇怪。
整合ApacheTomcat有好多种方法,这里我们使用Apache反向代理的方式。其他方法概览请查看这里

1. 服务器信息

下面以我们的云服务器为例说明这个过程,信息如下:

服务器域名 : www.nohup.cc
Apache端口号: 80
Tomcat端口号: 8888
Tomcat上的应用程序: hudson
Apache程序路径:/etc/apache2

当我们使用域名访问服务器时,服务器会默认把请求转发到80端口的服务器,在我们这里就是Apache服务器。因此访问Apache上的网站时,可以直接使用域名,比如可以通过 www.nohup.cc 直接访问本站点;而访问Tomcat上的网站时,需要添加端口号访问,比如访问hudson,需要使用 www.nohup.cc:8888/hudson才可以访问。

2. 配置Apache

我们通过配置Apache来使用反向代理,步骤如下:

  1. 进入Apache文件目录:

    cd /etc/apache2
    
  2. 修改Apache配置文件apache2.conf,在文件末尾添加一下几行:

    # 导入需要用到的两个模块
    LoadModule proxy_module modules/mod_proxy.so
    LoadModule proxy_http_module modules/mod_proxy_http.so
    
    # 下面的两条用来配置访问hudson的反向代理
    ProxyPass /hudson http://www.nohup.cc:8888/hudson/
    ProxyPassReverse /hudson http://www.nohup.cc:8888/hudson/
    
  3. 重启Apache

    sudo apache2ctl -k restart
    

    或者使用下面的命令重启

    sudo service apache2 restart
    

    前者报错信息比较直接,直接打印到了屏幕上。

    2.1 可能出现的问题

    • 如果你得到下面的信息:

        [Mon Oct 26 09:26:39.080381 2015] [so:warn] [pid 21343] AH01574: module proxy_module is already loaded, skipping
        [Mon Oct 26 09:26:39.080615 2015] [so:warn] [pid 21343] AH01574: module proxy_http_module is already loaded, skipping
      

      说明你的Apache服务器已经导入了这两个模块,不需要重复导入,删掉

        LoadModule proxy_module modules/mod_proxy.so
        LoadModule proxy_http_module modules/mod_proxy_http.so
      

      这两句就可以了。

    • 你也可能得到的是下面的报错信息:

        apache2: Syntax error on line 222 of /etc/apache2/apache2.conf: Cannot load modules/mod_proxy.so into server: /etc/apache2/modules/mod_proxy.so: cannot open shared object file: No such file or directory
        Action '-k restart' failed.
        The Apache error log may have more information. 
      

      错误信息很明确,查看/etc/apache2/mods-enabled文件夹,发现里面并没有proxy.*proxy_http.*文件,但是/etc/apache2/mods-available文件夹里面是有的,说明apache已经安装了这些模块,但是并未开启。
      使用下面的命令开启mod_proxymod_proxy_http模块:

        sudo a2enmod proxy
        sudo a2enmod proxy_http
      

      现在再看/etc/apache2/mods-enabled文件夹,里面已经有了proxy.*proxy_http.*文件,重启apache让配置生效,就不会再报错了。

现在就可以使用www.nohup.cc/hudson来访问Tomcat上的hudson程序了...

3. 配置Apache过滤特定的请求

考虑下面的需求:
我们想要把www.nohup.cc的请求全部转发到tomcat上处理,我们可以修改Apache配置文件apache2.conf

# 下面的两条用来配置所有请求的反向代理
ProxyPass / http://www.nohup.cc:8888
ProxyPassReverse / http://www.nohup.cc:8888

这样,当我们访问www.nohup.cc时,请求就会被自动转发到tomcat
假定我们在Apache上还有一个网站apache_website,可当我们访问www.nohup.cc/apache_website时,会被定位到http://www.nohup.cc:8888/apache_websitetomcat上没有这个网站,自然也就访问不到。我们通过修改Apache配置文件apache2.conf,告诉Apache,如果是/apache_website,就不要转发给tomcat了:

# 下面这条配置使得apache在转发请求时,过滤掉apache_website这个站点
ProxyPass /apache_website !

重启就可以生效。

实际上,在使用中,我们通常会使用Apache的过滤功能过滤掉静态资源(如html,图片,CSS等),这些请求由Apache处理的话速度要快很多,而动态的页面才转发给tomcat处理。
关于这些问题的配置请查看可以文末的链接。

参考网址

蛋白质序列特征提取方法之——PSSM

Overview

我在之前写的一篇博客中谈到整理那些混乱的数据源,发现有pssm fts文件夹中的子文件夹和文件并不清楚来龙去脉,这个问题困扰了我一段时间。最近在研究PSSM算法时,与Chris交流了一下,恍然大悟:这个文件夹中的t3pssm,t4pssm,t6pssm三个子文件夹中的形如t6_12.pssm的文件族,是由t3,t4,t6这三个文件夹中的形如t6_12.fasta的文件族经过BLAST或者PSI-BLAST等算法程序处理过后得到的数据。并且t3pssm,t4pssm,t6pssm三个子文件夹中的每个子文件中包含有两个矩阵,其中第一个即为我们在构建PSSM的步骤中提到的PSSM。这些PSSM将会在下一步特征提取的过程中做为比对(alignment)的模板。
由于并未拿到BLASTPSI-BLAST程序的源码,所以这里只介绍从PSSM开始到提取特征的步骤。

1.Python编程方案

我们的项目在提取特征这一阶段的第六个步骤中,运用的就是这个PSSM方法,这一步只有一个Linux命令

python t34pssm.py T4undrsmp.txt ./t4 ./t4pssm

其中,t34pssm这个Python文件调用了另一个Python文件pssm_ft.py,这两个文件是用来计算提取特征的;T4undrsmp.txt是输入文件,即被处理数据源;./t4./t4pssm是两个文件夹。
下面我们详细解析这个Python编程方案的每一步的代码和作用。在解析每一步的输出时,我在原代码的基础上引入了一个time

import time

这个包提供给我们一个很好的函数sleep(),这个函数输入的数基本时间单位为秒。这样我们在遇到大量循环的时候,可以在循环体内加上这么一句:

time.sleep(1)
print #所需要的打印内容

就可以观察每一步循环的输出内容了,睡眠时间可以根据需要调整,十分方便。在终止循环时,可以在命令行界面按下ctrl+c停止进程,这样就能记录下输出内容了。

1.1 解析t34pssm.py文件

引入相关的包:

import fileinput
import sys
from os import listdir
from os.path import isfile, join
import re
from pssm_ft import * #引入pssm_ft.py文件
import time

读取输入参数[1][2]

smplfasta = sys.argv[1]           #T4undrsmp.txt
spfasta = sys.argv[2]             #./t4
check_head = re.compile(r'\>')    #正则表达式匹配符号>

有一点需要说明:由于T4undrsmp.txt这个文件中原有的数据量太大,故输出量太大,耗时太大,我们仅取其第一条蛋白质序列作为例子说明这个算法的运算过程。第一条蛋白质序列如下:

>lpn:lpg0012 hypothetical protein (A)|1
MNTTEHTELGNGLLMDKIEGNPYLRIDEFGVLHLRLQRFGENNLPEPMELEMSAGEIIAMAGDYFTQANW
TMDLDLPKCELFNSPAELGKHLIRKPIDPKEENALITAYNNLAAPDVTRKEIDRIYSINNANYVPFSPTL
NFYAQQLMYYFRVKDYGEMLVRNQTHFTPWSIRVYILGHAIALRYARLSYELKQLATDRNYQSDNPDLST
LKISLQNKNETLSSNTLLDLANRYHAQAYSIELFTFHYYSDHFATGHMSMIGDLRVVLKERFGVWGNILA
NNLHDEVNRVGVYTVRPYDPTPNTTEAPSRARGDGKFDTCLNQFNRLACLNGMTASLKDINQVLSGGVIP
EQKKFGGLVHMPDVDFNSRQHQPLLVLSKGKVYYRNNLSQIHIISPSEYEALRENPEKHGYKELTSKWAA
FKLVTKLRLFPYLYDGSVMPVSDEKLAEIIADEKRRNPQRAPIPTPSCLPESEPTVFDWRTKASWRNNKD
SLDILDGLKKHSILAAKQTHSAIQEEEEVRLNLGV

T4undrsmp.txt文件中读取undersample fasta格式的蛋白质序列:

smplist = []              #存储蛋白质序列的名称序号信息
smpcnt = 0                #记录读取的蛋白质数量
for line, strin in enumerate(fileinput.input(smplfasta)):
if check_head.match(strin):
    smplist.append(strin.strip())
    #print smplist        #['>lpn:lpg0012 hypothetical protein (A)|1']
    #time.sleep(1)        #休眠一秒以观察输出
    smpcnt += 1

读取1232fasta格式的文件:

onlyfiles = [ f for f in listdir(spfasta) if isfile(join(spfasta,f)) ]  
#print onlyfiles    
#输出t4_1.fasta等1232个文件的列表,['t4_961.fasta', 't4_456.fasta', 't4_1179.fasta',……, 't4_954.fasta', 't4_567.fasta']
onlyfiles.remove('.DS_Store')        #移除干扰项

将蛋白质名和文件名对应起来:

fastaDict = {}                          #存储“蛋白质名:文件名”的键值对

for fi in onlyfiles:
    cntnt = ''
    for line, strin in enumerate(fileinput.input(spfasta+'/'+fi)):  
        if line == 0:
            cntnt += strin.strip()      #去除前后空格
    if cntnt in fastaDict:              #检查否重复(已经存在)
        print strin                     #打印重复(已经存在)
    fastaDict[cntnt] = fi               #生成键值对
print fastaDict     #{'>lpn:lpg0012 hypothetical protein (A)|1': 't4_1.fasta',……,'>gi|148360344|ref|YP_001251551.1| hypothetical protein LPC_2282 [Legionella pneumophila str. Corby]|-1': 't4_580.fasta'} 一共1232个键值对  

读取输入参数[3]

pssmdir = sys.argv[3]             #PSSM模板文件目录
finalist = []                     #存储PSSM文件名字的列表 
for smp in smplist:
    #print fastaDict[smp]        #t4_1.fasta
    finalist.append(pssmdir+'/'+fastaDict[smp].split('.')[0]+'.pssm')       #根据fasta文件的名字找出对应的PSSM文件,输出t4_1.pssm
    #print finalist               #['./t4pssm/t4_1.pssm']

for fi in finalist:               #读取PSSM文件列表
    pssm(fi, 't4pssm.txt')        #调用函数文件pssm_ft.py中的pssm(fi,output)函数

调用函数pssm(fi,output)这一步表明真正计算的步骤将要开始,且是t4pssm.txtPSSM之间的运算,与其他文件无关,更印证了我们最初的想法。

1.2 解析pssm_ft.py文件

这个pssm_ft.py文件的作用就是Retrieve PSSM features,也就是检索PSSM特征。将T4undrsmp.txt这个文件中的序列跟PSSM进行比对,得出蛋白质序列的特征向量。
引入相关的包:

import sys
import numpy as np
import math
import re
import fileinput 
import time

numpy这个包是用来处理大型矩阵运算的。
既然pssm(fi,output)这个函数首先被调用,我们就先分析一下它。

def pssm(fi,output):    #pssm(fi, 't4pssm.txt')     fi==./t4pssm/t4_1.pssm
    Amino_vec = "ARNDCQEGHILKMFPSTWYV"    #20个基本氨基酸
    
    PSSM = []            #存储PSSM的特征向量矩阵
    ACC = [ [0.0] * 20 ] * 20
    seq_cn = 0
    for line, strin in enumerate(fileinput.input(fi)):
        if line > 2:
        str_vec = strin.split()[1:22]    #切片截取该行的第1个至第21个字符存入向量str_vec中
        #print str_vec   #['M', '-2', '-3', '-4', '-5', '-3', '-2', '-4', '-4', '-3', '0', '3', '-3', '8', '-1', '-4', '-3', '-2', '-3', '-2', '0']等525个向量
        if len(str_vec) == 0:
            break
        PSSM.append(map(int, str_vec[1:]))    #将每个向量中的氨基酸符号去掉,其他字符转换为int类型,并全部存入PSSM特征向量矩阵    
        seq_cn += 1        #525
        ACC[Amino_vec.index(str_vec[0])] = map(sum, zip(map(int, str_vec[1:]), ACC[Amino_vec.index(str_vec[0])]))    #这句最长,从赋值号(=)右边开始逐步分析

"""
1. Amino_vec.index(str_vec[0]) 返回本行氨基酸符号在Amino_vec中的位置
2. ACC[Amino_vec.index(str_vec[0])] 初始值均为20维向量[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
3. map(int, str_vec[1:]) 将每个向量中的氨基酸符号去掉,其他字符转换为int类型
4. zip(map(int, str_vec[1:]), ACC[Amino_vec.index(str_vec[0])]) 返回一个tuple列表,形如[(-2,0.0),(-3,0.0),……,(-2,0.0),(0,0.0)]
5. map(sum, zip(map(int, str_vec[1:]), ACC[Amino_vec.index(str_vec[0])])) 将tuple列表中各元素转换为sum类型,即求和
6. 循环迭代,最后ACC(20维)为:
[[85.0, -27.0, -42.0, -58.0, -59.0, -14.0, -15.0, -25.0, -46.0, -53.0, -55.0, -20.0, -32.0, -75.0, -32.0, 10.0, -17.0, -88.0, -57.0, -28.0], [-23.0, 90.0, -3.0, -15.0, -80.0, 24.0, -3.0, -52.0, -10.0, -45.0, -42.0, 19.0, -31.0, -63.0, -44.0, -17.0, -10.0, -68.0, -32.0, -48.0], …… , [-5.0, -33.0, -36.0, -34.0, -36.0, -30.0, -26.0, -58.0, -49.0, 23.0, 11.0, -19.0, -12.0, -36.0, -36.0, -31.0, -5.0, -72.0, -35.0, 42.0]]

"""
    fileinput.close()
    
    ACC_np = np.array(ACC)    #转换成数组
    ACC_np = np.divide(ACC_np, seq_cn)    #将每一个元素除以525
    amcnt_max = np.amax(ACC_np)    #找出最大的元素0.28380952381
    amcnt_min = np.amin(ACC_np)    ##找出最小的元素-0.274285714286
    
    vfunc = np.vectorize(scale) #矢量化scale函数
    """
    numpy中的vectorize函数可以将scale函数转换成可以接受向量参数的函数
    """
    ACC_np = vfunc(ACC_np, amcnt_max,amcnt_min)    #调用scale函数,将矩阵范围控制在[-1,1]之间
    ACC_shp = np.shape(ACC_np)    #(20,20)
    PSSM = np.array(PSSM)
    PSSM_AC = pssm_ac_cal(PSSM, GROUP, seq_cn)      #调用函数pssm_ac_cal(PSSM, g, l)   
"""
pssm_ac_cal(PSSM, g, l)函数,三个参数:1.PSSM矩阵;2.10;3.蛋白质序列长度
返回10×20的矩阵PSSM_AC
"""               
    PSSM_shp = np.shape(PSSM_AC)    #(10, 20)
    file_out = file(output,'a')    #打开文件t4pssm.txt,如果没有则创建一个;a表示add,以追加的方式打开
    np.savetxt(file_out, [np.concatenate((np.reshape(ACC_np, (ACC_shp[0] * ACC_shp[1], )), np.reshape(PSSM_AC, (PSSM_shp[0] * PSSM_shp[1], ))))], delimiter=",")    #这句虽然比较长,但是逻辑比较简单清晰
"""
1. reshape函数将矩阵的行、列、维重新调整为200×1的矩阵.
2. concatenate函数将两个矩阵连接起来成为一个新矩阵.
3. savetxt函数将矩阵保存为txt格式文件,分隔符为",".
"""     

现在看一下另外两个函数scalepssm_ac_cal
scale函数:

def scale(x, max, min):
"""
把矩阵中最大的数和最小的数变为1和-1,其他数按以下比例缩放
"""
    try:
        return -(1-(x-min)/(max-min))+((x-min)/(max-min))
    except Exception:
        print "ZeroDivisionError, max-min:\s", max-min

pssm_ac_cal函数:

def pssm_ac_cal(PSSM, g, l):
    """
    输入PSSM、GROUP、蛋白质序列的长度,返回GROUP×20的PSSM_AC矩阵
    """
    PSSM_AC = np.array([ [0.0] * 20 ] * g)

    for pg in range(g):
        l_g = l - pg - 1
        for pj in range(20):
            sum_jl = 0.0
            for i in range(l):
                sum_jl += PSSM[i][pj]
            sum_jl /= l
        
            pssm_acjg = 0.0
            for i in range(l_g):
                pssm_acjg += (PSSM[i][pj]-sum_jl) * (PSSM[i+pg+1][pj]-sum_jl)
        pssm_acjg /= l_g
            PSSM_AC[pg][pj] = pssm_acjg
    return PSSM_AC

至此,Python编程方案已解析完毕。疏漏之处请与我们联系,一定尽快修改。

构建PSSM的步骤

Overview

PSSM算法是生物信息学领域中的一个常用算法,全名“位置特异性打分矩阵(position-specific scoring matrix)”,又称作"位置比重矩阵(position weight matrix)".有关该方法更多的细节,详见维基百科Position weight matrix.本文仅阐述其设计思想,实际项目的例子将在另一篇文章中进行介绍.

1.建模思想

一个PWM包含N行(列),当模型为DNA时,N=4;当模型为蛋白质时,N=20。维基百科给出的是DNA的核苷酸序列的例子,对于我们的项目来说,不同的也仅仅是矩阵行列数而已:组成DNA的基本核苷酸有A,C,G,T四种,故行列式有四行(列);组成蛋白质的基本氨基酸有二十种,故行列式有二十行(列)。同时PWM对于每个不同位置都对应一列(行)数据。而选择基本核苷酸(氨基酸)作为行或列,需要针对各自的需要进行调整。

1.1 构建位置频度矩阵(PFM)

位置频度矩阵(position frequency matrix)的构建的思路还是比较简单的,这里仍然以维基百科的例子介绍。DNA序列如下表所示:

GAGGTAAAC
TCCGTAAGT
CAGGTTGGA
ACAGTCAGT
TAGGTCATT
TAGGTACTG
ATGGTAACT
CAGGTATAC
TGTGTGAGT
AAGGTAAGT

共有109列,则相应的PFM

$$ M=\begin{matrix} A\\ C\\ G\\ T\end{matrix}\begin{bmatrix} 3 &6 &1 &0 &0 &6 &7 &2 &1 \\ 2 &2 &1 &0 &0 &2 &1 &1 &2 \\ 1 &1 &7 &10 &0 &1 &1 &5 &1 \\ 4 &1 &1 &0 &10 &1 &1 &2 &6 \end{bmatrix} $$

这个矩阵的列数也是9,而每一列的数各自加起来正好是10,也就是DNA序列的行数。这样,构建矩阵的原理就很清晰了:计算每一列中的各核苷酸的数量,然后存入矩阵的相应位置

1.2 构建位置概率矩阵(PPM)

通过PFMPPM,只需要下面的公式:

$$ M_{k,j}=\frac{1}{N}\sum_{i=1}^N I(X_{i,j}=k) $$

其中,i为行号,j为列号,即:

$$ i\in \begin{pmatrix} 1, &2, &\cdots, &N \end{pmatrix},j\in \begin{pmatrix} 1, &2, &\cdots, &l \end{pmatrix}. $$

I是指示函数,即:

$$ \mathbf{\mathit{I}}_(X_{i,j}=k) = \begin{cases} 1 &,\text{if } X_{i,j} = k ; \\ 0 &,\text{if } X_{i,j} \neq k . \end{cases} $$

则1.1中的PFM相应的PPM为:

$$ M=\begin{matrix} A\\ C\\ G\\ T\end{matrix}\begin{bmatrix} 0.3 &0.6 &0.1 &0.0 &0.0 &0.6 &0.7 &0.2 &0.1 \\ 0.2 &0.2 &0.1 &0.0 &0.0 &0.2 &0.1 &0.1 &0.2 \\ 0.1 &0.1 &0.7 &1.0 &0.0 &0.1 &0.1 &0.5 &0.1 \\ 0.4 &0.1 &0.1 &0.0 &1.0 &0.1 &0.1 &0.2 &0.6 \end{bmatrix} $$

在这里有一点需要特别注意

对于每一个匹配成功的核苷酸,我们计分为1,未匹配则记为0,这只是一个简单的思想。而在实际情况下,我们用BLAST或者PSI-BLAST等程序求序列的PSSM时,相应的打分规则就要复杂许多。这时候就需要根据自身的需要选择相应的各种参数,如:gap,λ。

关于上面一段话中提及的程序,更多的细节请点击:
https://en.wikipedia.org/wiki/BLAST
http://blast.ncbi.nlm.nih.gov/Blast.cgi
http://www.ebi.ac.uk/Tools/sss/psiblast/
http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE=Proteins&PROGRAM=blastp&RUN_PSIBLAST=on

我们假定这些矩阵在行列式的每个位置上都是统计独立的。非统计独立模型请参考Chris的一篇文章http://www.nohup.cc/article/111/,这里就不赘述了。那么对于DNA序列S=GAGGTAAAC来说,它出现的概率就是

$$ p\left ( S\mid M \right )= 0.1\times 0.6\times0.7\times1.0\times1.0\times0.6\times0.7\times0.2\times0.2=0.0007056 $$

1.3 构建位置比重矩阵(PWM)

这里引入一个参数bb=1/k,其中,当序列为蛋白质序列时,k=20;序列为DNA时,k=4.那么对于相同的位置上的PWMPPM的矩阵元素,其关系为:

$$ M_{PWM}=ln\left ( \frac{M_{PPM}}{b} \right ) $$

那么1.2中的PPM相对应的PWM为:

$$ M=\begin{matrix} A\\ C\\ G\\ T\end{matrix}\begin{bmatrix} 0.18 &0.87 &-0.91 &-\infty &-\infty &0.87 &1.02 &-0.22 &-0.91 \\ -0.22 &-0.22 &-0.91 &-\infty &-\infty &-0.22 &-0.91 &-0.91 &-0.22 \\ -0.91 &-0.91 &1.02 &1.38 &-\infty &-0.91 &-0.91 &0.69 &-0.91 \\ 0.47 &-0.91 &-0.91 &-\infty &1.38 &-0.91 &-0.91 &-0.22 &0.87 \end{bmatrix} $$

至此,我们已经介绍完构建PSSM的基本思想。

一种基于点互信息熵的特征提取算法

Overview

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

CKSAAP的核心原理

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

$$ (\frac{N_{LE}}{N_{Total}},\frac{N_{EE}}{N_{Total}},\frac{N_{EY}}{N_{Total}},……,\frac{N_{DL}}{N_{Total}},\frac{N_{LL}}{N_{Total}})_{400} $$

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

Motivation

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

设计

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

$$ I_{(x,y)} = log \frac{P_{(x,y)}}{ P_{(x)}P_{(y)} } $$

其中,

$$ P_{(x)} = \frac{ N_{(x)} }{ N_{(total)} }, P_{(y)} = \frac{ N_{(y)} }{ N_{(total)} }, P_{(x,y)} = \frac{ N_{(x,y)} }{ N_{(totalpair)} } $$

即 $P_{(x)}$,$P_{(y)}$,$P_{(x,y)}$ 分别为一条序列中x出现的概率,y出现的概率,以及xy序列对出现的概率。

细节讨论

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

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

    $P_{(A,C)} = N_{(A,C)} / N_{totalpair} =1/100$
    $P_{(A)} = N_{(A)} / N_{total} = 1/100$
    $P_{(C)} = N_{(C)} / N_{total} = 1/100$
    $\Longrightarrow$
    $P_{(A,C)} / (P_{(A)}P_{(C)}) = \frac{1}{100} / (\frac{1}{100} \ast \frac{1}{100}) = 100$

    $I_{(A,C)} = log \frac{P_{(A,C)}}{ P_{(A)}P_{(C)}} = log 100$

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

    $P_{(A,C)} = N_{(A,C)} / N_{totalpair} =50/100$
    $P_{(A)} = N_{(A)} / N_{total} = 50/100$
    $P_{(C)} = N_{(C)} / N_{total} = 50/100$
    $\Longrightarrow$
    $P_{(A,C)} / (P_{(A)}P_{(C)}) = \frac{50}{100} / (\frac{50}{100} \ast \frac{50}{100}) = 2$

    $I_{(A,C)} = log \frac{P_{(A,C)}}{ P_{(A)}P_{(C)}} = log 2$

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

    $P_{(A,C)} = N_{(A,C)} / N_{totalpair} =1/100$
    $P_{(A)} = N_{(A)} / N_{total} = 10/100$
    $P_{(C)} = N_{(C)} / N_{total} = 10/100$
    $\Longrightarrow$
    $P_{(A,C)} / (P_{(A)}P_{(C)}) = \frac{1}{100} / (\frac{10}{100} \ast \frac{10}{100}) = 1$

    $I_{(A,C)} = log \frac{P_{(A,C)}}{ P_{(A)}P_{(C)}} = log 1$

蛋白质序列特征提取方法之——CKSAAP

Overview

CKSAAP(Compositon of k-spaced Amino Acid Pairs)方法中,利用在蛋白质序列片断中k个间隔距离的残基对(residue pairs)在该序列中的组成比例,建立数学模型,提取出特征向量,从而达到预测泛素(Ubiquitin)的目的。
残基(residue)和泛素(Ubiquitin)信息详见维基百科:残基泛素,这里就不赘述了。
本文主要介绍该方法建立数学计算模型的思路及编程方案。

1.建模思路

以长度为48的序列LEEYRKHVAERAAEGIAPKPLDANQMAALVELLKNPPAGEEEFLLDLL为例(如果你不懂这些英文字母什么含义,请点击氨基酸简写符号),当k=0时,我们需要提取的残基对(residue pairs)为{LE,EE,EY,……,LD,DL,LL},即每个氨基酸和它相邻的下一个氨基酸组成一对提取出来,也就是说这两个氨基酸中间的间隔距离是k=0个氨基酸。以此类推,当k=1时,我们需要提取的残基对为{LE,EY,ER,YK,……,LD,LL,DL}……本方法中,k最大为5.
细心的你一定会发现k=0时候,结尾会有1个氨基酸没有配对,而提取的残基对数量为47k=1时,有2个氨基酸没有配对,而提取的残基对数量为46;所以,规律就是,当序列长度为N,间隔为k时,一共可以提取的残基对数量为N-k-1,记为

$$ N_{Total}=N-k-1 $$

这里要特别感谢一下Chris同学的提醒和点拨。
由于基本氨基酸数量为20,故而可以形成的残基对数量是20×20=400.我们统计的是这些残基对在这个蛋白质序列当中出现的概率,于是便产生了一个400维的特征向量,即

$$ (\frac{N_{LE}}{N_{Total}},\frac{N_{EE}}{N_{Total}},\frac{N_{EY}}{N_{Total}},……,\frac{N_{DL}}{N_{Total}},\frac{N_{LL}}{N_{Total}})_{400} $$

需要注意的是,本方法中设定的N为窗口值为window=27,即截取的片段长度为27(具体的截取方法将在下面给出)。可以想象,组成这个向量的400个比值当中,绝大部分为0.0,而这些比值总和为1.所以,当k分别取0,1,2,3,4,5等值时,一共会有400×6=2400个向量。

2.编程方案

encode这个文件夹中的算法,涵盖了数种主流的算法,这些算法介绍均会相继给出。

2.1 参数初始化

data.h这个文件中给出了许多初始值:

const int buffer=256;

// encode type
char encodeTypeInit[buffer]="cksaap";
char *encodeType=encodeTypeInit;

// input file name
char *inputFile;

// output file name
char outputFileInit[buffer]="output.txt";
char *outputFile=outputFileInit;
int outputFormat=0;

// pssm file directory
char pssmFileDirectoryInit[buffer]="PSSM";
char *pssmFileDirectory=pssmFileDirectoryInit;

// disorder file directory
char disorderFileDirectoryInit[buffer]="DisorderVSL2";
char *disorderFileDirectory=disorderFileDirectoryInit;

// aggregation file directory
char aggFileDirectoryInit[buffer]="Aggregation";
char *aggFileDirectory=aggFileDirectoryInit;

// AAindex parameter file
char aaindexFileInit[buffer]="AAindex.txt";
char *aaindexFile=aaindexFileInit;
char ifFeatureSelectionInit[buffer]="N";
char *ifFeatureSelection=ifFeatureSelectionInit;
char selectedFeatureFileInit[buffer]="SelectedAAindexFeature.txt";
char *selectedFeatureFile=selectedFeatureFileInit;

// KNN
char knnTrainFileInit[buffer]="train.txt";
char *knnTrainFile=knnTrainFileInit;
char topKFileInit[buffer]="topKValues.txt";
char *topKFile=topKFileInit;

以及截取蛋白质序列片段的参数(seqLen)和窗口值(window)

int seqLen=51;
int window=27;

这些参数中的大部分可以从程序入口函数修改设定,encode.cc这个c++文件是程序的入口。

2.2 参数输入

以我们的项目来说,

./encode -i <chen’s file format> -o ../t3_1.cksaap -t cksaap

这句Linux命令,提供了inputoutputencode type,而下面的语句则是判断inputoutputencode type的作用。

// input
if(strcmp(argv[i], "-i")==0){
    inputFile=argv[i+1];
}

// output
if(strcmp(argv[i], "-o")==0){
    outputFile=argv[i+1];
}

// encode type
if(strcmp(argv[i], "-t")==0){
    encodeType=argv[i+1];
}

而源代码里面还有判断其他选项的模块,诸如output format之类。当然,seqLenwindow的值也可以在命令行中直接修改为自己所需要的值。

// Fragment and encode window
if(strcmp(argv[i], "-L") == 0){
    seqLen=atoi(argv[i+1]);
}
if(strcmp(argv[i], "-W") == 0){
    window=atoi(argv[i+1]);
}

输入文件input加载进内存,以我们的文件T3toChen.txt的第一行内容为例,

ifstream icin(inputFile);
if(!icin){
    printMessage("Can not open the input file!");
}
string tmpStr;
while(!icin.eof()){//如果没有到文件结尾,则执行循环体
    getline(icin, tmpStr);//一行一行读取,相当于Java中的readline()
    if(!tmpStr.empty()){
        string str1, str2, str3;
        int tmpTag=0;//<chen’s file format>每一行均由三个string和一个int组成,这四个值用于存储这四个量
        istringstream iStream(tmpStr);
        iStream>>str1>>str2>>str3>>tmpTag;
        /*
        str1=="LEEYRKHVAERAAEGIAPKPLDANQMAALVELLKNPPAGEEEFLLDLLTNRVPPG
                VDEAAYVKAGFLAAIAKGEAKSPLLTPEKAIELLGTMQGGYNIHP"
        str2=="56479606"
        str3=="any"
        tmpTag==-1
        */
        str3.replace(0, 1, "");//str3=="ny"
        struct prot *p=new prot;
        p->next=NULL;
        p->protSeq=str1;
        p->protName=str2;
        p->position=atoi(str3.c_str());//p->position==0;
        p->tag=tmpTag=-1;
        if(head==NULL && tail==NULL){
            head=tail=p;
        }
        else{
            tail->next=p;
            tail=p;
        }
    }
}
icin.close();
icin.clear();

由于encodeType==cksaap

if(strcmp(encodeType, "cksaap") == 0){
    CKSAAPEncode_2(head, outputFile, ifFeatureSelection, selectedFeatureFile, outputFormat, seqLen, window);
}

故调用encordFunctions.h文件中的函数

void CKSAAPEncode_2(struct prot *data, char *sbjFile, char *featureSelection, char *selectedFeature, int outFormat, int fragLen, int encodeWindow){
    //函数体(较长),主要作用是统计残基对并求出出现的概率。
}

这几个参数中,

 1. data=head=p;
 2. sbjFile=outputFile=t3_1.cksaap;
 3. featureSelection=ifFeatureSelection="N";
 4. selectedFeature="SelectedAAindexFeature.txt";
 5. outFormat=0;
 6. fragLen=51;
 7. encodeWindow=27.

2.3 特征计算

char AA[]={'A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y'};

map<char, int> mymap;
for(int i=0; i<20; i++){
    mymap.insert(pair<char, int>(AA[i], i));
}

struct prot *head=data;
if(strcmp(featureSelection, "N") == 0){//满足条件
    if(outFormat==0){//满足条件
        while(head != NULL){//满足条件
            int m=1;
            if(head->tag == 1){
                ofssub<<"+1  ";
            }
            else{//满足条件
                ofssub<<"-1  ";
            }
            string seq=head->protSeq.substr((int)((fragLen/2)-int(encodeWindow/2)), encodeWindow);
            //这一句是从何处截取蛋白质序列片段的核心语句。
            /*
            对于本例子来讲,从第12个氨基酸开始,截取27个。即
            seq="AEGIAPKPLDANQMAALVELLKNPPAG".
            */
            for(int gap=0; gap<=5; gap++){
                float Num[400];//400维的数组,存储特征残基对。
                int sum=0;//存储残基对个数。
                for(int i=0; i<400; i++){ Num[i]=0; }

                for(int i=0; i<seq.length()-gap-1; i++){
                    if(seq.at(i)=='-'  || seq.at(i+gap+1)=='-') continue;
                    Num[mymap[seq.at(i)]*20 + mymap[seq.at(i+gap+1)]]++;//将20个氨基酸看做20进制的数,然后配对,构思巧妙。
                    sum++;
                }

                for(int i=0; i<400; i++){
                    ofssub<<m<<":"<<Num[i]/sum<<"  ";//输出特征向量到输出文件。
                    m++;//分量的维度。
                }
            }
            ofssub<<endl;
            head=head->next;
        }
    }
    else if(outFormat == 1){
        ofssub<<"@relation features\n\n";
        for(int i=1; i<=2400; i++){
            ofssub<<"@attribute f"<<i<<" real\n";
        }
        ofssub<<"@attribute play {yes, no}\n\n";
        ofssub<<"@data\n";

        while(head != NULL){
            string seq=head->protSeq.substr((int)((fragLen/2)-int(encodeWindow/2)), encodeWindow);
            for(int gap=0; gap<=5; gap++){
                float Num[400];
                int sum=0;
                for(int i=0; i<400; i++){ Num[i]=0; }

                for(int i=0; i<seq.length()-gap-1; i++){
                    if(seq.at(i)=='-' || seq.at(i+gap+1)=='-') continue;
                    Num[mymap[seq.at(i)]*20 + mymap[seq.at(i+gap+1)]]++;
                    sum++;
                }

                for(int i=0; i<400; i++){
                    ofssub<<Num[i]/sum<<",";
                }
            }
            if(head->tag == 1){
                ofssub<<"yes\n";
            }
            else{
                ofssub<<"no\n";
            }
            head=head->next;
        }
    }
    else{
        printMessage("Incorrect output file format.");
    }
}

这样,经过6次循环,这个序列就产生了一个特征向量,维度为2400,这个特征向量输出到t3_1.cksaap这个文件当中。
根据各种不同的需要,可以在命令行中或者直接在data.h文件中改变参数的初始值,以达到不同的效果。
这个算法更多细节,请点击cksaap_ubsite

由于水平和时间有限,难免存在理解的偏差及疏漏,不足之处请联系我们。