Jun 05 2008

使用python进行ClustalW格式和aligned FASTA格式转换

分类: programming, studyazalea @ 5:49 pm 367次阅读

ClustalW格式是多序列比对程序CLUSTAL W的默认输出格式,大概是这个样子:

Aligned FASTA格式是多序列比对程序MUSCLE的默认输出格式,与一般的FASTA格式很相似,只是根据比对结果在序列中填充“-”,使所有序列长度相同,例如:

Aligned FASTA格式更便于机器阅读,而ClustalW格式是比较人性化的 ,使用者可以直观地看到相同和不同的区域。

由于MUSCLE程序的profile alignment功能需要aligned FASTA格式的输入,因此我写了以下一段简单的python代码,把ClustalW格式转换成aligned FASTA格式:

PYTHON:
  1. #!/usr/bin/python
  2. #Usage: python clw2afa.py input.clw
  3. #Default output is STDOUT, you can redirect using '>'
  4.  
  5. import sys
  6. import re
  7.  
  8. pat=re.compile('(\S+)\s+(\S+)')
  9.  
  10. f=open(sys.argv[1])
  11. header=f.readline()
  12. if not header.startswith('CLUSTAL') and not header.startswith('MUSCLE'):
  13.     sys.exit("Wrong input file type.")
  14. ln=f.readline()
  15. while not ln.strip():
  16.     ln=f.readline()
  17. seqdict={}
  18. EOF=False
  19. while ln:
  20.     m=pat.search(ln)   
  21.     if m and not '*' in ln:
  22.         seqdict[m.group(1)]=m.group(2)
  23.         ln=f.readline()
  24.     elif not ln:
  25.         EOF=True
  26.         break
  27.     else:
  28.         break
  29. if not EOF:
  30.     while True:
  31.         ln=f.readline()
  32.         if not ln:
  33.             break
  34.         m=pat.search(ln)
  35.         if m and not '*' in ln:
  36.             seqdict[m.group(1)]+='\n'+m.group(2)
  37. f.close()
  38. for k,v in seqdict.items():
  39.     print ">%s\n%s"%(k,v)

输入文件

写得有点混乱,希望大家如果使用的话把BUG报告一下^^

相关文章



Jun 05 2008

收集一堆DNA和蛋白质序列比对软件

分类: studyazalea @ 2:46 am 475次阅读

最近好像爱上收集软件了(其实不是的。。),今天就来收集一下DNA和蛋白质序列比对软件(Sequence Alignment Tools).

最最最著名的就是BLAST了,BLAST虽然叫做the Basic Local Alignment Search Tool(基本局部比对搜索工具),但是研究中经常用BLAST来寻找同源序列。网络版的BLAST非常强大,http://www.ncbi.nlm.nih.gov/blast/Blast.cgi,想必大家都很熟悉,就不多说了。BLAST也可以在本地进行,方法详见如何在本地BLAST, 另外有个本地版的程序叫bl2seq,用来进行2个序列的比对。

FASTA是和BLAST类似的工具,可以在一个数据库中搜索某序列的相似序列,也可以进行2个序列的比对。FASTA比BLAST古老,似乎现在很少有人用,弗吉尼亚大学有FASTA的web server.

一个很古老的序列全局比对(global alignment)工具叫Needle这里有个web版的Needle比对工具。 与Needle对应的是Water,是个很古老的序列局部(local alignment)比对工具,全局比对和局部比对的区别在于,前者尽可能使2个序列在全部长度上可以排列比对,而后者追求达到局部最优比对。这里同时提供了needle和water两种算法进行比对。needle和water只能用于2个序列的比对。题外话:今天得知实验室的师兄用的是Needle,真是怀旧啊。。

另外一个很古老但至今仍然很流行的多序列比对(multiple sequence alignment)工具叫Clustal http://www.clustal.org/. Clustal的图形用户界面版软件叫ClustalX,而命令行模式的版本叫ClustalW,Ubuntu下可以轻松安装CulstalW,只需sudo apt-get install clustalw. Clustal的原理是渐进比对:“在比对过程中,先对所有的序列进行两两比对并计算它们相似性分值,然后根据相似性分值将它们分成若干组,并在每组之间进行比 对,计算相似性分值。根据相似性分值继续分组比对,直到得到最终比对结果。在比对过程中,相似性程度较高的序列先进行比对而距离较远的序列添加在后面。”(抄自这里)因此Clustal不仅可以用来进行多序列比对,也可以构建邻接树(NJ tree).

