您正在查看: 生物信息 分类下的文章

python计算smoothed PSSM(二)

Overview

上一篇文章python计算smoothed PSSM(一)当中,介绍了以当前氨基酸残基为基点,左右取相同数目的序列,然后叠加计算。Chris介绍,这样的算法有特定的用场:蛋白质后修饰。但是,普通的蛋白质序列提取特征就不太适用了:因为窗口值(smoothed window)只能取奇数,而如果有偶数长度的序列片段包含有特征,这种算法就会漏掉。于是决定写一个新的python脚本,把所有特征全部包含进去。
想法很简单:以当前残基为基点,直接向后连续取w_smth个氨基酸,并叠加计算最后存入新矩阵的当前位置。我做了一个循环,将窗口值不大于w_smth的矩阵全部存入一个相同的新矩阵当中,这样特征就全了。

1 python编码

1.1 t34pssm.py

这部分代码跟前面的代码只有一处不同:不用判断窗口值是否为奇数。这部分内容可参考上一篇内容。

1.2 pssm_smoothed_head.py

这部分代码完成的功能如下:

  1. 将每条序列的pssm矩阵的左半部分截取,存入矩阵PSSM-orig
  2. 对矩阵PSSM_orig进行叠加操作,生成矩阵PSSM_smth_head_full
  3. 根据需要截取PSSM_smth_head_full的前n个序列,并存入PSSM_smth_head_final
  4. PSSM_smth_head_final合并为一行写入文件,每条序列占一行。

代码如下:

    #! /usr/bin/env python
    # -*- coding: utf-8 -*-
    # vim:fenc=utf-8
    
    """
    Retrieve smoothed_head PSSM features 
    """
    import sys
    import numpy as np
    import math
    import re
    import fileinput 

    def pssm_smth_head_n(fi,output_smth_head,w_smth_head,n):
        # 0-19 represents amino acid 'ARNDCQEGHILKMFPSTWYV'
        w_smth_head=int(w_smth_head)
        n=int(n)
        Amino_vec = "ARNDCQEGHILKMFPSTWYV"

        PSSM = []
        PSSM_orig = []  
        seq_cn = 0

        # 读取pssm文件
        for line, strin in enumerate(fileinput.input(fi)):
        
            if line > 2:
                str_vec = strin.split()[1:22]
                if len(str_vec) == 0:
                    break
                PSSM.append(map(int, str_vec[1:]))
                seq_cn += 1
        print seq_cn            
        fileinput.close()

        #original PSSM
        #将每条序列的`pssm`矩阵的左半部分截取,存入矩阵`PSSM-orig`
        PSSM_orig=np.array(PSSM)
        #print PSSM_orig
    
        PSSM_smth_head_final=np.array([[0.0]*20]*(n*w_smth_head))

        #section for PSSM_smth_head features
        for k in range(1,w_smth_head+1):
            PSSM_smth_head = np.array([[0.0]*20]*seq_cn)
            #print PSSM_smth_head

            #对矩阵`PSSM_orig`进行叠加操作,生成矩阵`PSSM_smth_head_full`。                                                                                                            
            PSSM_smth_head_full=pssm_smth_head(PSSM_orig,PSSM_smth_head,k,seq_cn)    
            #print PSSM_smth_head_full
            #print np.shape(PSSM_smth_head_full)

            #根据需要截取`PSSM_smth_head_full`的前`n`个序列,并存入`PSSM_smth_head_final`。
            for i in range(n):
                PSSM_smth_head_final[i+n*(k-1)]=PSSM_smth_head_full[i]          
        #print PSSM_smth_head_final
        PSSM_smth_head_final_shp=np.shape(PSSM_smth_head_final)     
        file_out_smth_head=file(output_smth_head,'a')

        #将`PSSM_smth_head_final`合并为一行写入文件,每条序列占一行
        np.savetxt(file_out_smth_head, [np.reshape(PSSM_smth_head_final, (PSSM_smth_head_final_shp[0] * PSSM_smth_head_final_shp[1], ))], delimiter=",")

    def pssm_smth_head(PSSM_orig,PSSM_smth_head,w_smth_head,l):
        for i in range(l):
            if i <=l-w_smth_head:
                for j in range(i,i+w_smth_head):
                            
                    PSSM_smth_head[i]+=PSSM_orig[j]             
            else:
                for j in range(i,l):
                    PSSM_smth_head[i]+=PSSM_orig[j]     
        return PSSM_smth_head

