您正在查看: 2016年1月

python转换文件格式的一处bug

Overview

在转换fasta格式的文件为chen's format文件时,发现前人的程序有些bug,会将最后一条正样本的class标记为-1于是将程序改了一下,这下就没有问题了。

解决方案

fasta格式如下,第一行为头信息,第二行为氨基酸序列:

>sp|Q2YIT7|VIRB3_BRUA2 Type IV secretion system protein virB3 OS=Brucella abortus (strain 2308) GN=virB3 PE=3 SV=1|1
MTTAPQESNARSAGYRGDPIFKGCTRPAMLFGVPVIPLVIVGGSIVLLSVWISMFILPLIVPIVLVMRQITQTDDQMFRLLGLKAQFRLIHFNRTGRFWRASAYSPIAFTKRKRES

现在需要将其转换为如下格式,以供cksaap以及aaindex算法之用:

MTTAPQESNARSAGYRGDPIFKGCTRPAMLFGVPVIPLVIVGGSIVLLSVWISMFILPLIVPIVLVMRQITQTDDQMFRLLGLKAQFRLIHFNRTGRFWRASAYSPIAFTKRKRES    Q2YIT7  any 1

修改后的python代码如下:

import sys
import re
import fileinput  

check_head = re.compile(r'\>')
seqs = ''
head_line = ''
cls = ''
for line, strin in enumerate(fileinput.input(sys.argv[1])):

    if check_head.match(strin):
        strin_vec=strin.split('|')
        cls = strin_vec[len(strin_vec)-1].strip()
        head_line = strin_vec[1].strip()                
    else:
        seqs = strin.strip()    
        print seqs+'\t'+head_line+'\t'+'any'+'\t'+cls

有些程序处理类似的格式时,有可能会把下一行或者上一行的信息和本行信息混淆,所以处理完之后最好手动检查一下正负样本交界处的数据。

R中randomForest包与ggplot2的一个不兼容问题

Overview

在做实验的时候因为要用到随机森林,所以使用了R中的randomForest包,但在画图的时候报了一个非常诡异的错误。

1. 错误描述

下面是我引入randomForest包之后的代码,这里省略了一些细节,只保留跟错误有关的代码:

    ## randomforest
    library("randomForest")
    # randomforest调参
    bestmtryMatrix <- tuneRF(train.bal[,-1],train.bal[,1], ntreeTry=100, 
     stepFactor=1.5,improve=0.01, trace=TRUE,plot=FALSE, dobest=FALSE)
    bestmtryMatrix
    # 找到最佳参数
    bestmtry <- bestmtryMatrix[which.min(bestmtryMatrix[,2]),1]
    # 训练模型
    rfNormal <-randomForest(Class~.,data=train.bal, mtry=bestmtry, ntree=1000, 
     keep.forest=TRUE, importance=TRUE)
    # 预测
    rfNormal.predictions <- predict(rfNormal,type="prob",newdata=test)

在画图时的代码:

g <- ggplot(data=rocdata,aes(x=ejeX,y=ejeY))+geom_line(aes(colour=Metodo),size=0.8)+scale_colour_manual(values=c("red","blue","green4"))+xlab("False positive rate")+
ylab("True positive rate")+theme_custom()+theme(legend.position=c(0.73,0.1))+
#设定x,y轴坐标刻度
scale_x_continuous(limits=c(0, 1),breaks=seq(0,1,0.2))+
scale_y_continuous(limits=c(0, 1),breaks=seq(0,1,0.2))
# 添加一条对角线
g <- g+geom_abline(intercept=0,slope=1) 

上面代码中,theme_custom()就是我自定义的样式函数,可以在 使用R语言ggplot2包画图时的一个不兼容问题 中查看这个函数的详细信息。

在运行时报了下面的错误:

错误于UseMethod("margin") : "margin"没有适用于"unit"目标对象的方法
Calls: theme_custom ... inherits -> theme -> element_text -> structure -> margin
停止执行

又是element_textmargin的问题,可见更新之后这里有点不稳定,我们在 使用R语言ggplot2包画图时的一个不兼容问题 讨论的问题也与margin相关。

2. 解决方法

刚开始我以为是这个函数名在randomForest中有同名函数导致冲突,所以修改了theme_custom(),改为theme_custom_chris(),结果依然报错。但是我试着删除了randomForest的相关代码(library("randomForest")这句也一定要删掉),则不再报错。

实际上,正是library("randomForest")引入randomForest包导致的问题,那在使用完randomForest之后,手动将它从当前的程序中去掉是不是就可以了呢?

R中使用完randomForest包的代码之后,ggplot画图代码之前,添加下面的代码:

