您正在查看: 标签 bioinformatics 下的文章

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
    

参考资料

配置Hibernate解决MYSQL连接失效问题

Overview

之前将SecretEPDB部署到了云服务器上之后,再打开需要连接数据库的网页时总是会出现莫名其妙的错误,之前一直没管它,主要是因为这个错误不是每次都出现,出现之后刷新几次又可以访问了。

1. 错误描述

每次打开需要连接数据库的网页,就很有很大概率出现下面的错误信息:

Struts Problem Report

Struts has detected an unhandled exception:

Messages:   
Broken pipe
The last packet successfully received from the server was 144,004,781 milliseconds ago. The last packet sent successfully to the server was 144,004,781 milliseconds ago. is longer than the server configured value of 'wait_timeout'. You should consider either expiring and/or testing connection validity before use in your application, increasing the server configured values for client timeouts, or using the Connector/J connection property 'autoReconnect=true' to avoid this problem.
No operations allowed after connection closed.
could not execute query

第一反应是Struts的问题,但很快就可以排除,因为:

  1. 在本地开发时,不会有这个问题,只有部署到服务器上了才会出现。
  2. 错误信息很明确,这个服务器链接失效了。

2. 分析

查了一些资料,定位了这个错误,是因为MYSQL的链接失效问题,在 MySQL相关问题 中,我们配置了wait_timeout=31536000这句话,所以MYSQL31536000毫秒之后会失效,尽管我们已经将wait_timeout的值改得很大了(默认是8小时),但是即使再大,也会有失效的时候,这时候服务器会得到一个失效的数据库连接,但是Hibernate并不知道它失效了,继续交给Struts使用,就会报上面的错误。

重启tomcat可以解决这个问题,因为所有的连接都重新初始化了嘛,但是时间久了肯定还会出问题,这也是为什么在本地开发时不会出问题,因为tomcat总在不停地重启,连接就没等到失效的时候。

3. 解决方法

查了一下原因,为了偷懒,选择了最简单的解决方法,在hibernate配置文件hibernate.cfg.xml中添加下面的配置:

<property name="connection.autoReconnect">true</property> 
<property name="connection.autoReconnectForPools">true</property> 
<property name="connection.is-connection-validation-required">true</property> 

第二天发现又有问题了,看来这样不是很可行。

3.1 使用c3p0数据库连接池

其实也不是很麻烦,Hibernate默认就支持c3p0,先去你使用的Hibernate版本的文件包的lib文件夹中,找到c3p0-***.jar,版本号与当前的Hibernate有关。从Hibernate包中去找是一个最佳实践,可以让你避免一些不兼容性问题。我们使用的是Hibernate 3,所以我去下了一个Hibernate3的包,官网上已经不再提供Hibernate 4以下的包了,可以去 sourceforge上下载这个版本。

下载时候,将lib文件夹里面的c3p0-***.jar加入项目的引用,并在hibernate.cfg.xml中添加下面的配置(其实默认hibernate.cfg.xml中有这部分代码,只是是注释掉的,取消注释就可以):

<!-- configuration pool via c3p0-->      
<property name="hibernate.connection.provider_class">org.hibernate.connection.C3P0ConnectionProvider</property> 
<property name="c3p0.min_size">5</property> 
<property name="c3p0.max_size">30</property> 
<property name="c3p0.time_out">1800</property> 
<property name="c3p0.max_statement">50</property> 
<property name="c3p0.acquire_increment">1</property> 
<property name="c3p0.idle_test_period">120</property> 
<property name="c3p0.validate">true</property> 

重新发布项目,在服务器上部署以后,不仅解决了链接失效问题,数据库访问速度快了很多。因为c3p0维护了数据库连接池,并在使用时验证链接是否可用(<property name="c3p0.validate">true</property>),如果失效了,则重新打开。

4. 堆栈错误补充

还有一种情况可以导致MySQL链接失效,错误信息如下:

… Stacktraces
org.hibernate.exception.JDBCConnectionException: could not execute query
…
com.mysql.jdbc.exceptions.jdbc4.MySQLNonTransientConnectionException: No operations allowed after connection closed.
…
java.lang.OutOfMemoryError: Java heap space
… Struts Problem Report

Struts has detected an unhandled exception:

Messages:   
Java heap space
No operations allowed after connection closed.
could not execute query

我将这个错误分成两条显示,是为了表示下面的部分才是核心问题所在:由于堆栈溢出,导致java虚拟机崩溃,从而数据库连接中断,而且之后每次查询都会报数据库连接错误,根本原因是堆栈溢出。好了,知道了问题所在,就容易解决了:更改java虚拟机内存。
找到tomcatbin目录下的catalina.sh,打开它 (以tomcat7为例,如果是用apt-get命令安装的tomcat,这个bin目录在/usr/share/tomcat7下面)。
cygwin=false之前,加上下面这句话:

