生物信息学常用软件 李蔚 Co-head, Bioinformatics Department Beijing Genomics Institute liw@genomics.org.cn

目录 什么是生物信息学？ 大规模基因组测序中的信息分析 基因组注释（annotation）软件介绍 其他常用软件介绍 应用实例 文献查阅方法概述 大规模基因组测序中的信息分析 基因组注释（annotation）软件介绍 其他常用软件介绍 应用实例

广义地说，生物信息学从事对生物信息的获取、加工、储存、分配、分析和释读；并综合运用数学、计算机科学和生物学工具，以达到理解数据中的生物学含义的目标。与此相应，生物信息学具有三方面的科学基础 它需要发达的、复杂的、可相互交流的数据库系统。 其次，生物信息学需要强有力的创新算法和软件。 自动化的大规模高通量的生物学研究方法与平台技术。

DNA(基因组学, exon/intron, SNP) Bio-technology (HGP、RGP、PGP) Info-technology BIT 生物信息学/BioX Database Data analysis Bioinformatics software Data service 转录 mRNA (cDNA, EST, Gene, SNP) 表达谱，生物芯片 翻译 基因芯片 蛋白芯片 Protein (蛋白质学，2D system,Sequence) 3D structure、 Folding、Drug design

专业文献管理软件EndNote 5.0 Search bibliographic databases on the Internet Organize references and images in a database Watch the bibliography and figure list appear as you write!

搜索、管理与课题相关的所有文献 http://www.ncbi.nlm.nih.gov/entrez/query.fcgi?db=PubMed，

B1 Related Articles References B2 cited by other articles http://wos.isiglobalnet.com Times Cited

搜索、管理与课题相关的所有文献 http://www.facultyof1000.com/ B-3 BioMail 每周定期自动检索美国NCBI的文献资料库, 即PubMed® Medline®数据库，你可以自定义你的检索词，我们定期将检索结果用email提供给你。 这样你就可以定期收到最近发表在专业期刊上的最新文章摘要。我们相信这对于国内科研工作者是很有帮助的。

搜索、管理与课题相关的所有文献 B-4 Bioinformatics and Systems Biology Electronic Library http://genome.ucsd.edu/library/index.php

搜索、管理与课题相关的所有文献 B-5 SCI咨询中心 http://www.lib.tsinghua.edu.cn/NEW/home810.htm http://159.226.159.166/factor/

搜索、管理与课题相关的所有文献 B-6 Email list of Journal TOC

免费获取全文 (1) 免费全文杂志 http://china.sciencemag.org/ http://www.pnas.org/ http://www.dna-res.kazusa.or.jp/

免费获取全文 (2) 全文信息网站 http://intl.highwire.org/

免费获取全文 (3) http://www.freemedicaljournals.com/htm/index.htm 提供小语种、多种学科的免费外文原文，我觉得和上面介绍的High Wire Press 不相上下。都是巨好的免费外文全文提供网站。

免费获取全文 (4) 根据作者E-mail地址，向作者索要。这是一种看似麻烦，实际非常有效的办法。有人写了索取原文模板，您只要把作者姓名、文章标题等等加上就成了。天天是不是很体贴人呀？ 呵呵~~~~~~~ 拜托您使用时稍加改动，避免千篇一律。信的内容如下：

免费获取全文 (5) 如果作者有自己的主页，可以去作者的主页看看

免费获取全文 (6) 北大图书馆 http://162.105.138.23/reference/reference.asp

免费获取全文 (7) Search as you wish: 直接用文章名作关键词搜索。http://www.google.com

免费获取全文 (8) 虫http://www.scienceworm.com/paper.asp

Sequence Data Distribution 大规模基因组测序中的信息分析 Phred Phd2fasta Crossmatch Phrap Consed Repeatmasker Glimmer Blastx Blastn Clastal tRNAscan Sequencing Basecalling Vector Quality Control QualCal QualDraw QualStat Vector Mark Assemble Finishing Repeat Mark ORF Prediction Gene Annotation COGs Swiss-port, PIR, GDB, GenBank Sequence Data Distribution