detach("package:randomForest", unload=TRUE)

意义也非常直观明了,在运行过程中detach randomForest包,然后重新运行程序,一切正常。

使用R语言ggplot2包画图时的一个不兼容问题

Overview

昨天写了个R语言脚本,主要是借助ROCR包画ROC曲线,只是最后画图时,没用ROCR提供的plot函数画ROC曲线,而是用ROCR处理了数据后,提取了画图的数据,用ggplot2包画了ROC曲线。因为ggplot2可以提供强大的自定义绘图功能,没想到也正是这个自定义样式函数,在移植的时候出现了一些兼容性问题。

1. 问题描述

我使用的ggplot2样式函数,是一个我已经用了一年多的样式,专门写在了一个名字叫style.R的R文件中,内容如下:

library(grid)

theme_custom <- function (base_size = 10, base_family = "serif") {
    theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
          line =               element_line(colour = "black", size = 1, linetype = 1, lineend = "butt"),
          rect =               element_rect(fill = "white", colour = "black", size = 1, linetype = 1),
          text =               element_text(family = base_family, face = "plain", colour = "black", size = base_size, hjust = 0.5, vjust = 0.5, angle = 0, lineheight = 0.9),
          axis.text =          element_text(size = rel(0.8), colour = "grey50"),
          strip.text =         element_text(size = rel(0.8)),
          ## used to delete the text up the picture    
          #strip.text = element_blank(),

          axis.line =          element_blank(),
          axis.text.x =        element_text(family = base_family, size = base_size * 1.2, lineheight = 0.8, vjust = 1.2),
          axis.text.y =        element_text(family = base_family, size = base_size * 1.2, lineheight = 0.8, hjust = 1.2),
          axis.ticks =         element_line(colour = "black", size=0.2),
          axis.title.x =       element_text(family = base_family, size = base_size*1.5, vjust = 0.5),
          axis.title.y =       element_text(family = base_family, size = base_size*1.5, angle = 90, vjust = 0.5),
          axis.ticks.length =  unit(0.15, "cm"),
          axis.ticks.margin =  unit(0.1, "cm"),
      
          legend.background =  element_rect (fill="white",colour="grey50", size=0.1, linetype=1),
          legend.margin =      unit(0.5, "cm"),
          legend.key =         element_rect(fill = "grey95", colour = "white"),
          legend.key.size =    unit(1.2, "lines"),
          legend.key.height =  NULL,
          legend.key.width =   NULL,
          legend.text =        element_text(family = base_family, size = base_size * 1.2),
          legend.text.align =  NULL,
          ##legend.title =       element_text(family = base_family, size = base_size * 3, face = "bold", hjust = 1),
          legend.title =element_blank(),
          legend.title.align = NULL,
          legend.position =    c(0.7,0.15),
          legend.direction =   NULL,
          legend.justification = "center",
          legend.box =         "horizontal",

          panel.background =   element_rect(fill = "white", colour = NA),
          panel.border =       element_rect(fill = NA, colour = "grey50"),
          panel.grid.major =   element_line(colour = "grey30", size = 0.5, linetype="dotted"),
          panel.grid.minor =   element_blank(),
          #panel.grid.minor =   element_line(colour = "grey30", size = 0.5, linetype="dotted"),
          ## panel.margin =       unit(c(0.1, 0.1, 0.1, 0.1), "lines"),

          strip.background =   element_rect(fill = NA, colour = NA),
          strip.text.x =       element_text(family = base_family, size = base_size * 0.8),
          strip.text.y =       element_text(family = base_family, size = base_size * 0.8, angle = -90),

          plot.background =    element_rect(colour = NA, fill = "white"),
          plot.title =         element_text(family = base_family, size = base_size),
          ##plot.title = element_blank(),
          plot.margin =        unit(c(0.5, 1, 0.5, 0.5), "lines") #top,right,bottom,left
    )
}

这个样式我在MacUbuntu 12.04上都使用过,都可以正常使用,但是这两天在windowsUbuntu 15.04上使用时都出了问题,报错如下:

警告信息:
`axis.ticks.margin` is deprecated. Please set `margin` property  of `axis.text` instead
错误于FUN("text"[[1L]], ...) :
  Theme element 'text' has NULL property: margin, debug
Calls: print ... element_render -> calc_element -> lapply -> FUN -> lapply -> FUN
停止执行

这其实是两个问题,一个警告,外加一个错误。

第一个警告信息很好理解,我把代码中的

axis.ticks.margin =  unit(0.1, "cm"),

直接注释掉就没有了。但是在windowsUbuntu 15.04上使用时,依然报错:

错误于FUN("text"[[1L]], ...) :
  Theme element 'text' has NULL property: margin, debug