一个很新的多序列比对工具叫MUSCLE (multiple sequence comparison by log-expectation) http://www.drive5.com/muscle/. MUSCLE对Clustal的算法进行了改进,具体我也不懂,作者测试的结果是与Clustal相比,运行速度和准确性都有所提高。

最近一直在用ClustalW和MUSCLE,下面对二者的优劣谈谈我的看法。ClustalW可以方便的设定Gap Penalty,即对序列比对出现的缺口如何罚分,这个参数对比对的结果有很大影响,但是ClustalW只能输出clw格式,不方便程序的处理。而MUSCLE可以输出多种格式的文件,总有一款适合你-.- 而且MUSCLE可以在已经进行好的比对里增加新的序列,再次比对,而不必从头对这些序列一起比对。

以上软件除了可以设定一些参数外,不能人为对比对结果进行修改。那么如果我们已知这些序列的某段子序列是比对在一起的,而软件计算的结果却不在一起,我们怎么办呢?换句话说,如果我想预先强制确定某个序列的某个位点与另外一个序列的某位点是比对在一起的,那么就需要这个软件:DiAlign, http://bibiserv.techfak.uni-bielefeld.de/dialign/ DiAlign可以输入一个文件,包含每条序列的锚定位点(即强制比对的区域),在此基础上比对其余序列。上述网站提供了web server,本地安装的方法在这里

此外今天看到一个满可爱的web序列比对软件,Base by Base, http://athena.bioc.uvic.ca/workbench.php?tool=basebybase&db=。它是专门为病毒基因和蛋白质序列比对而开发的,不过推而广之也可以用来比对任何序列。使用者可以手动拖动序列,改变比对的位点,也可以选中想要重新比对的子序列,使用T-Coffee, ClustalW或MUSCLE算法重新进行比对。具体使用方法详见the OpenHelix Blog,有视频演示,非常赞。

罗嗦了一大堆,本来还想写一下各个软件适用的情况,不过我也不是很清楚,暂时算啦,留待以后补充。在MUSCLE的文档里看到,如果不能确信比对结果,就去搞几个不同的比对软件,所有结果中共同的比对区域可信度就高。。 我个人的偏好是2个序列比对或数据库搜索用BLAST,多序列比对用MUSCLE,因为选项丰富,可以根据需要改变一些参数,另外就是可以大规模处理数据。

相关文章



Jun 03 2008

收集一堆生物分子三维结构浏览软件

分类: studyazalea @ 11:41 pm 280次阅读

Crystal Structure of Circadian Clock Protein KaiA

上图是某海洋细菌Synechococcus elongatus的生物钟蛋白质KaiA的分子结构。 可惜这幅图是二维的,如果想从其他角度看这个蛋白质分子,就需要三维结构浏览器啦。

首先有一堆网页版的浏览器,比如KiNG, JMol, 这些软件基于web浏览器,不用本地安装,但是都需要Java支持。在http://www.rcsb.org/pdb/explore.do?structureId=1r8j的右侧图片下面有"Display Options",点击下面的各个brower就可以使用了。

以下介绍本地版的分子三维结构浏览器。

RasMol

http://openrasmol.org

RasMol是最著名的分子三维结构浏览器,由Roger Sayle在1992年开发,之后易手数人进行维护更新。免费开源,可以用鼠标旋转观看,可以以多种模式显示(如球棍模式,空间填充模式,卡通模式等),并可以导出为平面图形文件。

下载在这里,有Windows版本和源代码。Ubuntu可以直接 sudo apt-get install rasmol

下面就用上述蛋白质演示一下。