Sequencing 在文库建立完成以后，就需要对短克隆片断文库进行高通量的测序了。当然，也要适当的加入长克隆片断的双向测序，以利于下一阶段的工作。一般而言，随着测序量的增加，通过软件组装得到的连接群的数目将呈现一个先增加后减小的变化过程（如图5.1A），这就是著名的Lander-Waterman曲线（Lander et al. 1988）。因此，理论上说，在测序量达到大约4～5倍的覆盖率（converage）时，也就是说获得基因组大约95%的序列，一般就可以通过计算机软件对它们进行拼接，从而获得此基因组序列的工作草图（working draft）。然后，继续随机测序的过程，到大约8～9倍的覆盖率左右（一般可以获得基因组99%的序列），此时连接群的数目已经比较少了，而且从Lander-Waterman曲线上可以看出，再增加测序量时，连接群的数目下降的将非常慢。因此，这意味着随机测序的阶段已经不能继续进行了，必须转入下一阶段。这里虽然已经获得了基因组99％以上的序列，但是由于此时得到的还只是一些彼此孤立的连接群，它们之间还存在着缺口（gap），因此距离序列拼接的最终目标——“完成图”还有一定的距离。

Basecalling Phred是一个采用快速傅利叶变换（fast Fourier transform）分析技术以及动态规划算法（dynamic programming algorithm）从DNA测序所得到的图形数据中提取DNA碱基排列顺序信息（即Base-Calling），从而得到DNA序列的软件。Phred 对序列中的每一个数据产生一个被广泛接受的带有质量控制标准（quality scores）的“Base Call”。Phred质量指标x就相当于约10-x/10的误差概率。因此，Phred质量指标20就相当于在原始数据中一个Base Call的精确度为99%。具体的说，Phred可以读取DNA测序仪生成的色谱图（chromatogram）文件（二进制格式），通过分析每个碱基的“质量”信息而输出每个测序序列的碱基序列和质量信息文件（文本格式）。 Phred程序在使用前需要将软件包中Phred的参数文件phredpar.dat拷贝到一个所有的Phred的用户都可以访问的目录中（比如说/usr/local/etc/PhredPar）。然后设置一个环境变量“PHRED_PARAMETER_FILE”来指向此文件的位置。比如，在C shell中我们可以使用如下命令进行设置： > setenv PHRED_PARAMETER_FILE /usr/local/etc/PhredPar/phredpar.dat

phred程序常用参数表 （所适用版本0.990722j） 类型 参 数 说 明 输入参数 -id <dir_name> 读取并处理由<dir_name>指定的目录中的所有文件 -if <file_name> 读取并处理列在<file_name>文件中的所有文件 -zd <dir_name> 指定压缩程序所在的位置（目录） -zt <dir_name> 指定解压缩时临时文件的存放目录 处理参数 -trim_alt <enzyme_seq> 使用所给定的序列对输入序列进行“修整”，即根据质量指标去除掉输入序列两端与给定序列一致的部分。并且将低于相应阈值（可由-trim_cutoff指定）的数据质量指标直接设置为0。可以用””将指定序列设置为空 -trim_cutoff <value> 为“-trim_alt”选项设置进行“修整”时的阈值（即错误率），默认为0.05。并且此阈值将被写入Phd输出文件中 输出参数 -st fasta 设置输出序列文件格式为Fasta(默认) -sd <dir_name> 设定输出的序列文件名为原来输入文件名后面加“.seq”后缀，并将他们全部保存在由<dir_name>指定的目录中 -sa <file_name> 设定将此次运行Phred所处理的全部输入对应的序列以Fasta格式写入由<file_name>指定的单个输出文件中 -qt fasta 设定输出的质量文件格式为Fasta(默认) -qd <dir_name> 设定输出的质量文件名为原来输入文件名后面加“.qual”后缀，并将他们全部保存在由<dir_name>指定的目录中 -qa <file_name> 设定将此次运行Phred所处理的全部输入对应的质量以Fasta格式写入由<file_name>指定的单个输出文件中 -qr <file_name> 产生一个关于所有输入的“高质量”碱基数目的柱状图统计数据 -pd <dir_name> 设定输出的Phd文件名为原来输入文件名后面加“.phd”后缀，并将他们全部保存在由<dir_name>指定的目录中 -dd <dir_name> 设定输出的多态性文件名为原来输入文件名后面加“.poly”后缀，并将他们全部保存在由<dir_name>指定的目录中 其他参数 -v <n> 显示详细的处理过程，可以通过设置n的不同大小来决定屏幕输出信息的多少（默认为0，即最少的屏幕输出） -h, -help 显示Phred的一个命令行参数的摘要 -doc 显示Phred的详细说明文档 -V 显示Phred的版本信息