1.3 总结

这个程序可以得到比较全的特征。比如取窗口值为10,那么窗口值为1,2,3,4,5,6,7,8,9,10的特征将会全部被包含在内,并合成一行。
linux中进入文件所在的目录,然后终端运行如下命令:

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

需要保存的结果文件,自己修改。

python计算smoothed PSSM(一)

Overview

最近几天,Chris和我看了很多论文,对PSSM有了更深的认识。但是,鉴于PSSM本身包含单个位置的信息更明显,而几乎没有包含蛋白质序列片段信息,我们两人思考如何将蛋白质序列片段信息编码,终于找到了一种PSSM的处理方式,这种方式叫做smoothed window,特此记录一下。
该算法原理,请参考这篇论文:Predicting RNA-binding sites of proteins using support vector machines and evolutionary information,并在此感谢该论文作者!并感谢Chris对我的鼓励和帮助!

1 python编码

1.1 t34pssm.py

这部分代码是主函数,是之前这篇文章中的代码,可以参考这里蛋白质序列特征提取方法之——PSSM,就不详加解释了。

    #! /usr/bin/env python
    # -*- coding: utf-8 -*-
    # vim:fenc=utf-8
    """
        python t34pssm.py total_train_60.fasta ./total_train_60 ./total_train_60_pssm  w_smth n 
        param:
            1. 总的fasta格式蛋白质序列
            2. 分开的fasta格式蛋白质序列的文件夹 
            3. 分开的fasta格式蛋白质序列对应的pssm的文件夹  
            4.smooth-window 值,要求是奇数。
            5.截取的序列长度,一般为25,30,50。本例为30.
    """
    import fileinput
    import sys
    from os import listdir
    from os.path import isfile, join
    import re
    from pssm_smoothed import *
    
    smplfasta = sys.argv[1]  
    spfasta = sys.argv[2]   
    check_head = re.compile(r'\>')  
    
    #read from undersample fasta, store 
    smplist = []
    smpcnt = 0
    for line, strin in enumerate(fileinput.input(smplfasta)):
        if check_head.match(strin):
            smplist.append(strin.strip())
            smpcnt += 1
    onlyfiles = [ f for f in listdir(spfasta) if isfile(join(spfasta,f)) ]

    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

    pssmdir = sys.argv[3]

    w_smth = sys.argv[4]
    #如果窗口值不是奇数,退出程序
    if int(w_smth)%2 ==0:
        print 'Please change your input argument ' + w_smth + ' to an odd smoothing-window number!!!'
        sys.exit(1)
    n=sys.argv[5]

    for smp in smplist:
        finalist.append(pssmdir+'/'+fastaDict[smp].split('.')[0]+'.pssm')

    for fi in finalist: 
        pssm_single(fi,'total_train_60_pssm_smth',w_smth,n)