JAVA_OPTS="-Xmx512m"

注意:双引号要加上。如果512MB还是不行,可以加大。

问题起源:我们在secretEPDB项目线上版本测试时,发现查询第1962条数据时,查询总是中断,之后整个程序都会崩溃,即使第一次就查询第1962条数据仍会程序全部崩溃,我们推断,这不是数据库连接池的问题,应该是堆栈过小导致虚拟机崩溃的问题。结果和我们的推测一样,这条数据特别大,而我们的服务器默认的虚拟机内存为128MB,改为512MB之后,程序一切正常。

5. 参考资料

SCRATCH的预测结果格式

Overview

最近在预测蛋白质序列的二级结构,结构性区域,水溶性等特征时,使用了不同的软件,发现不同软件预测结果中对同一特征的表示方式略有不同,所以在这里一并总结。

1. SCRATCH中的输出格式

我们在 SCRATCH的安装和使用 介绍了SCRATCH的安装和使用,直接使用

./run_SCRATCH-1D_predictors.sh input_fasta  out_prefix  [num_threads]

就可以得到四个结果,分别是SSpro(后缀为.ss)SSpro8(后缀为.ss8)ACCpro(后缀为.acc)ACCpro20(后缀为.acc20)

以下面一段氨基酸为例:

MQIFVKTLTGKTITLEVEPSDTIENVKAKI

产生的结果格式如下:

SSpro预测产生的二级结构序列(.ss文件):
CEEEEEEECCCEEEEEECCCCCHHHHHCCC

SSpro8预测产生的二级结构序列 (有8个种类,.ss8文件):
CEEEEEEEESEEEEEEECCCSHHHHEECCC

ACCpro预测产生的水溶性 (exposed threshold为25%,.acc文件):
ee---ee-eeee-e-e-eeeee-ee-eeee

ACCpro20预测产生的水溶性 (.acc20文件):
0%   eeeeeeeeeeeeeeeeeeeeeeeeeeeeee
5%   eeeeeeeeeeeeeeeeeeeeeeeeeeeeee
10%  eeeeeeeeeeeeeeeeeeeeeeeeeeeeee
15%  eee--eeeeeee-e-eeeeeeeeeeeeeee
20%  eee--ee-eeee-e-eeeeeeeeee-eeee
25%  eee--ee-eeee-e-e-eeeeeeee-eeee
30%  ee---ee-eeee-e-e-eeeee-ee-eeee
35%  ee---ee-eeee-e-e-eeeee-ee-eeee
40%  ee---ee-eeee-e-e-eeeee-ee-eeee
45%  ee---e--eee----e---ee--ee-eeee
50%  ee--------e--------e---ee-eeee
55%  e----------------------e---eee
60%  e--------------------------eee
65%  e---------------------------ee
70%  -----------------------------e
75%  ----------------------------e
80%  -----------------------------e
85%  -----------------------------e
90%  ------------------------------
95%  ------------------------------
100% ------------------------------

上面是 官方的帮助文档 给出的输出格式,实际上我算出来的结果.acc20文件与这个不一样,其他都一致。

2. SCRATCH输出格式说明

这里重点说明SSpro(后缀为.ss)SSpro8(后缀为.ss8)ACCpro(后缀为.acc)这三个,以及 DISpro(没在SCRATCH安装包中,需要单独装)的输出格式,并列出了其他同类软件的不同表示。

2.1 SSpro(后缀为.ss)

输出的蛋白质二级序列有三种类别:

  • H = helix
  • E = strand
  • C = the rest

PSIPRED(参见 PSIPRED的安装和使用 )产生的二级序列类别一致。

2.2 SSpro8(后缀为.ss8)

输出的蛋白质二级序列有八种类别:

  • H: alpha-helix
  • G: 310-helix
  • I: pi-helix (extremely rare)
  • E: extended strand
  • B: beta-bridge
  • T: turn
  • S: bend
  • C: the rest

2.3 ACCpro(后缀为.acc)

水溶性的表示如下:

  • - : the residue is buried
  • e : the residue is exposed

也见过如下形式的表示:

  • b : the residue is buried
  • e : the residue is exposed

2.4 DISpro

非结构性区域的表示方式为:

  • O : the residue is ordered
  • D : the residue is disordered

DISOPRED(参见 DISOPRED的安装和使用)中的表示方式为:

  • . : the residue is ordered
  • * : the residue is disordered

3. 参考资料