Vector Mark 负责标记载体序列（mark vector）的cross_match程序 cross_match和phrap程序部分常用参数表（所适用版本0.990329） * 表示只phrap程序可用的参数 参 数 说 明 -penalty n 设置使用Smith-Waterman比对算法（以下简写为SWAT）时的错配（即出现碱基替代）罚分为n（默认值为-2） -minmatch n 设置SWAT比对时的最小片断长度为n（默认值为14） -minscore n 设置SWAT匹配时的最低分值为n（默认值为30） -screen 主要在cross_match中使用，此时将生产一个“.screen”的Fasta格式的输出序列文件。且在此文件中seq1原来序列中与载体序列（seq2等）匹配的部分已经被用X替代 -trim_qual n* 定义在组装时所使用的“高质量”的标准为n（默认值为13） -forcelevel n* 设置在拼接生成contig时的严格程度为n（默认值为0，即最严格。最大可以去10） -revise_greedy* 将最初的“贪婪”拼接在连接较“弱”处“扯开”，而重新考虑全局的最优匹配 -view* 生产一个“.view”的文件（为便于使用PhrapView软件来查看全局拼接情况） -new_ace* 设置用新的“.ace” 文件格式输出，它包含更多信息，且只有原来大小的约60％，可用于高版本（1998年以后）的Consed软件

Assemble (1) Phrap是一个用于将鸟枪法测序的原始序列拼接成连接群的软件，它的核心仍是Smith-Waterman算法。但是，有所不同的是，Phrap是结合了相应的质量控制标准值对Phred程序Base-Calling所得到的带有质量信息的DNA序列进行拼接的。它首先搜索具有匹配单字的反应序列，删除掉明显是完全一致的序列，计算出有匹配的两两序列的匹配分值，标记出可能是一致的反应序列、自身有匹配的反应序列、在有匹配的序列对中可信度不高或缺少高分值匹配片段的序列，并给出可能的嵌合体或有缺失的序列。然后根据以上对序列的标记，利用一致性好的匹配片段来构建一致序列的框架，并利用标记出的匹配来确认序列的可信部分，更新原有的质量分值。

Assemble (2) 另外，使用phrap进行拼接，对输入文件还有一些特殊的要求：①Reads的命名规则，phrap需要从文件名了解此反应序列所属的亚克隆、测序方向和测序所用的化学试剂等信息，因此为了更有效的使用它，必须遵循其所设定的命名规则[LW1]；②序列（sequence）文件，一般phrap使用单个的输入序列文件，其中包含了所有的reads的序列，而cross_match则使用两个输入序列文件，前一个是query序列（一般就是所有reads的序列），后一个是subject序列（一般是载体序列），同时所有序列文件都必须是Fasta格式；③质量（Quality）文件，质量文件必须与输入的序列文件放在同样的路径中，并且命名应该与之相同，而只是增加了“.qual”的后缀，同时每个碱基的质量指标值应该在0到99之间；④载体（Vector）序列文件，此文件只是用于cross_match程序（作为第二个输入参数），其中应该包含了序列扩增时所用载体的序列，以便在进行正式的序列组装之前将其屏蔽（screen）。

