This site is used to record and share the informations & useful tips during the coding and research.
Overview1. separatePosNeg.py 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
Overview1. 按顺序合并文件2.比较两个文件的区别参考资料 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命令实现按顺序合并文件
Overview1. t34pssm.py2. pssm_ft.py Overview 昨天跟Chris讨论SVM分类预测准确性的时候，知道PSSM_AC的作用比PSSM作用更明显，于是决定将以前的python脚本改进一下，输出PSSM和PSSM_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=",")
Overview1. 安装R1.1 添加源1.2 更新源1.3 安装2. 安装 R packages参考资料 Overview 虽然已经用了很久的R语言，但一直没整理过，正好需要在我们的云服务器上安装R，所以一并记录下来了。下面的过程虽然是在Ubuntu 14.04上安装的，但是对于其他版本的系统，R和R 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。 如果你想查看一个函数的用法，比如mean函数（R中用来求均值的函数），可以使用 help(mean) 或者 ?mean .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/library是R自带包的所在目录。 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包安装的时候会讲到。 install.packages("mypackage")，用来安装的函数，默认会安装在.libPaths()列出来的第一个目录下。这个函数也提供了参数lib用来指定你要把包安装的位置。 library("mypackage")，载入一个已经安装的包，如果还没有安装这个包，则会报错，程序中断。 require("mypackage")，载入一个已经安装的包，如果还没有安装这个包，则显示警告信息，然后程序继续向下执行。（与library("mypackage")相比，require强烈不推荐）。 remove.packages("mypackage")，卸载package。 update.packages( )，更新package，可以定期执行。 安装Bioconductor 包： Bioconductor提供了一个函数用来安装Bioconductor 包： 首先，载入这个安装函数： source("http://bioconductor.org/biocLite.R") 之后，就可以用这个函数安装 Bioconductor 包了： biocLite("mypackage") 下面还有一些方法，不太常用，列出来参考： 下面的内容主要摘自：【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"中的所有对象。 参考资料 UBUNTU PACKAGES FOR R Ubuntu-14.04 安装 R 语言开发环境
Overview1. 常规的shell命令1.1 查看文件的前/后n行1.2 查看文件的行数1.3 查看文件的列数1.4 查看序列文件中的序列数1.5 文件合并1.6 单个文件去掉重复的行1.7 文件的交集，并集1.7 查看及关闭端口参考资料 Overview 一直想把常用的命令搜集起来，以便平时用到的时候查阅，可惜一直没抽出来时间专门整理下。最近在做序列的特征提取和多个特征文件合并时，频繁使用到了一些命令，干脆一并整理到这里，以后边用边添加整理新的命令。 这里的linux命令主要在MAC 10.9.5和Ubuntu 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 参考资料 做生物信息常用到的linux命令 Shell或命令行计算数组或文件的行数和列数 linux文件合并，去重，分割 mac下查找端口是否占用及关闭占用端口进程