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

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命令实现按顺序合并文件

生物信息中常用的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
    

参考资料

蛋白质序列特征提取方法之——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编程方案已解析完毕。疏漏之处请与我们联系,一定尽快修改。

Ubuntu 12.04下R的安装

Overview

很久以前,安装R以及R的程序库时,遇到了一些问题,当时做了笔记,现在整理一下。

1.安装R

直接在ubuntu 12.04上安装的R版本是2.14.2,安装ggplot2总是失败。需要在软件源里添加第三方软件源
命令如下:

sudo sh -c "echo deb     http://mirror.bjtu.edu.cn/cran/bin/linux/ubuntu precise/ >>/etc/apt/sources.list"

然后

sudo apt-get update 

可能会报一个公共密钥不存在的错,无法更新,原因在于第三方软件源未得到ubuntu服务器验证,需要一条命令解决:

暂时忘了是什么了,遇到这个错误了在解决吧~

然后重新

sudo apt-get update 

此时

apt-get install r-base

就可以装上最新版本的R软件。

ubuntu14.04上应该不用这么麻烦,ubuntu14.04上软件源比较新,直接安装就可以了。

2.升级R

R 2.x 升级 3.x 需要重新(编译)安装所有包:

  1. 运行R

  2. 执行:

    update.packages(checkBuilt = TRUE, ask = FALSE)
    

3. 安装ggplot2

  1. 运行R

  2. 执行命令:

    install.packages("ggplot2")

4.安装其他R程序库

我在运行figure-10-twenty-flow-congestion-preprocess.R 对数据进行预处理时,发现doBy包不存在,安装doBy:

install.packages("doBy")