Finishing (1) 这一阶段通常又称为“Finishing阶段”，是基因组测序组装中最困难的阶段之一。因为，此时由于前面所说的种种原因，已经无法通过增加随机测序片断的数目来减少连接群的数目了。这时，就必须充分的利用长克隆片断文库中已经进行了双向测序的片断的每一对正反向的reads之间的链接（link）关系。也就是说，考察每一对这样的正反向reads，看它们是否都参与了拼接？并且是否存在于目前的同一个连接群内部？如果是，则可以通过分析它们两者之间的“距离”是否在10kb左右，而用以在一定程度上验证拼接的正确性；否则，就可以初步确定出它们分别所在的这两个连接群之间的位置关系了（如图5.2）。当然，如果一对短克隆片断的正反向reads恰好落在两个不同的连接群内（一般应该在末端），按照同样的原理，也可以用来确定这两个连接群之间的位置关系。

Finishing (2) 在确定了连接群之间的位置关系以后，通常对于不太大的缺口，可以通过Primer Walking的办法来将其补上。所谓Primer Walking，即按照上一个测序结果在其中400-450碱基左右的位置设计引物继续测序（因为引物区后必须还有足够的准确序列以便拼接，所以设计引物的位置就必然比较靠前），并将测序结果拼接（如图5.3）。另外，要补一些比较大的缺口，有时还可以采用以基因组DNA序列做模板而进行的长程PCR（long-range PCR）的办法来进行（这可能需要使用一些现成的试剂盒）。

Finishing (3) 在这一阶段的很多工作，通常也需要借助一些软件工具来完成，而且事实上这一类的工具现在也有很多种。但是，如果在前面的拼接阶段使用的是Phred-Phrap软件包，那么这里最佳的选择必然是与之配套的Consed软件（Gordon et al. 1998）。它是一个基于Unix/Linux系统的专门用来查看和编辑phrap的拼接结果的可视化工具，目前的最高版本是12.0，支持的操作系统包括Solaris、Redhat Linux、Irix、DEC-Alpha和HP-UX等。由于Consed软件也是华盛顿大学开发的，所以它从开发阶段起就紧密的与phrap结合——Phil Green（Phrap的作者）直接指导了它的设计。因此，只有使用Consed软件才能充分发挥phrap拼接的所有功能（比如有效的读取phrap所产生的一系列标记等）。而且Consed的图形界面使得phrap的拼接结果变得很直观，用户可以直接查看测序色谱图，并可以很容易的检查一些低质量区域、错拼或由于Base-Calling的误差而造成的不匹配等等，甚至还可以直接根据连接群的序列设计PCR引物等等。另外，在后期的Consed版本中还加入了一个新的称为Autofinish的功能，可以很大程度的简化补洞的工作量并降低Finishing阶段的成本（Gordon et al. 2001）。

基因预测软件介绍 基于人工神经网络和模式识别的算法和软件 基于语言学方法 基于隐含马尔可夫模型 GRAIL（http://compbio.ornl.gov） GeneParser （http://beagle.colorado.edu/~eesnyder/GeneParser.html） GENEID http://www1.imim.es/software/geneid/index.html 基于语言学方法 Genlang http://www.cbil.upenn.edu/genlang/genlang_home.html 基于隐含马尔可夫模型 Genie http://www.cse.ucsc.edu/~dkulp/cgi-bin/genie） HMMgene http://WWW.cbs.dtudk/services/HMMgene/） GeneFinder 软件由美国华盛顿大学的Colin Wilson，LaDeana Hilyer,和Phil Green研发 GENSCAN 软件由斯坦福大学数学系Chris Burge和Samuel Karlin所研发

基因预测软件介绍 对于细菌等原核生物的可读框的识别，目前理论手段已经发展的比较成熟，而且特异性（specificity）和灵敏度（sensitivity）都相对较高。其中使用最为广泛的软件是由TIGR中心开发的Glimmer（Salzberg et al. 1998a），这个名称是“基因定位和内插马尔可夫模型”（Gene Locator and Interpolated Markov Modeler）的缩写，它是专门用于在微生物基因组中定位基因的软件，对于细菌（bacteria）和古细菌（archaea）特别有效。同样的，Glimmer也是对学术研究机构免费的，详细的信息和如何获取可以从TIGR的网站上的相关网页中了解到（http://www.tigr.org/software/glimmer/）。目前可以得到的是它的改进版本2.02（Delcher et al. 1999）的源代码包，解开编译（make）以后将产生四个独立的可执行程序：①long-orfs，用于根据启动子和终止子的位置，给出可能的范围较长且重叠较少的阅读框；②extract，用于根据指定的信息从序列文件（FASTA格式）中提取特定的子序列；③build-icm，用于建立并输出内插马尔可夫模型（IMMs）；④glimmer2，用于根据给定的序列文件和内插马尔可夫模型（由build-icm程序产生的）列出所有可能的阅读框。 具体的参数和用法可以参阅其软件包中提供的针对每个程序的名为“*.readme”的文档。在实际使用时可以直接运行其提供的简单的shell脚本run-glimmer2，也可以根据各自的需要调节参数或单独运行。