首先在这个页面点下载图标 ,下载该蛋白质的pdb文件。然后打开RasMol,Ubuntu下自动安装在Applications->Education目录下。或者直接在终端输入rasmol运行。然后打开下载的1R8J.pdb文件,开始胡搞。Display菜单可以选择显示模式,Colours菜单可以选择着色方式,Options菜单可以选择是否显示氢原子、是否显示镜像等。然后Export菜单允许用户把图像导出为平面图形文件。

比如下图是根据温度着色的KaiA分子结构的卡通显示模式。

此图和本文最前面的图显示的是同一个分子的结构,但因为把分子旋转了一下,所以视角不同了。

ICM Browser

http://www.molsoft.com/icm_browser.html

是Molsoft公司开发的免费三维结构浏览软件,可以直接从网上获取分子结构数据,不需要自行下载数据。也提供一些显示选项,但不如RasMol丰富。此外有自动旋转和摇摆显示功能。wikipedia上的介绍比较详细。 下载在这里,如果是Linux下运行,需要在终端输入:安装目录/icmbrowserpro -g,比如我的Ubuntu上是/home/azalea/icm-browser-pro-3.5-1n/icmbrowserpro -g

Cn3D

http://www.ncbi.nlm.nih.gov/Structure/CN3D/cn3d.shtml

没用过,据说名字的意思是See in 3D。以下完全抄自http://bio.med.stu.edu.cn/download/soft.htm

“观看蛋白质三维结构的软件,是由NCBI开发的。其设计的主要目的是为NCBI在其站点中的蛋白质结构数据库MMDB提供专业的结构观察软件。与其他的类 似的软件,如Rasmol,Weblabview等相比,其在结构观察方面主要功能上基本相似,但是图形格式上比Rasmol和Weblabview要差 一些。而在与网络连接上,该软件能依托NCBI所建立的所建立的MMDB结构数据库,能直接根据输入的序列号从数据库中利用其内嵌的Entrez搜索引擎 调出蛋白结构来进行观察,比其他软件要简便。而Cn3D主要的特点是能够将两个蛋白放在一起直观地进行三维结构上的比较,如下图中是将两种核酸外切酶的三 维结构通过VAST对准得到的结构比较图:同样,Cn3D在结构比较方面也能利用其内嵌的Blast搜索引擎直接访问Genbank数据库找到具有局部相 似性的结构数据,并在三维结构图中显示出二者具有相似性的结构区域。”

其他一些流行的分子三维结构浏览软件:

PyMol

http://pymol.sourceforge.net

VMD

http://www.ks.uiuc.edu/Research/vmd

Protein Explorer

http://www.umass.edu/microbio/rasmol/index.html

都不了解,就不多说了,希望用过的人评价一下。

相关文章



Apr 30 2008

Linux安装libsbml

分类: studyazalea @ 7:46 pm 255次阅读

libSBML是用来读写和操作Systems Biology Markup Language (SBML)的库。根据 sbml.org 的描述,"SBML是一种描述生物化学反应模型的机器可读格式,可以应用到新陈代谢、细胞信号转导等模型"。直观地说,SBML是一种类似XML的格式,用于表示系统生物学模型。

libSBML可以在http://sourceforge.net/project/showfiles.php?group_id=71971&package_id=71670下载,下载后解压,进入解压后的文件夹。

Linux下的libSBML下的安装需要以下步骤:

第1步:

./configure

注意事项1:默认情况下libSBML尝试使用libxml2 XML library的2.6.16或更新版本,如果没有这个library,可以用命令:

./configure --with-expat

./configure --with-xerces

Expat, Apache Xerces-C++和Libxml2都是常用的XML parser library.

注意事项2:

默认情况下,libSBML只创建C和C++的API,如果想通过其他软件使用libXBML,需要在configure时加上参数--with-java, --with-lisp, --with-python, --with-perl, --with-matlab, --with-octave, 或 --with-ruby

比如希望Python和Ruby使用libSBML,就用命令:

./configure --with-python --with-ruby

第2步:

make

第3步:

sudo make install

以上不用多说

第4步也是最后一步:

sudo ldconfig

不熟悉的可以man ldconfig (这话不是我说的。。) 这步的目的是使第3步的软件可以在运行时找到libSBML的库文件。