Calls: print ... element_render -> calc_element -> lapply -> FUN -> lapply -> FUN
停止执行

2. 错误复现

首先想到的就是:为什么在MacUbuntu 12.04却没有报错,查了一些资料,终于在ggthemes的github issues 上面找到了原因,原因在于ggplot2总是在不停地更新,而更新之后,对于旧版本的兼容性不是很好,我的MacUbuntu 12.04使用的是旧版本,因此还可以正常使用,但是使用新版本的电脑就不能很好的兼容这段代码。因此我把mac下的ggplot2更新一下。在这里我并没有更新ggplot2,而是更新了ggthemes

2.1 ggthemes

ggthemes包:丰富ggplot2的表现力 这个链接的名字已经解释了ggthemes

ggplot2包的最新0.93版本允许自定义主题,这样ggplot的表现力可以通过各种不同的主题获得提升。
ggthemes包就是ggplot2的主题扩展包,提供了供ggplot2使用的新主题,尺度,几何对象和一些新函数。

其实这一两节是我后面加上的,原因是我google的时候正好在ggthemesgithub issues中找到了错误和复现的解决方法。 然后又理了一下思路,发现其实我想浮现这个问题,应该更新ggplot2包就可以了,单纯更新ggthemes应该会顺带更新它的依赖项,所以ggplot2也顺带更新了。

这里也只是一个逻辑上的分析,标记一下,以后再仔细确认。

2.2 安装devtools

在R交互命令行下,运行下面的命令安装devtools包:

install.packages("devtools")

2.3 从github上更新ggthemes

安装devtools包之后,就可以使用devtools包直接从github上更新最新的ggthemes了,而不是从CRAN上更新。 命令如下:

devtools::install_github("jrnold/ggthemes")

安装之后,重新运行ROC画图代码,谢天谢地,也报了同样的错误。

3. 解决方法

很明显,这个是由于ggplot2版本更新引起的一个兼容问题,在 github ggthemes issue=57可以找到答案:

It does look like the documentation for ggplot2 may not have been updated to account for changes made with ggplot2 2.0.0. ggplot2 version 2.0.0 added the elements margin and debug to element_text, so you need to change the text element of the theme...

新版本的ggplot2中,text = element_text()中,多了两个参数margindebug,我们上面style.R中,并没有使用这两个参数,所以才会报错,把style.R中:

text =               element_text(family = base_family, face = "plain", colour = "black", size = base_size, hjust = 0.5, vjust = 0.5, angle = 0, lineheight = 0.9),

修改为

text =               element_text(family = base_family, face = "plain", colour = "black", size = base_size, hjust = 0.5, vjust = 0.5, angle = 0, lineheight = 0.9,margin = margin(), debug = FALSE),

重新运行画图的R代码,就可以正常绘图了。

R语言学习笔记(一)

写在前面的后记

昨晚写完这一篇的时候,Chris看了之后跟我讨论道:“R语言开篇讲述的基本数据类型为什么和C语言等编程语言不同?”后来得出结论:每种语言都有其最适用的领域,以R为例,它主要运用在统计学领域,处理大量数据,基本单位就是向量(vector),故将其他语言中的int类型视为最简单的一维向量,如c(5)。我后来验证了从矩阵中取出一个值,类型为vector。同样的道理,Perl语言的基本数据类型也和其他语言不同,因为它用处不同。
学习新东西的时候,不能盲目学,要勤思考,发现不同,查阅资料,融会贯通,这才是正确的学习方法。感谢Chris的指点!

Overview

R语言是一种跨平台的开源编程语言,主要用于统计分析、绘图、数据挖掘、机器学习领域。安装教程可以参考Chris写的这篇文章:Ubuntu 12.04下R的安装。由于最近开始新项目,需要用到这方面的知识,所以学习之余将这些知识点记录下来,以便不时查阅。

1. R语言概览

R语言主要应用在三个应用场景下(按照处理数据的顺序排列):数据查询(Data Exploration)、数据挖掘(Data Mining)、数据展现(Data Presentations)。在学习的过程中发现该语言的许多用法和以前自己用过的另一门统计语言matlab有很多相近的地方。

2. R语言用法

学过多门编程语言的人一定会发现,基本上编程语言的内置函数名,都是英文单词或者其缩写,当然R语言也不例外,所以R语言的语法还是比较简单的。我们从创建向量(矩阵)开始。

2.1 创建向量和矩阵