本地使用Genscan也很简单，其用法如下： > genscan <parameter_file> <input_file> [options] 其中，genscan是可执行文件的名称，parameter_file是所使用的参数文件的路径和名称（不可省略！），input_file是要预测的输入文件的路径和名称（可以是Fasta 或 GenBank 格式），options是输入选项。Genscan 自身附带了三个参数文件：HumanIso.smat用于人和脊椎动物的序列（果蝇也可以）；Arabidopsis.smat用于拟南芥的序列；Maize.smat用于玉蜀黍属的序列。另外，Genscan只有四个输入选项：① -v，设置在文本输出中附加上解释信息；② -cds，设置在输出预测出的氨基酸链的同时，输出对应的核酸序列；③ -subopt <cutoff>，设置同时显示那些预测的概率大于指定阈值cutoff（最小值是0.01），但未达标准的（suboptimal）外显子；④ -ps <out.ps> <N>，设置以每行N个碱基的比例生成一个PS格式的输出文件out.ps（任选，缺省的文件名是在输入的文件名后加“.ps”后缀），其中显示了预测出的基因和外显子的分布情况。

基因功能预测 一般可以分为至少两个层次进行。 首先是根据一个基因整个序列与数据库中已知基因的序列的相似性来确定其功能； 其次，当无法找到与其完整序列具有相似性的已知基因时，还可以根据其中的一些结构域（domain）的特征来预测它们的功能。

a) BLAST（Basic Local Alignment Search Tool） 将在新基因组中预测出的所有的基因在公共数据库，比如GenBank所提供的非冗余蛋白数据库（non-redundant protein database，简称NR），进行检索。 NR蛋白质数据库可以从NCBI的Ftp服务器上直接下载（ftp://ftp.ncbi.nih.gov/blast/db/nr.Z），而BLAST软件包则需要根据所要使用操作系统而选择的相应的编译过的可执行程序下载（ftp://ftp.ncbi.nih.gov/blast/executables/）

Blast-like Tools MuSeqBox （blast 结果处理程序，可以设定 e值，score值等等） http://bioinformatics.iastate.edu/bioinformatics2go/mb/MuSeqBox.html PatternHunter http://www.bioinformaticssolutions.com/products/ph.php SSAHA http://www.sanger.ac.uk/Software/analysis/SSAHA/ Blat http://genome.ucsc.edu/

b) Domain search 用HMMER软件包对Pfam数据库进行检索。 http://hmmer.wustl.edu/ Pfam数据库（Protein families database of alignments and HMMs）是由英国的Sanger中心（The Wellcome Trust Sanger Institute）在1998年建立并维护至今的(Sonnhammer et al. 1998; Bateman et al. 2002)，其中收集了大量通过多序列比对产生的在一定范围具有隐马尔可夫模型（Hidden Markov Models，简称HMMs）的蛋白质结构域和蛋白质家族

Sequence Viewer Artemis http://www.sanger.ac.uk/Software/Artemis/

