PHYLIP

1. PHYLIP简介

PHYLIP,the Phylogeny Inference Package, is a package of programs for inferring phylogenies (evolutionary trees). It can infer phylogenies by parsimony, compatibility, distance matrix methods, and likelihood. It can also compute consensus trees, compute distances between trees, draw trees, resample data sets by bootstrapping or jackknifing, edit trees, and compute distance matrices. It can handle data that are nucleotide sequences, protein sequences, gene frequencies, restriction sites, restriction fragments, distances, discrete characters, and continuous characters.

2. PHYLIP的下载和安装

下载PHYLIP,然后解压缩,得到phylip的安装文件包,其中包含doc、exe和src三个文件夹以及一个phylip.html文件。html文件结合doc中的文件在浏览器中打开是phylip详细的说明文件。当然,官网也提供了phylip 3.6的PDF版的Manual文件

$ wget http://evolution.gs.washington.edu/phylip/download/phylip-3.695.tar.gz
$ tar zxf phylip-3.695.tar.gz -C /opt/biosoft/
$ cd /opt/biosoft/phylip-3.695/src/
$ cp Makefile.unx Makefile
$ make install

OrthoMCL介绍

1. OrthoMCL的用途

基于序列的相似性,OrthoMCL能将一组proteins(比如全基因组的proteins)归类到ortholog groups、in-paralogs groups和co-orthologs。

2. OrthoMCL-DB

OrthoMCL-DB包含了很多proteins,这些proteins来自一些已经完全测序的真核或原核生物的基因组。OrthoMCL-DB将这些proteins进行了聚类,分成很多的ortholog groups。
2010.5.31,发布了OrthoMCL-DB第4版,包含 116,536个ortholog groups、1,270,853个proteins、88个真核生物基因组、16个古菌基因组、34个细菌基因组。
2011.5.31,发布了OrthoMCL-DB第5版,包含 124,740个ortholog groups、1,398,546 个proteins、150个基因组
2013年末即将发布OrthoMCL-DB第6版。

3. OrthoMCL的两种使用方法

1. OrthoMCL-DB的官网已经将数据中的proteins进行了ortholog的聚类,其网站提供了一个工具,用于接收上传的基因组proteins,再将这些proteins group到相应的ortholog groups中。官网提供的工具Assign your proteins to OrthoMCL Groups用于进行分析。
2. 如果要对多个基因组的proteomes进行聚类,则可以使用OrthoMCL单机版的软件来进行运算。其用法详见:OrthoMCL的使用

4. OrthoMCL算法

1. 将多个proteomes转换成orthomcl兼容的FASTA文件。
2. 移除低质量的序列。
3. All-versus-All BLASTP with 1e-5 cutoff。即使用这些proteomes的protein sequences构建blast数据库,再将所有的这些序列和数据库进行BLASTP比对,取evalue小于1e-5的比对结果。
4. Filter by percent match length。计算比对结果的percent match length ( 所有hsp中比对上序列的长度之和 / 两条序列中短的那条序列的长度 )。取50%的cutoff值。
5. 寻找不同物种间potential ortholog pairs(两两物种的protein序列相互是best hits);寻找同一物种内in-paralog pairs(相互之间是better hits,即对于2个序列之中的任意一条序列,和其in-paralog序列之间的evalue值 <= 这条序列和其它物种比对的evalue值). 6. 根据上一步结果寻找co-ortholog pairs(pairs connected by orhthology and in-paralog,并且pairs之间的evalue值低于1e-5). 7. 对所有的pairs进行E-values的Normalization,以利于下一步MCL的计算。见下一部分内容,或参考OrthoMCL Algorithm Document
8. 将所有的ortholog,in-paralog和co-ortholog pairs,以及它们的标准化后的weight值输入到MCL程序中,来进行聚类分群。MCL documentation

5. pairs的evalue计算和标准化

pairs的evalue计算:pairs的两条序列相互blast后有两个evalue值,这两个值常常不相等。但是为了计算需要,于是pairs的之间的两个evalue值要进行一个计算,得到pairs weight,weight= ( -log10(evalue1) + -log10(evalue2) ) / 2 。
pairs的evalue的标准化:1. 对于in-paralog pairs,在某一个基因组中,取两条序列中任意一条序列有ortholog的in-paralog pairs为有效in-paralog pairs。若在这个基因组没有这样的pairs,则该基因组所有的in-paralog pairs都为有效的in-paalog pairs。最后得到所有基因组所有有效的in-paralog pairs。然后取这些有效in-paralog pairs的weight的平均值。最后,每个in-paralog pair的evalue标准化后的值为其weight除以average weight。 2. 对于ortholog或co-ortholog pairs则简单很多,求所有weight的平均值,然后使用各个pair的weight除以average weight,则将其标准化了。

6. 网络版的OrthoMCL的使用