1.2 pssm_smoothed.py

    #! /usr/bin/env python
    # -*- coding: utf-8 -*-
    # vim:fenc=utf-8
    """
    Retrieve smoothed PSSM features 
    """
    import sys
    import numpy as np
    import math
    import re
    import fileinput 
    def pssm_single(fi,output_smth,w_smth,n):
        # 0-19 represents amino acid 'ARNDCQEGHILKMFPSTWYV'
        w_smth=int(w_smth)
        n=int(n)
        Amino_vec = "ARNDCQEGHILKMFPSTWYV"

        PSSM = []
        PSSM_orig = []  
        seq_cn = 0 
        # 读取pssm文件
        for line, strin in enumerate(fileinput.input(fi)):
            if line > 2:
            str_vec = strin.split()[1:22]
            if len(str_vec) == 0:
                break
            PSSM.append(map(int, str_vec[1:]))
            seq_cn += 1
        print seq_cn            
        fileinput.close()
        PSSM_smth = np.array([[0.0]*20]*seq_cn)
        #print PSSM_smth
        PSSM_orig=np.array(PSSM)
        #print PSSM_orig
        #section for PSSM_smth features
        PSSM_smth_full=pssm_smth(PSSM_orig,PSSM_smth,w_smth,seq_cn)
        PSSM_smth_final=[[0.0]*20]*n

        #截取PSSM_smth_full矩阵的前n行,作为输出内容
        for i in range(n):
            PSSM_smth_final[i]=PSSM_smth_full[i]
            #print PSSM_smth_final[i]
        PSSM_smth_final_shp=np.shape(PSSM_smth_final)
        # for i in range(seq_cn):
        #   print PSSM_smth_final[i]
        file_out_smth=file(output_smth,'a')
        np.savetxt(file_out_smth, [np.reshape(PSSM_smth_final, (PSSM_smth_final_shp[0] * PSSM_smth_final_shp[1], ))], delimiter=",")

    #这个函数会求出整条序列的smoothed pssm矩阵
    def pssm_smth(PSSM_orig,PSSM_smth,w_smth,l):
        for i in range(l):
            #smooth窗口超过pssm上边界
            if i <(w_smth-1)/2:
                for j in range(i+(w_smth-1)/2+1):
                    #print i,j              
                    PSSM_smth[i]+=PSSM_orig[j]
                    #print PSSM_smth[i]
            #smooth窗口超过pssm下边界
            elif i>=(l-(w_smth-1)/2):
                for j in range(i-(w_smth-1)/2,l):   
                    #print i,j          
                    PSSM_smth[i]+=PSSM_orig[j]
                    #print PSSM_smth[i]
            else:
                for j in range(i-(w_smth-1)/2,i+(w_smth-1)/2+1):
                    #print i,j
                    PSSM_smth[i]+=PSSM_orig[j]
                    #print PSSM_smth[i]     
        return PSSM_smth

1.3 总结

该算法以行为单位进行运算。以本行为中心,上下扩展,扩展的上下长度为smooth窗口值,将这些行的值加起来,存入新的矩阵的相同行位置。这样新的矩阵就包含了多个连续氨基酸序列片段的信息,会为特征提取提供新的思路。

python分离正负样本

Overview

机器学习算法在项目中应用的时候,有时候会用到分离fasta格式的正负样本。于是就写了一个python脚本,效果不错,记录一下。

1. separatePosNeg.py

    #! /usr/bin/env python
    # -*- coding: utf-8 -*-
    # vim:fenc=utf-8
    import fileinput
    import sys
    import re

    ls=fileinput.input(sys.argv[1])
    check_head = re.compile(r'\>')
    cls = ''
    smplist = ''
    for line, strin in enumerate(ls):
        if check_head.match(strin):
            strin_vec=strin.split('|')
            cls = strin_vec[len(strin_vec)-1].strip()
            if cls == str(sys.argv[2]):
                smplist+=strin.strip()
                smplist+='\n'
                smplist+=ls.next()
    print smplist

记录下在linux命令行中的用法:

    #提取正样本时:
    python separatePosNeg.py original.fasta '1' >> pos.fasta
    #提取负样本时:
    python separatePosNeg.py original.fasta '-1' >> neg.fasta

Ubuntu 14.04 安装R和R packages

Overview

虽然已经用了很久的R语言,但一直没整理过,正好需要在我们的云服务器上安装R,所以一并记录下来了。下面的过程虽然是在Ubuntu 14.04上安装的,但是对于其他版本的系统,RR packages的安装都大同小异。

1. 安装R

1.1 添加源

Ubuntu 14.04中的R版本比较旧,默认安装可能会出很多问题(我试过了)。所以最好添加一个新的源。

etc/apt/source.list中添加如下信息:

deb http://<my.favorite.cran.mirror>/bin/linux/ubuntu trusty/

这里的<my.favorite.cran.mirror>换成一个适合你的镜像地址,可以从https://cran.r-project.org/mirrors.html找到所有的镜像,比如我选择了厦门大学的镜像:http://mirrors.xmu.edu.cn/CRAN/,那么在etc/apt/source.list中添加:

deb http://mirrors.xmu.edu.cn/CRAN/bin/linux/ubuntu trusty/

1.2 更新源

sudo apt-get update

不出意外会遇到一个错误:

W: GPG 错误:http://mirrors.xmu.edu.cn trusty/ Release: 由于没有公钥,无法验证下列签名: NO_PUBKEY 51716619E084DAB9

运行下面的命令:

sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys 51716619E084DAB9

之后重新更新源:

sudo apt-get update

就可以顺利更新了。

1.3 安装

sudo apt-get install r-base

如果要从源码安装R,使用下面的命令:

sudo apt-get install r-base-dev

2. 安装 R packages

直接在命令行输入R,就可以进入R语言的命令行交互模式,使用几个命令可以快速熟悉一下R。

  1. 如果你想查看一个函数的用法,比如mean函数(R中用来求均值的函数),可以使用

     help(mean)
    

    或者

     ?mean
    
  2. .libPaths(),查看R中包安装的目录,在我的Ubuntu14.04下,有三个输出下面的内容:

     [1] "/usr/local/lib/R/site-library" "/usr/lib/R/site-library"
     [3] "/usr/lib/R/library"
    

    说明上面三个目录都可以作为R packages安装的目录,其中,如果你在安装包时不指定,默认会安装在第一个目录/usr/local/lib/R/site-library,而/usr/lib/R/libraryR自带包的所在目录。

  3. library(),可以查看R中以及安装的所有包,这个函数会去.libPaths()中列出的所有目录查看,并按照目录将各种包列出来:

     Packages in library ‘/usr/local/lib/R/site-library’:
    
     digest                  Create Compact Hash Digests of R Objects
     entropy                 Estimation of Entropy, Mutual Information and Related Quantities
     FSelector               Selecting attributes
     randomForest            Breiman and Cutler's Random Forests for Classification and Regression
     RWeka                   R/Weka interface
     RWekajars               R/Weka interface jars
     Packages in library ‘/usr/lib/R/site-library’:
     rJava                   Low-level R to Java interface
     Packages in library ‘/usr/lib/R/library’:
     base                    The R Base Package
     boot                    Bootstrap Functions (originally by Angelo Canty for S)
     class                   Functions for Classification
     cluster                 Cluster Analysis Extended Rousseeuw et al.
     codetools               Code Analysis Tools for R
     compiler                The R Compiler Package
     ...
    

    这里省略了一些包,并没有全部列出来,你可能很好奇,rJava这个包为什么会比较特别,安装在了一个其他包都不在的地方。下面在介绍FSelector包安装的时候会讲到。

  4. install.packages("mypackage"),用来安装的函数,默认会安装在.libPaths()列出来的第一个目录下。这个函数也提供了参数lib用来指定你要把包安装的位置。

  5. library("mypackage"),载入一个已经安装的包,如果还没有安装这个包,则会报错,程序中断。

  6. require("mypackage"),载入一个已经安装的包,如果还没有安装这个包,则显示警告信息,然后程序继续向下执行。(与library("mypackage")相比,require强烈不推荐)。

  7. remove.packages("mypackage"),卸载package。

  8. update.packages( ),更新package,可以定期执行。

  1. 安装Bioconductor 包

    Bioconductor提供了一个函数用来安装Bioconductor 包

    • 首先,载入这个安装函数:

       source("http://bioconductor.org/biocLite.R")
      
    • 之后,就可以用这个函数安装 Bioconductor 包了:

       biocLite("mypackage")
      
  2. 下面还有一些方法,不太常用,列出来参考:

    下面的内容主要摘自:【R】Linux安装R语言包(Installing R packages on Linux)

    • getOption("defaultPackages"):查看启动R时自动载入的包。
    • vignette('mypackage'):有的包,特别是bioconductor的包有vignette,用函数查看
    • openVignette('mypackage'):这个函数也可以查看vignette,更好用一些
    • RSiteSearch("helpinfor"):搜索R网站上的“helpinfor”相关信息
    • help.start():查看已经安装包的详细HTML文档
    • search():查看当前载入的包
    • sessionInfo():查看R中载入的包
    • methods():查看某个S3泛型函数中所有的方法或者一个类中所有的方法(S3:S version 3)
    • showMethods(class = "myClass"):查看S4类的方法
    • findMethods("myMethods"):查看method的代码
    • class(myObject):查看某个对象的类
    • getClass(“class/package”):查看某个class或者包的具体内容
    • getSlots("class"):查看某个class的slot
    • slotNames(MyObject):查看某个对象的slot。
    • Myobject@slotNames:访问对象Myobject的slot值,这个@可以连续用。

    查询包内信息:

    • ?function/method:查看某个“函数”或者“方法”的详细内容
    • class?graph::graph:查看“组”的详细内容的一个例子。这个例子的来源是查询graph包时候,查看其中class的信息,输入??graph后出现一个graph::graph-class。
    • ls("package:mypackage"):查看"mypackage"中的所有对象。