比较基因组研究 “近亲”的物种之间的差异非常小，有些在基因组序列水平的相似性就高达90％以上。但尽管如此，那些不足10％的序列差异却导致了它们之间的很多本质的变化，如种属的差异或致病性的有无等等。因此，对它们进行比较基因组的研究可能可以为人们理解物种间差异的本质、基因的横向转移在进化中的作用和特定微生物致病的机制等提供一条捷径。 例如，威斯康星基因组中心通过对大肠杆菌的致病株O157:H7株和非致病株K-12株进行基因组序列比较发现，它们之间共有了一个约4.1Mbp的高度同源且几乎线性一致的“主干”基因组序列（只在复制的终止点附近出现了一个约422Kbp的倒置），而这个主干序列则被几百个具有株特异性（stain-specific）的片断序列所“切断”,并且在O157:H7株的这些特异性的序列片断中，他们发现了很多可能的毒力相关的基因和一些前噬菌体的基因等，从而为研究O157:H7株的致病性机理提供了很大的帮助

比较基因组常用软件 ACT 可视化的显示功能 http://www.sanger.ac.uk/Software/ACT/ PipMaker(Schwartz et al., 2000) - http://bio.cse.psu.edu MUMer(Delcher et al., 1999) Using the following algorithmic techniques: Suffix tress, Special variant of sequence alignment known as longest increasing subsequence, Regular Smith-Waterman alignment DBA (Jareborg et al., 1999) GLASS (Batzoglou et al., 2000) WABA (Kent and Zahler, 2000) Dialign (Morgenstern et al., 1998; G¨ottgens et al., 2001) Alfresco (Jareborg and Durbin, 2000) LAJ (Wilson et al. 2001)

物种进化分析 华盛顿大学（University of Washington）的Joseph Felsenstein等开发的PHYLIP（PHYLogeny Inference Package）软件包是目前在进化分析中使用最广泛的工具之一。它从1980年就开始发布，并不断的增加修改完善，当前的版本（v3.6）共包含了34个程序，这些程序基本上覆盖了进化分析的所有方面。另外，PHYLIP是一个免费的软件，从其主页上可以直接下载源代码（http://evolution.genetics.washington.edu/phylip.html），

应用实例 我们使用pUC18作为克隆载体，选取大小在1.5kb和10kb的两种插入片断。同时为保证文库的随机性，我们使用了两种不同的方法构建此文库——超声和酶切（Sau3AI）。 大规模随机测序阶段，共产生了约75 000个有效的reads（这里我们设定了phred程序的过滤质量为20，并只保留有效读长大于50bp的reads），其中包含了约2000个长片断克隆的正反向测序reads，覆盖率接近10倍（这里要测序的微生物预计基因组大小约2.7Mbp）。基因组测序的数据的Base-Calling和拼接，我们使用的是Phred-Phrap软件包，主要命令及参数如下：

> phred –id chromat_dir –pd phd_dir –trim_alt "" –trim_cutoff 0 > phred –id chromat_dir –pd phd_dir –trim_alt "" –trim_cutoff 0.01 > phd2fasta –id phd_dir –os assemble/name.seq –oq assemble/name.seq.qual > cross_match name.seq vector.seq –minmatch 12 –minscore 20 –screen > screen.out > phrap name.seq.screen –trim_qual 20 –revise_greedy –forcelevel 1 –minmatch 17 –minscore 40 –view –new_ace > phrap.out 对于基因的预测，我们主要使用了Glimmer2.0（基本为默认参数），命令如下： > long-orfs genome.fna | get-putative > train.coord > extract genome.fna train.coord > train.seq > build-icm < train.seq >train.model > glimmer2 genome.fna train.mode | get-putative > g2.coord 随后，我们使用如下的命令将它们与一个本地的NR蛋白质序列库进行比对检索： > blastall –p blastp –i orfs.faa –d nr –e 0.01 –o Orfs_Vs_NR.out –F F –a 2 –b 1 –v 1 而对于那些没有找到可信比对的序列，我们还将它们单独挑选出来用如下的命令在本地的Pfam数据库进行了检索： > hmmpfam -A 1 -E 0.01 Pfam NomatchWithNR_orfs.faa 对基因组上的tRNA，我们采用如下的命令和参数进行识别： > tRNAscan-SE –B –H –o tRNA.out –f tRNA_struc.out –d genome.fna 使用NCBI提供的专门用于序列提交的工具Sequin将其转化为ASN.1格式后提交到了GenBank公共数据库

Circular representation of the genome of T. tengcongensis MB4