#create v1,v2
#c()表示将数值combine,生成一个向量或列表
v1=c(1,2,4,6,2) #1 2 4 6 2
v2=c(2,1,3,5,7) #2 1 3 5 7
#row bind to a matrix
m1=rbind(v1,v2) 
#   [,1] [,2] [,3] [,4] [,5]
#v1  1    2    4    6    2
#v2  2    1    3    5    7
#column bind to a matrix
m2=cbind(v1,v2) #这一步和rbind()就是行列互换了一下,就不展示了。

v1v2的数据类型感兴趣的话,可以用mode()函数查看一下

mode(v1) #结果是数字类型numeric

我们也可以创建字符类型的向量,只需将各元素加上双引号""。当然,如果是数字和字符混合类型,则统一处理为字符类型(character)
我们也可以用m=c(1:20)这样的用法创建向量,其中m20个元素,从120,也就是说,默认的步长step为1,当然,这个是可以修改的,具体用法如下。

c(1:10-1) #将从1到10这个向量每一项都-1
c(1:10*2) #将从1到10这个向量每一项都×2

还可以进行组合:

c(1:10*2+1) #四则运算有优先级,如果不确定,最好加上括号

matrix()这个函数就可以直接创建矩阵了:

matrix(1:12,nrow=3,ncol=4)#默认按列创建矩阵

默认按列创建矩阵,如果把第四个参数byrow加上,效果就会不同了

matrix(1:12,nrow=3,ncol=4,byrow=T)

其中byrowboolean类型的参数,可取TF。需要注意的是:并没有bycol这个参数。
上面的语句和下面的语句效果是一样的:

matrix(1:12,3,4,T)

也就是说,参数名不用写。
有了矩阵,就可以进行矩阵之间的运算了,当然需要满足矩阵运算的条件,具体条件可以参考线性代数,这里就不细说了。
矩阵加减:

m1+m2 #需要两矩阵行数相等,列数相等

矩阵转置:

t(m1)

矩阵相乘:

a %*% b #须a的列数与b的行数相等

对于方针,有下面的用法:

a=matrix(1:16,nrow=4,ncol=4)
diag(a) #取对角线元素1    6    11    16
diag(diag(a)) #这个用法很有意思:以1    6    11    16为对角线,其余位置补0
diag(4) #创建单位矩阵,维度为4×4

矩阵求逆:

solve(a) #a须为方阵

上式其实是solve(a,E)的简写,其中E为单位矩阵。
solve(a,b)的普遍用法就是解方程组。
求特征值和特征向量:

eigen(a,symetric)

其中a为方阵,symmetricboolean类型,可以不写。我建议不写,因为系统会自动去判断a是否为对称阵。

我们可以查看向量的长度:length(v1).那如果对m×n的矩阵用length呢?答案是m×n
还有一点:R语言是没有size函数的。
还有一个类似的函数:

seq(1,10) #产生序列1    2    3    4    5    6    7    8    9    10

默认这个等差数列公差为1,自行设置公差方法如下:

seq(1,10,by=2) #1    3    5    7    9

上面的写法中,对于我的R语言版本R version 3.2.2可以省略by=2,其他版本的可以尝试一下。
可以指定数列长度:

seq(1,10,length=11) #1.0  1.9  2.8  3.7  4.6  5.5  6.4  7.3  8.2  9.1 10.0

需要注意,公差和长度只能存在一个,不能同时使用。

2.2 从向量或矩阵中取值

我们要是想取出来向量或者矩阵当中某个元素只需v1[3]或者m1[2,3],若矩阵取值时,方括号中仅有1个数字,则按照首行首列的位置,沿列方向,若到达列尾,则继续沿着下一列的首行位置查询取值。需要注意的是,R语言取值的下标和其他编程语言从开始不同,它是从1开始的。
如果方括号中为负数,例如:

v1[-2] #将v1[2]从向量中剔除,剩余的部分将是结果
# 1    4    6    2

也可以这样连续取值:

v1[1:3] #1    2    4 

非连续取值:

v1[c(1,3,5)] #1    4    2

还有一种常见的错误写法:

v1[1,3,5] #提示“量度数目不对”

按条件取值:

v1[v1>1 & v1<5] #取出大于1且小于5的数字

以上,如果语法正确但是没有符合条件的数据,则返回:

numeric(0) #数目为0

查询下标值的函数:

which()

用法如下:

v=c(2,5,1,4,7,2,5,6,8)
which.max(v)
which.min(v)
which(v>3)
which(v==2)

以上写法返回的都是下标。
还可以排序以及逆序:

sort() #默认从小到大排列
rev() #逆序排列

2.3 求统计数据

x=c(1:100)

mean(x) #平均值mean
sum(x) #和sum
max(x) #最大数maximum
min(x) #最小数minimum
var(x) #方差variance
prod(x) #积product
sd(x) #标准差standard deviation

配置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. 参考资料