参考资料

生物信息中常用的Linux命令

Overview

一直想把常用的命令搜集起来,以便平时用到的时候查阅,可惜一直没抽出来时间专门整理下。最近在做序列的特征提取和多个特征文件合并时,频繁使用到了一些命令,干脆一并整理到这里,以后边用边添加整理新的命令。

这里的linux命令主要在MAC 10.9.5Ubuntu 14.04下测试,涉及到平台差异性的时候,会尽量指出来,没有区分的话就表示两种平台下都可以使用。如果仍有没涉及到的问题,欢迎补充~

1. 常规的shell命令

1.1 查看文件的前/后n行

  • Linux

    下面的命令显示文件名为filename的文件的前/后5行:

      head filename -n5
      tail filename -n5 
    

    设置-n后面的数字,可以指定要查看的行数。如果省略-n参数,则默认显示10行。
    另外,-n参数的位置是比较灵活的,而且后面的数值可以加一个空格区分,所以,下面几种形式都是可以的:

      head filename
      head -n5 filename
      tail -n5 filename
      head filename -n 5
      tail filename -n 5
    
  • Mac
    Mac下的head/tail命令是BSD平台的实现,-n参数必须在文件名前面,放在文件名后面会被当做一个新的文件名。如果省略-n参数,则默认也显示10行。
    下面几种形式都是可以的:

      head -n5 filename
      tail -n5 filename
      head -n 5 filename
      tail -n 5 filename
    

1.2 查看文件的行数

wc -l filename

上面的命令会把把名为filename文件的行数以及文件名本身一起列出来。有时候我们想要同时看很多个文件的行数。可以使用:

wc -l filename1 filename2

这条命令输出为:

281 filename1
281 filename2
562 total

分别为filename1,filename2的行数,以及他们行数的和。当然,也可以直接查看当前文件夹下所有的文件行数:

wc -l *

1.3 查看文件的列数

在提取完特征之后,检查特征文件的列数是必须要做的,所以查看文件列数的命令也非常常见:

cat filename | awk -F ',' '{print NF}' | head -n1

这里假定filename是以,分割的文件,如果是以空格或者别的标记分割的文件,可以自行修改
','中的符号。

另外,awk -F ',' '{print NF}'中输出的结果,实际上是每一行的列数,我们假定这里处理的是特征文件,因此每一行的列数都是相同的,所以我们只需要使用head -n1查看第一行的列数就可以了。

1.4 查看序列文件中的序列数

grep -c '>' filename

因为序列文件中每条序列都是以>开头的,因此统计>的个数就可以知道序列数目。

1.5 文件合并

  • 上下合并:

      cat file1 file2 > file3
      cat file1 file2 file3 > newfile
    

    很容易扩展到多个文件合并

  • 左右合并

      paste file1 file2 > file3
      paste file1 file2 file3 > newfile
    

    一样可以扩展到多个文件合并。

    默认情况下,paste合并文件之间使用空格或者tab分开,如果你合并的是csv文件,就需要显示指定用,分开:

      paste -d "," file1 file2 file3 > newfile
    

1.6 单个文件去掉重复的行

  • 重复的多行只留一行

      sort filename | uniq > newfile
    
  • 重复的行全部去掉,只留下文件中的非重复行

      sort filename | uniq -u > newfile
    

1.7 文件的交集,并集

下面的操作很容易可以扩展到多个文件。

  • 文件的并集(重复的行只保留一份)

      cat file1 file2 | sort | uniq > newfile
      cat file1 file2 file3 | sort | uniq > newfile
    
  • 取出两个文件的交集(只留下同时存在于两个文件中的文件)

      cat file1 file2 | sort | uniq -d > newfile
      cat file1 file2 file3 | sort | uniq -d > newfile
    
  • 删除交集,留下其他的行

      cat file1 file2 | sort | uniq -u > newfile
      cat file1 file2 file3 | sort | uniq -u > newfile
    

1.7 查看及关闭端口

通常在tomcat没能正常退出的情况下,8080端口就会被占用,这时就无法重新启动tomcat。因此经常需要查看端口使用情况,并在端口被占用时关闭端口。

Mac

  • 查看端口8080

      lsof -i:8080
    
  • 关闭8080

      kill -9 8080
    

参考资料