OrthoMCL-DB已经对150个proteomes进行了OrthoMCL的分析,对orthologs进行了聚类。这个过程由于数据量大,因此,在好几百的CPU资源下也需要好几个星期才能做完。
在OrthoMCL-DB上上载 a set of proteins,服务器则会将所有上传的proteins比对到OrthoMCL-DB中所有的proteins上;选取evalue 1e-5和50% match的cutoff;然后将protein归类到其top hit所对应的protein的类上;如果top hit所对应的proteins没有group,则该protein归类到NO_GROUP。
然后,再对上一步cutoff掉的proteins来使用OrthoMCL-DB的in-paralog算法来创建in-paralogs pairs,然后再进行MCL的聚类。
使用该方法,最后将a set of proteins进行了同源基因的聚类,但是缺点如下:
1. 这种方法是单向最佳,根据protein比对的最佳结果去归类到已有的group中去。但是反过来,最佳比对结果对应的protein不一定和query protein是最佳的。这和OrthoMCL的算法是有出入的,所以该方法省了时间,但是结果和真正的结果是有一定差别的。
2. 只使用cutoff后剩下的proteins进行in-paralog分析,而没有进行所有query proteins之间的in-paralog分析。
3. 没有ortholog pairs和co-ortholog pairs的信息,没法进行单拷贝同源基因的提取与分析。

7. 本地OrthoMCL的使用

对指定的a set of proteomes进行同源基因分析,则使用本地的OrthoMCL进行分析。而官网不提供这种服务,因为消耗的计算机资源过大。

8. 注意事项

1. 序列都是使用protein序列,而不是nucleotide序列,是因为protein序列更精确。
2. proteomes中的序列要去除可变剪切,只留取alternative proteins中长度最长的。否则在有alternative proteins存在的情况下,则会造成pseudo-in-paralogs(即alternative proteins称为in-paralogs),给后续的分析造成麻烦。
3. paralog分为in-paralog和out-paralog。in-paralog是指同一个物种的paralogs的分化发生在物种分化之后,这样的话,代表in-paralogs之间的序列相似性比其orthologs的相似性要高;通过OrthoMCL的原理可以看出,很好得到分析。而out-paralog是指paralogs的分化发生在物种分化以前,这代表out-paralogs的序列之间的相似性比某一个物种的orthologs的相似性要低;这样是很不好分析的,因为不好定阈值,或者得到的结果不易得不到大众的认可;OrthoMCL也没进行out-paralog的分析;当然,也可以只将一个query proteome输入到orthomcl来进行分析,得到的是所有paralog分析结果,包含了out-paralog。

OrthoMCL的使用

分13步进行,如下:

1. 安装和配置数据库

Orthomcl可以使用Oracle和Mysql数据库,而在这里只介绍使用Mysql数据库。
修改配置文件/etc/my.cnf,对Mysql进行如下配置:

1. 设置myisam_sort_buffer_size的值为可用内存的一半;
2. 设置myisam_max_sort_file_size为orthomclBlastParser程序生成文件similarSequences.txt的5倍大小;
3. 软件的说明文档中设置read_buffer_size的值为???,但是设置为3个问号或1个问号,则mysql启动不了。我将其设置为2000M。

2. 安装mcl软件

mcl,即Markov Clustering algorithm,其最新的软件下载地址:http://www.micans.org/mcl/src/mcl-latest.tar.gz。下载后使用’./configure && make && make install’安装即可。

3. 安装并配置OrthoMCL软件