搞定啦。详细的安装说明以及Windows和Mac OS的安装说明请看http://sbml.org/Software/libSBML/docs/python-api/libsbml-installation.html

可以简单测试一下:

python

>>>import libsbml

相关文章



Apr 27 2008

生物学家能修好收音机么

分类: fun, studyazalea @ 1:18 am 353次阅读

这个问题是Yuri Lazebnik (以下简称Y)在这篇paper: "Can a biologist fix a radio?—Or, what I learned while studying apoptosis"中探讨的。

Y是研究细胞凋亡的生物学家,他刚出道时,向著名的David Papermaster请教,D说,一个生物领域的兴起好比Gold rush,而他进入细胞凋亡领域时就如同真相大白大家发现金子没那么容易淘的低迷时期。一个悖论就是,虽然许多许多人在淘金,而且金矿确实存在,但是并不保证能淘到金。D说,如说你想在举世沮丧的环境下继续有意义的研究,就要学会使用好的工具,保持清醒的头脑。

Y于是听了D的话,结果研究了一段时间,又迷惑了。Y不明白没啥研究了这么久,数万论文发表了,细胞凋亡或细胞周期的途径,甚至一个蛋白质p53的功能却似乎越来越深奥。

于是Y又请教了Joe Gall,J说事物发展的规律就是这样,停滞一段时间就会继续发展,J举例说明道,19世纪就有对细胞死亡的研究,沉寂了一个世纪,而今中兴,10年内发表了60000篇论文。Y听了虽然如释重负,但仍觉得,细胞凋亡研究的悖论说明生物学家解决问题的方法有问题。

那么问题在哪里呢?前面都不用看,现在才进入正题:

Y的高中数学老师讲过,检测一个方法可行性的最佳手段就是把这个方法应用到一个答案已知的问题上。于是Y就想,不如让生物学家来修理收音机,如果能修好就说明生物学的研究方法不错。因为收音机和生物学里的信号转导途径相似,都是把信号从一个形式转换为另一种形式。于是开始想像之旅:

首先,生物学家发现一台收音机挂了(简单起见,就比如上图的老式收音机,不能再播放电台的音乐了),于是他们通过资助搞来一大堆相同的但是正常工作的收音机,肢解后与挂掉的收音机比较。生物学家终于学会打开收音机,并发现一堆形状颜色大小各异的零件。

然后,生物学家每次移除一个零件,看看有什么表现型(即收音机哪里不再正常工作)。 虽然移除很多零件并无大碍,但是一个幸运的博士后偶然发现一条电线如果挂了,收音机就不播放音乐了。于是这个喜出望外的博士后把这条电线命名为"偶然发现组件"Serendipitously Recovered Component (SRC) ,然后发现SRC是连接一个长长的零件和收音机其他部分的唯一零件,因此这个长零件理所当然被命名为"最重要组件" the Most Important Component (MIC)"。然而,一个研究生改变了MIC的长度,发现并不显著影响播放音乐的质量,而他发现了一个重要组件被称为"真正重要组件" Really Important Component (RIC)。于是关于MIC和RIC谁重要的争论此起彼伏,直到一个聪明的博士后发现一个开关,它的状态决定播放音乐依赖于MIC还是RIC。于是自然这个开关就得名为"无疑最重要组件" Undoubtedly Most Important Component (U-MIC)。。。

通过切断一个组件与其他的关系进行研究的方法会提供无数关于连接性的信息,但是如果收音机有可调控组件,这个方法就未必有效了。收音机不正常工作可能是因为一些组件没有被调好(比如把音量调小到你听不到的程度),而对此收音机组件的连接性质不能提供任何有价值的信息。

然而我们知道几乎任何一个修理工都能修好收音机。为啥区别这么大呢?Y认为是由于生物学和工程学上使用的语言不同,下图表示得很明白:

恩,最后强烈建议大家去读一下原文,绝对不会睡着的,而且发人深省

还要感谢Y的paper里的图,还望原谅我未经允许就贴过来了。bow~~

相关文章



第 1 页,共 3 页123>