您正在查看: 2016年2月

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

生物信息中常用的Linux命令(二)

Overview

做项目时出现过某些蛋白质序列出现O或者X等情况,导致计算出的PSSM矩阵也有问题。今天又遇到这种情况,在比对文件的时候,用到了两条文件操作的linux命令,记录一下。其他更多的内容参考之前Chris写的另一篇文章生物信息中常用的Linux命令

1. 按顺序合并文件

普通的合并文件可以直接用一个cat命令,而按顺序合并多个文件必须遍历这些文件,逐个合并。命令如下:

for ((i=1;i<=k;i++))do echo fileName$i;done | xargs -i cat {} >> newFileName

2.比较两个文件的区别

diff的不加参数的基本命令最常用。

diff file1 file2

参考资料

Shell脚本中使用for循环和cat命令实现按顺序合并文件

改进计算PSSM的python脚本

Overview

昨天跟Chris讨论SVM分类预测准确性的时候,知道PSSM_AC的作用比PSSM作用更明显,于是决定将以前的python脚本改进一下,输出PSSMPSSM_AC这两个文件,方便观察。该脚本包括两部分,本文将按顺序记录下来。
以前的脚本可以参考我之前的文章蛋白质序列特征提取方法之——PSSM

1. t34pssm.py

#! /usr/bin/env python
import fileinput
import sys
from os import listdir
from os.path import isfile, join
import re
from pssm_ft 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]
    finalist = []
    for smp in smplist:
        finalist.append(pssmdir+'/'+fastaDict[smp].split('.')[0]+'.pssm')

    for fi in finalist:
        pssm(fi, 'total_train_60.pssm','total_train_60.pssm_AC') #此处多了一个参数,多了一个输出文件

2. pssm_ft.py

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

GROUP = 10

def scale(x, max, min):
"""
scale x to [-1,1]
"""
    try:
        return -(1-(x-min)/(max-min))+((x-min)/(max-min))
    except Exception:
        print "ZeroDivisionError, max-min:\s", max-min

def pssm_ac_cal(PSSM, g, l):
"""
:param PSSM: ['1','-1'...]
:param g: group
:param l: sequence length
:return PSSM_AC: g * 20 matrix
"""
    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

def pssm(fi,output,output_AC):
    # 0-19 represents amino acid 'ARNDCQEGHILKMFPSTWYV'
    Amino_vec = "ARNDCQEGHILKMFPSTWYV"
    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]
            if len(str_vec) == 0:
                break
            PSSM.append(map(int, str_vec[1:]))
            seq_cn += 1
            #test1_1 += float(str_vec[1]) test1_1 = -103
            ACC[Amino_vec.index(str_vec[0])] = map(sum, zip(map(int, str_vec[1:]), 
                ACC[Amino_vec.index(str_vec[0])]))  
    fileinput.close()
    # section for ACC features 
    ACC_np = np.array(ACC)
    ACC_np = np.divide(ACC_np, seq_cn)

    amcnt_max = np.amax(ACC_np)
    amcnt_min = np.amin(ACC_np)

    vfunc = np.vectorize(scale)
    ACC_np = vfunc(ACC_np, amcnt_max,amcnt_min)
    ACC_shp = np.shape(ACC_np)

    #section for PSSM_AC features
    PSSM = np.array(PSSM)
    PSSM_AC = pssm_ac_cal(PSSM, GROUP, seq_cn)
    PSSM_shp = np.shape(PSSM_AC)

    file_out = file(output,'a')
    file_out_AC = file(output_AC,'a')
    #此处将保存文件的代码拆分开
    np.savetxt(file_out, [np.reshape(ACC_np, (ACC_shp[0] * ACC_shp[1], ))], delimiter=",")
    np.savetxt(file_out_AC, [np.reshape(PSSM_AC, (PSSM_shp[0] * PSSM_shp[1], ))], delimiter=",")

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
    

参考资料