下载OrthoMCL软件(http://orthomcl.org/common/downloads/software/)后,解压缩后,其中包含文件夹:bin、config、doc、lib四个文件夹。
在一个工作目录中运行OrthoMCL,该目录包含数据文件和结果文件。将doc/OrthoMCLEngine/Main/orthomcl.config.template复制到工作目录中。该文件为OrthoMCL的配置文件,以使用mysql数据库为例,其中的内容如下:

dbVendor=mysql   #使用的数据库为mysql
dbConnectString=dbi:mysql:orthomcl   #使用mysql数据库中名为orthomcl的数据库
dbLogin=test    #数据库的用户名
dbPassword=123  #相应的密码
similarSequencesTable=SimilarSequences #
orthologTable=Ortholog
inParalogTable=InParalog
coOrthologTable=CoOrtholog
interTaxonMatchView=InterTaxonMatch
percentMatchCutoff=50
evalueExponentCutoff=-5
oracleIndexTblSpc=NONE

4. 安装orthomcl数据库的表

首先,进入mysql数据库,新建一个名为orthomcl的数据库;然后,使用orthomclInstallSchema命令在数据库中创建一些表,这些表的名字则是orthomcl.config.template配置文件中指定的5个名称。

$ mysql -u test -p
mysql> create database orthomcl;
$ cp /opt/biosoft/orthomclSoftware-v2.0.9/doc/OrthoMCLEngine/Main/orthomcl.config.template .
$ orthomclInstallSchema orthomcl.config.template [log species]

orthomclInstallSchema命令后面不接参数则给出帮助文档。以上命令使用orthomcl.config.template配置文件中的设置生成了数据库中相应的表,如果加入方括号中的内容,则会生成日志文件log,生成的表后缀都是species。

5. 创建orthomcl的输入文件

orthomcl的输入文件为fasta格式文件,但是fasta文件的序列名称要满足这样的要求:

>taxoncode|unique_protein_id
MHDR...
>hsa|sequence_1
MHDR...
>led|scaffold_1.1
MHDR...

序列名第一列是物种的代码,一般是3到4个字母;中间使用’|’符号隔开;第二列是蛋白质序列独一无二的id。
一般输入文件是fasta格式,其序列名由空格或’|’隔开,使用orthomclAdjustFasta程序,将fasta文件转换出兼容orthomcl的fasta文件

$ redun_remove protein.fasta > non_dun_protein.fasta
$ mkdir compliantFasta; cd compliantFasta
$ orthomclAdjustFasta led ../non_dun_protein.fasta 1

上述命令去除可变剪切的蛋白质序列;创建了文件夹compliantFasta;然后使用orthomclAdjustFasta命令选取了protein.fasta序列名的第一列作为输出的fasta文件的序列id;输出的文件为led.fasta.

6. 过滤序列

对compliantFasta文件夹中的序列进行过滤,允许的最短的protein长度是10,stop codons最大比例为20%;生成了两个文件goodProteins.fasta和poorProteins.fasta两个文件。

$ orthomclFilterFasta compliantFasta/ 10 20

compliantFasta只能过滤低质量的序列。而实际上最好还需要过滤掉可变剪切,只留取可变剪切中最长的蛋白质序列,这个需要自行解决。

7. 对goodProteins.fasta中的序列进行BLAST

下载最新版本的Blast+,和最新版本的OrthoMCL DB的protein序列,将OrthoMCL DB的protein序列加上gooProtein.fasta中的序列合到一起做成一个blast+的数据库。然后对基因组的蛋白质序列进行比对。

$ /opt/biosoft/ncbi-blast-2.2.28+/bin/makeblastdb -in orthomcl.fasta -dbtype prot -title orthomcl -parse_seqids -out orthomcl -logfile orthomcl.log
$ /opt/biosoft/ncbi-blast-2.2.28+/bin/blastp -db orthomcl -query goodProteins.fasta -seg yes -out orthomcl.blastout -evalue 1e-5 -outfmt 7 -num_threads 24

生成orthomcl的blast DB需要97秒左右;使用-outfmt 7生成带注释的表格结果,这一步需要很长时间了,取决于电脑的运算性能。我使用24个线程,每分钟运行约27.75条序列,大约7.2个小时,运行1.2万条protein序列的比对。
blast中使用了-seg yes表示使用seg程序来进行过滤,将那些影响比对结果的低复杂度区域过滤掉。blast生成的文件结果,从第1列到第12列分别是:query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, q. end, evalue, bit score。

8. 处理Blast的结果

对上一步blast的结果进行处理,从而得到序列的相似性结果,以用于导入到orthomcl数据库中。compliantFasta文件夹中包含下载下来的OrthoMCL DB的所有蛋白质数据的文件orthomcl.fasta.

$ grep -P "^[^#]" orthomcl.blastout > blastresult
$ orthomclBlastParser blastresult compliantFasta > similarSequences.txt
$ perl -p -i -e 's/\t(\w+)(\|.*)orthomcl/\t$1$2$1/' similarSequences.txt
$ perl -p -i -e 's/0\t0/1\t-181/' similarSequences.txt

第一条命令将orthomcl.blastout中的注释行去掉,生成新的文件blastresult,不然再下一个命令中会报错的。
第二条命令生成文件similarSequences.txt,从第1列到第8列分别是:query_id, subject_id, query_taxon, subject_taxon, evalue_mant, evalue_exp, percent_ident, percent_match。值得注意的是subject_taxon是orthomclBlastParser读取的在compliantFasta文件夹中fasta文件的前缀,在此结果中,这一列则全是orthomcl。
第三条命令将subject_taxon修改为正确的分类名。
第四条命令修改evalue_mant, evalue_exp,将evalue为0修改为1e-181,这在后续步骤寻找pairwise relationships时候有要求。

9. 将similarSequences.txt载入到数据库中

生成的similarSequences.txt文件大小为83M,则修改/etc/my.cnf文件,在myisam_sort_buffer_size这一行上加一行‘myisam_max_sort_file_size = 424960’。数值是83M的5倍。然后运行:

$ orthomclLoadBlast orthomcl.config.template similarSequences.txt

10. 寻找成对的蛋白质

$ orthomclPairs orthomcl.config.template orthomcl_pairs.log cleanup=no

输入为数据库中的表SimilarSequences,和数据库的空表InParalog, Ortholog, CoOrtholog tables;输出为对这些空表的操作。故配置文件中的用户要有 update/insert/truncate权限。

11. 将数据从数据库中导出

$ orthomclDumpPairsFiles orthomcl.config.template

生成了一个文件mclInput和一个文件夹pairs;文件夹中包含3个文件coorthologs.txt,inparalogs.txt,orthologs.txt。

12. 使用mcl进行对pairs进行聚类

$ mcl mclInput --abc -I 1.5 -o mclOutput

13. 对mcl的聚类结果进行编号

$ orthomclMclToGroups led 1 < mclOutput > groups.txt

对聚类结果进行编号,依次为led1,led2, led3…

注意事项

第7步Blast,是整个过程中最关键的一步。有以下2点需要注意:
1. 数据库中的蛋白质序列数量:在OrthoMCL DB中选取和要分析的物种亲缘关系较近的几个物种的基因组,或下载其它公布的基因组,加上要分析的物种的基因组;使用这些基因组总体的蛋白质序列来构建Blast数据库。如果只是使用要分析的物种的蛋白质序列建数据库,则inparalogs文件中成对的序列实际上是paralogs,数目比真正的inparalogs要多很多。使用所有的OrthoMCL DB中的序列,第5版版含150个基因组,信息量太大,不使用几百个核的超算或计算机集群去运行,是很不现实的。
2. 对数据库中所有的蛋白质序列来使用blast比对到该数据库中得到结果。如果只是对要分析的物种进行Blast,则只能得到inparalogs的信息,而没有orthologs和coorthologs。

Linux下PHP的mcrypt模块安装

参考:http://www.linuxidc.com/Linux/2009-05/20108.htm

1. 下载并安装libmcrypt-2.5.8.tar.gzmhash-0.9.9.tar.gzmcrypt-2.6.8.tar.gz

2. 安装php-devel,只有安装了php-devel后才会有命令phpize。

3. 下载相应的php版本的源码包,并解压缩后进入源码包的ext/mcrypt目录,执行:

# phpize
# ./configure && make && make install
# echo 'extension=extname.so' >> /etc/php.ini

CentOS安装php5.2

捣鼓了好久才搞定。

当前,CentOS6默认安装的php是5.3版本,使用默认源更新着会将php更新到5.5版本;CentOS5默认安装的php是5.1.6版本,使用默认源则不会更新。而当要使用最新的wordpress时候,需要的php最低版本是5.2.4;而要使用很多电子商务网站的php源码的时候,需要的php版本则要低于5.3版本,同时需要安装ZendOptimizer。

1.使用本地源来安装mysql和php及其组件

为此,于是重新安装额CentOS5.9的系统,然后通过系统盘创建了本地源,并使用本地源更新了php和mysql的其它模块。

# alias yumlocal='yum --disablerepo=\* --enablerepo=c5-media'
# yum install php* *mysql*  --skip-broken

只有更新php和mysql的其它模块,电子商务等php源码才能正常挂到网上。

2.更新php到5.2版本

参考该博客的一篇文章:http://www.webtatic.com/packages/php53/,提供了一个源用于更新php到5.2。

# rpm -Uvh http://mirror.webtatic.com/yum/centos/5/latest.rpm
# echo 'xclude=php*5.3*' >> /etc/yum.conf
# yum --enablerepo=webtatic install php
# yum --enablerepo=webtatic update php
# yum --enablerepo=webtatic install php php-xml php-ldap php-devel php-common php-gd php-mysql php-cli php-mbstring php-dba php-bcmath php-pear

以上第1条命令添加源;第2条修改yum配置文件,设置php不要更新到5.3;第3条安装php;第4条更新php及相关模块。

3.安装ZendOptimizer

现在的很多网站需要ZendOptimizer来进行网站加速加密。但是ZendOptimizer仅运用php5.3版本以下。当然ZendOptimizer有相应的升级版本Zend Guard Loader,支持php5.3及以上版本。但是很多网站的运行在Zend Guard Loader和php5.3及以上版本中无法运行。因此,需要安装ZendOptimizer。

使用搜索引擎搜索“ZendOptimizer3.3.9 linux x86_64”来下载ZendOptimizer3.3.9版本,然后进行安装:

# tar zxf ZendOptimizer-3.3.9-linux-glibc23-x86_64.tar.gz
# cp ZendOptimizer-3.3.9-linux-glibc23-x86_64/data/5_2_x_comp/ZendOptimizer.so /usr/lib64
# echo 'zend_extension=/usr/lib64/ZendOptimizer.so' >> /etc/php.ini
# /etc/init.d/httpd restart

我的第一届生物信息学培训班总结

1 培训的准备

在今年4月份的时候就萌生了一个想法,想将自己学到的一些生物信息学知识分享给它人,同时也赚点生活费。后来偶然在一个生物信息的群里面,有人说很想在暑期去上生物信息学班,学习到实用的知识;同时又抱怨华大等举办的生物信息学培训班不仅学费贵,同时也学不到自己动手的东西。于是我就写到考虑在暑期办一个培训班。于是当时一冲动,就在华中农业大学的南湖论坛上发了一个帖子,说想要开一个生物信息学的班。后来很多人回复并感兴趣,于是这成了我举办生物信息学培训班的最大动力。后来在一些回复的要求下建立了一个qq群,于是将培训班的事业拉倒正轨上。

开班培训班,则必须要有培训资料,于是就开始整理培训资料。其实这些培训资料很多都是取自我零零散散写到我的blog上的内容。在今年暑假7.14-7.24日,我参加了复旦大学举办的一个进化分析的培训班。那个时候白天上课,晚上就在宾馆写培训资料,同时准备上课的数据。培训资料和数据的准备绝大部分是在那个时候完成的,当然,由于比较匆忙,在最后一章的内容就写得粗糙了,同时培训资料中还有很多的错别字等。

通过qq群的通知,培训时间确定在7.27日下午3点进行一个试讲,讲述培训的大纲;而从7.28日开始正式的讲课。

2 培训内容

这次的生物信息培训主要是NGS数据的分析,是现在比较成熟又热门的知识。讲解了基因组的组装、基因预测、SNP分析、转录组分析、基因功能注释等,还包含了网站搭建,mysql使用和perl编程等。全程课程主要偏重实践运用于上机操作,而原理方面则不会深究。

本次培训的大纲如下:

1.CentOS 6.4 x86_64的安装
2.Linux系统入门
3.Next Generation Sequencing Techonology
4.基因组组装
5.Genom Repeat Sequence Prediction
6.短序列比对
7.Perl入门
8.Variants calling
9.转录组分析
10.基因预测
11.基因组浏览器Gbrowse
12.基因功能注释与富集分析

培训讲义资料有151面,内容还是很丰富的。

3 培训情况

在7.27日下午3点开始培训前试讲,当时天气很热,空调刚开,但是大出我的预料,来了很多人,将会议室都挤满了。现场报名,报名费相对便宜,500元。由于场地不够大,于是随机抽取了部分有意向的人报名,最后总共报名人数27人。当然,有一个高级博后从上海专程到这儿来听课的,像这样的需要优先考虑。报名时,当时收钱的感觉还是很爽的,嘿嘿…这算是我人生的第一桶金,但是紧接而来则是要多学员负责,要教给他们有用的知识了。

培训,准确说是集训,从7.28日开始,每天早上9:00-12:00,下午2:00-5:00,晚上7:00-9:00上课。这样总共持续了10天,其中包含放假了一天。刚开始两天,早中晚都讲课,后来晚上就专门用来练习,手把手答疑了。最后一天,人人都很累,同时培训资料最后一章也写得比较简略,数据也准备的不是很充分,于是结束课程也算比较匆忙。想起来,当初我做刚生物信息学,就是从最后一章基因功能注释的内容开始的,对这一章的内容,其实了解得非常深入。

最后,所有的学员在实验室门前合了一张影,作为留恋:
陈连福的第一届培训班学员合影1
陈连福的第一届培训班学员合影2

4 我的培训收获

通过这次培训,我认识了一些人,很多人给我留下了深刻的印象。比如:从上海特异过来的高级博后,本身对NGS分析很感兴趣,学得很快;基因楼的爱学习,喜欢辩论的小伙子;在测序公司工作几年后回来继续读书的,问题很多的同学;还有刚发表了Nature techonology马上去瑞士读博士的牛人;在空调口吹风又很活跃的一群人;也有喜欢钻研的女生等。

培训课上也有一些同学退学,也有一些同学加入,也有人苦苦支撑到最后。

作为老师,我也学到了很多上课经验。在培训课上,上课是讲一会儿,上机一会儿,造成了课间休息不规律;上课节奏要慢,吐词要清晰等。当然,对培训课程进度的把握,现在有了经验。原本准备在3~5天内讲完,但是在很快的节奏下也需要了9天的时间。

5 总结

本次培训,通过手把手教导怎么进行NGS分析与上级操作,学员们最后能独自自主进行高质量的:基因组组装,SNP和INDEL检测,基因预测及其可视化,转录组分析,基因功能注释等。

培训资料最后只是发放了纸质版的,电子版的由于版权等问题暂时保留在我手上。以后要是有时间,我考虑要出一本NGS分析的生物信息学书籍。

培训班在8.6日结束,随后8.13-8.25参加了在哈尔滨医科大学举办的生物信息学培训。于是原本预算的第二届生物信息学培训班就没有举办了。这个时候群里面的王志文在群里面打了下广告,貌似在洪山体育馆那儿举办了一期生物信息学培训班。

慢慢筹备我的下一次的培训吧,估计内容会更加丰富…

CentOS安装zend

PHP 5.3 下,Zend Optimizer 已经被全新的 Zend Guard Loader 取代。

1. 下载对应版本的 Zend Guard Loader 压缩包,使用”php –version”参看php版本,然后下载对应php版本和系统32位或64位版本的安装包。(http://www.zend.com/en/products/guard/downloads)

2. 解压缩得到ZendGuardLoader.so(Linux)或 ZendLoader.dll(Windows)。将之复制到指定的系统路径下,比如 /usr/lib64/

3. 修改 /etc/php.ini 文件,在首行以下加入如下几行:

; Add the following line to your php.ini file for loading the ZendGuardLoader:
zend_extension=/usr/lib64/ZendGuardLoader.so
; Add an aditional line to your php.ini for enabling ZendGuardLoader
; Enables loading encoded scripts. The default value is On
zend_loader.enable=1
; Disable license checks (for performance reasons)
zend_loader.disable_licensing=0
; The Obfuscation level supported by Zend Guard Loader. The levels are detailed in the official Zend Guard Documentation. 0 - no obfuscation is enabled
zend_loader.obfuscation_level_support=3

4. 重启web服务,使用“php -v”查看安装情况。

ssh使用80端口进行登录

校园网只对外开放80端口,当多天出差的时候需要登录服务器。此时,需要将ssh默认的22端口改为80端口。同时,若占用了80端口,服务器对外的www服务则不能访问了,为了解决此矛盾,按如下方式解决:

1. 修改端口

# cd; pwd
/root
# cat > chport.sh
#!/bin/bash
perl -p -i -e 's/^Port 22/Port 80/' /etc/ssh/sshd_config
perl -p -i -e 's/^Listen 80/Listen 8080/' /etc/httpd/conf/httpd.conf

/etc/init.d/sshd stop
/etc/init.d/httpd stop

/etc/init.d/sshd start
/etc/init.d/httpd start

clear
ctrl + d
# chmod 755 chport.sh
# /root/chport.sh

以上命令建立了一个脚本文件/root/chport.sh。该文件用于修改sshd和http服务的端口,将sshd端口改为80,将httpd的端口改为8080.

这样,修改端口后,即可在外网使用如下命令登录服务器:

$ ssh -X -p 80 user@ip

2. 连接服务器后再改回端口

为了正常进行www服务器的访问,需要在登录进服务器后,将端口改回。这样通过ssh连接到了服务器,同时也能进行www访问了。

# cd; pwd
/root
# cat > rtport.sh
#!/bin/bash
perl -p -i -e 's/^Port 80/Port 22/' /etc/ssh/sshd_config
perl -p -i -e 's/^Listen 8080/Listen 80/' /etc/httpd/conf/httpd.conf

/etc/init.d/sshd stop
/etc/init.d/httpd stop
ctrl + d
# chmod 755 rtport.sh
# /root/rtport.sh

3. 登录到root时自动还原端口,退出root时自动修改端口

登录到了root,则不需要修改后的端口了,这是需要正常的www服务了。这时可以修改.bashrc文件;而退出root的时候,则需要修改端口,否则,下次登录就登录不上去了。这时,可以修改.bash_logout文件。

# cat >> /root/.bashrc
/root/rtport.sh
ctrl + d
# cat >> /root/.bash_logout
/root/chport.sh
ctrl + d

4. crontab定时修改或还原端口

若在还原端口后,突然断网导致中断链接,则再也不能登录服务器了,那不就悲剧了,特别是出差在外,不能回校的时候。此时就需要使用crontab来定时运行上述脚本:

# crontab -e
0       *       *       *       *       /root/chport.sh
30      *       *       *       *       /root/rtport.sh
:wq

上述命令在每小时0分钟的时候修改端口,于是能使用ssh进行远程登录了;每小时30分钟时候,则还原端口,可以正常访问www服务。

当然,在使用ssh登录后,在sudo到root权限,则自动还原端口,正常访问www服务。必要的时候手动运行以上脚本来开启www服务以维持www服务的运行。

5. 解决httpd不能重启的问题

这样反复重启服务,可能造成httpd重启的时候出错,提示程序已死但是pidfile存在的错误。则需运行如下命令来解决:

# ipcs -s | grep apache | \
perl -e 'while(<>){@a=split /\s+/, $_;print "ipcrm sem $a[1]\n"}' | \
sh
# rm -f /var/lock/subsys/httpd
# /etc/init.d/httpd start

神秘学:轮回,灵魂,通灵…

1. 灵魂与轮回

今天一QQ好友说信主已有1年多,出于好奇,我于是对其信仰和他讨论了一下。

为什么信主,他说:“耶稣是永生神的儿子,他为了世人的罪被订在十字架上死了;人是有灵魂的,信了耶稣,死后灵魂就可以去天堂和神在一起。”

他说其小时候的经历让他相信了灵魂的存在:“以前小时候我大伯带我去游泳,我想从游泳池边上跳到水中游泳圈里,结果跳歪了,没进游泳圈旁边去了。然后我就看见我在水里挣扎,看见哦,就是灵魂出窍了。”

我想估计是他小时候的事情让他相信了灵魂的存在,现在他在国外接触了一些教义,顺其自然地信主。他推荐我读这本书:Many Lives Many Masters。他说通过这本书能知道灵魂确实是存在的。

《Many Lives Many Masters》这本书的中文名是《前世今生,生命轮回的前世疗法》,作者Brian L. Weiss。可以在百度百科中查看本书的介绍

本书中描述了人在死后,灵魂进入一道温暖的白光中,可能会与灵性大师交流,然后进行下一轮命运。作者作为心理医生,对一名女性患者使用催眠术进行治疗,揭露了其很多轮前世的记忆片段;而作者重点关注了其死亡时的片段。

2. 通灵

在《Many Lives Many Masters》书最后,作者说其女性患者又通过一名通灵占星师对其前世的预测,和催眠得到的结果一致。

于是我有又想找一些通灵相关的东西,于是找到了介绍嘎玛仁波切传奇事件的博客:湖心亭看雪客新浪博客和有关轮回的《西藏生死书》

基因组Repeat Sequence的寻找

1. 重复序列简介和相关软件

参考自文献:Review:Identifying repeats and transposable elements in sequenced genomes: how to find your way through the dense forest of programs。E Lerat。Heredity (2010) 104, 520–533; doi:10.1038/hdy.2009.165; published online 25 November 2009

1.1 Repeats的分类

基因组中的repeats依据其序列特征分成2类:串联重复(tandem repeats) 和 散在分布在基因组中的重复序列(interspersed repeats).其中第二类主要是transposable elements(TEs).

第一类串联重复包含:microsatellites 或 simple sequence repeats(1-6个碱基为一个重复单元) 和 minisatellites(10-60个碱基的长序列为一个重复单元).

TEs包含2种类型:class-I TEs通过RNA介导的(copy and paste)机制进行转座;class-II TEs通过DNA介导的(cut and paste)机制来转座. 前者称为retroelements,后者称为DNA transposons。

class-I TEs中主要由LTR(long terminal repeat)构成。LTR的部分序列可能具有编码功能。而non-LTR则包含2个子类:LINEs(long interspersed nuclear elements)和SINEs(short interspersed elements),其中前者可能具有编码功能,后者则没有。

class-II TEs中加入了一个子类 MITEs(miniature inverted repeat transposable elements),基于DNA的转座因子,但是确通过”copy and paste”的机制来转座(Wicker et al., 2007)。

1.2 鉴定重复序列的软件

对于不同的重复序列,需要使用不同的软件来进行鉴定。而鉴定的方法可以分为:基于library,基于重复序列的特定结构 或 重头预测。文献中给出的软件很多:http://www.nature.com/hdy/journal/v104/n6/fig_tab/hdy2009165t1.html#figure-title

1.2.1 基于library-based的软件

library-based的软件,需要构建library,该library中包含很多来自不同物种某一重复序列的一致性序列,然后通过相似性比对来鉴定repeats。这种方法能对所有的种类的重复序列进行鉴定。此方法最经典最流行的软件是RepeatMasker;此方法中CENSOR和MASKERAID两个软件可以用于改良RepeatMasker的结果;此外,用于基因组的重复序列鉴定的还有GREEDIER(Li et al.,2008),该软件在其文章中表明该软件性能还不错,在repeats鉴定的敏感性上稍微比RepeatMasker高一点,但是repeats的鉴定率只有RepeatMasker的一半左右.

1.2.2 基于signature的软件

基于signature的方法主要用于TEs的鉴定。
1. 鉴定LTR逆转座子: LTR_STRUC (McCarthy and McDonald, 2003), LTR_PAR (Kalyanaraman and Aluru, 2006), FIND_LTR (Rho et al., 2007), RETROTECTOR (Sperber et al., 2007), LTR_FINDER (Xu and Wang, 2007) and LTRHARVEST (Ellinghaus et al., 2008)。文献中,这些软件中LTR_STRUC的敏感性最高(96%),但是LTR的鉴定率只有30%;而LTRharvest的鉴定率最高(42%),敏感性67%.总体上,作者依次推荐的软件是LTRHARVEST和FIND_LTR(敏感性83%,鉴定率37%)。
2. 鉴定non-LTR retrotransposons: TSDFINDER (Szak et al., 2002), SINEDR (Tu et al., 2004) and RTANALYZER (Lucier et al., 2007)。其中,第一个软件用于验证RepeatMasker检测到的L1 insertions;第二个软件用于检测侧翼有TSDs(target site duplications 当重复序列插入到基因组上时,其两侧会带入短核酸序列的重复)的SINEs;第三个软件通过一些特征,比如TSDs,polyA尾和5’端核酸内切酶位点等来通过打分来检测L1逆转座子。
3. 鉴定MITEs:FINDMITE (Tu, 2001), TRANSPO (Santiago et al., 2002), MITE Analysis Kit (MAK) (Yang and Hall, 2003) and MITE Uncovering SysTem (MUST) (Chen et al., 2009)。文献中作者使用第一个软件报错,使用第三个软件却下载不到。第二个软件不能寻找新的MITEs,看来最好是使用最新的第四个软件。
4. 鉴定helitrons: HELITRONFINDER。该软件(Du et al., 2008)用来寻找玉米基因组中的helitrons(在动物和植物中有发现)。

1.2.3 重头预测的软件

1.2.3.1 自我比较的方法

通过BLAST、PALS等方法,将序列和自身进行比较,从而找出重复序列。软件有:REPEAT PATTERN TOOLKIT (Agarwal and States, 1994), RECON (Bao and Eddy, 2002), PILER (Edgar and Myers, 2005) and the BLASTER suite (used in Quesneville et al., 2005).其中RECON软件使用最广泛。

1.2.3.2 k-mer and spaced seed approaches

一定长度的k-mer出现了多次,可以鉴定为重复序列;spaced seed则是k-mer的一种衍生,seed上允许有一定的差异。

软件有:REPUTER (Kurtz and Schleiermacher, 1999), VMATCH (Kurtz, unpublished), REPEAT-MATCH (Delcher et al., 1999), MER-ENGINE (Healy et al., 2003), FORREPEATS (Lefebvre et al., 2003), REAS (Li et al., 2005), REPEATSCOUT (Price et al., 2005), RAP (Campagna et al., 2005), REPSEEK (Achaz et al., 2007), TALLYMER (Kurtz et al., 2008) and P-CLOUDS (Gu et al., 2008).

1.2.4 其它重复序列鉴定软件

其它一些用于鉴定非TEs的软件:Tandem Repeats Finder (TRF) (Benson, 1999), Tandem Repeat Occurrence Locator (TROLL) (Castelo et al., 2002), MREPS (Kolpakov et al., 2003), TRAP (Sobreira et al., 2006) and Optimized Moving Window Spectral Analysis (OMWSA) (Du et al., 2007) have been developed specifically to detect tandem repeats. The Inverted Repeat Finder (IRF) program (Warburton et al., 2004) was designed to search for inverted repeats.

1.3 多个软件整合的pipeline程序

REPEATMODELER pipeline (Smit, unpublished http://www.repeatmasker.org/RepeatModeler.html) includes the programs RECON, REPEATSCOUT, REPEATMASKER and TRF. It uses the output of the RECON and REPEATSCOUT programs to build, refine and classify consensus models of putative interspersed repeats.

当然,在此文献中,也讲述了其它很多专门用途的其它pipeline软件。而REPEATMODELER pipeline是现在运用于基因组的重复序列鉴定最主流的软件。以下将讲述该软件运用。

2. RepeatModeler的安装与使用

2.1 软件的安装

RepeatMaskerRepeatModeler是ISB(Institute for Systems Biology)的软件。ISB is located in the South Lake Union neighborhood。

根据RepeatMasker说明,其安装与使用需要:Perl 5.8.0以上版本,序列比对Engine,TRF和Repeat Database。其中序列比对engine可以安装多个,但每次只能使用其一,可以使用Cross_match,RMBlast,HMMER和ABBlaast/WUBlast等。

根据RepeatMideker说明,其安装与使用需要:Perl 5.8.8以上版本,RepeatMasker,Repeat Database,RECON,RepeatScout,TRF和序列比对engine。其中序列比对engine可以安装多个,但每次只能使用其一,可以使用RMBlast和ABBlaast/WUBlast。

再安装完毕需要的软件后,对RepeatMasker和RepeatModeler进行configure的时候填入相应软件的路径即可安装。

2.2 RepeatModeler的使用

2.2.1 使用RepeatModeler来通过基因组序列构建library

$ $RepeatModelerHome/BuildDatabase -name species \
  -engine ncbi species.genome.fasta
$ $RepeatModelerHome/RepeatModeler -database species

结果生成了一个文件夹,名称为RM_[PID].[DATE] ie. “RM_5098.MonMar141305172005″。该文件夹中的”consensi.fa.classified”即为library,用于RepeatMasker的输入。值得注意的是,可能fasta文件需要序列间不能有空(空格和换行等),否则会程序出错。

2.2.2 使用RepeatMasker来进行重复序列掩盖和重复序列计算

cp RM_5098.MonMar141305172005/consensi.fa.classified .
mkdir Repeat_result
$ $RepeatMaskerHome/RepeatMasker -pa 8 \
  -e ncbi -lib consensi.fa.classified \
  -dir Repeat_result -gff species.genome.fasta

则生成的结果文件位于Repeat_result文件夹中。