——本笔记来自山东大学巩晶教授生物信息学课程的整理

第一章 生物数据库

生物数据库首先分成三大类。核酸数据库，蛋白质数据库和专用数据库。核酸数据库顾名思义，是与核酸相关的数据库。蛋白质数据库是与蛋白质相关的数据库。专用数据库是专门针对某一主题的数据库，或者是综合性的数据库，以及无法归入其他两类的数据库。 核酸数据库和蛋白质数据库又分为一级和二级。一级数据库存储的是通过各种科学手段得到的最直接的基础数据。比如测序获得的核酸序列，或者 X 射线衍射法等获得的蛋白质三维结构。一级核酸数据库，他主要包括三大核酸数据库和基因组数据库。 蛋白质的一级数据库还可以再具体分为蛋白质序列数据库和蛋白质结构数据库。二级数据库是通过对一级数据库的资源进行分析、整理、归纳、注释而构建的具有特殊生物学意义和专门用途的数据库。

一、 文献数据库PubMed

PubMed是拥有超过两百六十万生物医学文献的数据库。这些文献来源于 MEDLINE，也就是生物医学文献数据库、生命科学领域学术杂志以及在线的专业书籍。他们大部分提供全文链接。 PUBMED 搜索结果页面，在显示内容格式这个下拉菜单里，除了总结，摘要，还有个叫 MEDLINE 的项目。你可以把它简单的理解为数据库中文献记录的内部结构。比如 TI 代表题目，AB 代表摘要，AU 代表作者。比如搜索 Down 这个词。我们在 Down 的后面加上空格，中括号 AU（Down [AU]），就会返回所有作者名里有 Down 这个词的文献。如果加上[TI]，则返回题目中有 Down 的文献。中括号 AD 是搜索发表单位。如果什么限制都没有，只写 Down 的话是在任意地方搜索。另外，PubMed 还有几点要说明的。在搜索时，可以使用引号，引号里的词会被当作一个整体来看待，而不会被拆开；可以使用逻辑词AND、OR 、NOT。比如你可以规定，题目里有 dUPTase，并且题目里有 bacteria，但是作者里不要有 Smith（dUTPase [TI] AND bacteria [TI] NOT Smith [AU]）。每篇文章都有自己唯一的 PubMed ID（PMID）。通过这个号，可以直接找到某一篇文章。最后，不得不说的是，有的时候 PubMed 也帮不了你。比如，搜索 1995 年以前的文献中排名十位以后的作者是白费力气。搜索 1976 年以前的文献是没有摘要的。搜索 1965 年以前的文献，就别想了。

二、 一级核酸数据库

三大核酸数据库包括NCBI 的 Genbank，EMBL 的 ENA 和DDBJ，它们共同构成国际核酸序列数据库。 NCBI隶属于美国国立卫生研究院 NIH； 欧洲核苷酸序列数据集 ENA 由欧洲分子生物学研究室 EMBL 开发并负责维护。 亚洲的核酸数据库 DDBJ由位于日本静冈的日本国立遗传学研究所 NIG 开发并负责维护。 Genbank，EMBL 与 DDBJ共同构成国际核酸序列数据库合作联盟 INSDC。

（一）GenBank 原核生物核酸序列

从 NCBI 的主页（http://www.ncbi.nlm.nih.gov/）选择Genbank 数据库。Nucleotide数据库就是 Genbank 数据库。以一条原核生物的 DNA 序列为例解读序列信息。这个序列是编码大肠杆菌 dUTPase 的基因，在Genbank 里的数据库编号是 X01714。搜索结果解释说明： LOCUS 这一行里包括基因座的名字，核酸序列长度，分子的类别，拓扑类型，原核生物的基因拓扑类型都是线性的，最后是更新日期。 DEFINITION 是这条序列的简短定义，也就是前面看到的标题。 ACCESSION 就是在搜索条中输入的那个数据库编号，也叫做检索号，每条记录的检索号在数据库中是唯一且不变的。即使数据提交者改变了数据内容，Accession 也不会变。你会发现，这条记录里，Accession 和 Locus 是一样的。这是因为这个基因在录入数据库之前并没有起名字，因此录入数据的时候便将检索号作为了基因的名字。但是有些基因，在录入数据库之前已经有了自己的名字，那么这些基因所对应的 Accession 和 Locus 就不一样了。 Version 版本号和 Locus，Accession 长得差不多。版本号的格式是“检索号点上一个数字”。主要用于识别数据库中一条单一的特定核苷酸序列。在数据库中，如果某条序列发生了改变，即使是单碱基的改变，它的版本号都将增加，而它的 Accession 也就是检索号保持不变。比如，版本号由 U12345.1 变为 U12345.2，而检索号依然是 U12345。版本号后面还有个 GI 号。GI 号与前面的版本号系统是平行运行的。当一条序列改变后，它将被赋予一个新的 GI 号，同时它的版本号将增加。 KEYWORDS 提供能够大致描述该条目的几个关键词，可用于数据库搜索。 SOURCE 基因序列所属物种的俗名。他下面还有一个子条目ORGANISM，是对所属物种更详细的定义，包括他的科学分类。 REFERENCE 是基因序列来源的科学文献。有时一条基因序列的不同片段可能来源于不同的文献，那样的话，就会有很多个 REFERENCE 条目出现。REFERENCE 的子条目包括文献的作者、题目和刊物。刊物下面还包括 PubMed ID 作为其子条目。 COMMENT 是自由撰写的内容，比如致谢，或者是无法归入前面几项的内容。 FEATURES 是非常重要的注释内容，它描述了核酸序列中各个已确定的片段区域，包含很多子条目，比如来源，启动子，核糖体结合位点等等。 source 说明了核酸序列的来源，据此可以容易的分辨出这条序列是来源于克隆载体还是基因组。可以看到，当前序列来源于大肠杆菌的基因组 DNA。 promoter 列出了启动子的位置。细菌有两个启动子区，-35 区和-10 区。-35 区位于第286 个碱基到第 291 个碱基 ，-10 区位于第 310 个碱基到第 316 个碱基。 misc_feature 列出了一些杂项，比如，这条说明了从第 322 个碱基到第 324 个碱基是一个推测的，但无实验证实的转录起始位置。 RBS 是核糖体结合位点的位置。 CDS： Coding Segment，编码区。对于原核生物来讲，CDS 记录了一个开放阅读框，从第 343 个碱基开始的起始密码子 ATG 到第 798 个碱基结束的结束密码子 TAA。除了位置信息，还包括翻译产物的诸多信息。翻译产物蛋白的名字是 dUTPase，这个编码区编码该蛋白的第 1 到第 151 个氨基酸。翻译的起始位置和翻译所使用的密码本，以及计算机使用翻译密码本根据核酸序列翻译出的蛋白质序列。需要强调的是，这不是生物自然翻译的，而是计算机翻译的。事实上，蛋白质数据库中的大多数蛋白质序列都是根据核酸序列由计算机根据翻译密码本自动翻译出来的。中间部分是翻译出的蛋白在各种蛋白质数据库中对应的检索号。通过这些检索号可以轻松的链接到其他数据库。 ORIGIN 作为最后一个条目记录的是核酸序列，并以双斜线作为整条记录的结束符。至此整条记录就浏览完了。 

（二）GenBank 真核生物核酸序列 mRNA

真核生物的基因与原核生物不同，是非线性排列的，也就是基因里有外显子和内含子。因此真核生物核酸序列的数据库记录要要比原核生物复杂。 我们以编码人 dUTPase的成熟 mRNA 序列为例。成熟 mRNA 是已经剪切掉内含子，只剩外显子的序列，所以这条成熟mRNA 序列和我们之前看到的原核生物的 DNA 序列从拓扑结构上看是几乎一样的，都是线性的。输入这条成熟 mRNA 序列的检索号 U90223，搜索。 打开数据库记录，基本的注释内容和原核生物的差不多，这里只挑两点特别的地方说一下。大家看到 KEYWORDS 后面只有一个点。这个点提示我们，数据库并不是完美的，所有数据库都存在数据不完整的问题。再有，JOURNAL 后面我们看到是写的是未正式发表。但事实上，这篇文章早在 1997 年就已经发表在 JBC 上了。因此，忠言逆耳：别指望 Genbank或任何一个数据库能够百分百做到数据无误且实时更新。 Features 里的注释内容与原核生物的数据库记录相似，CDS 指出了从 63 到 821 是一段编码区，在这段编码区里基因是连续的，因为是经过剪切后的成熟 mRNA，它将被翻译成线粒体型 dUTPase 蛋白。下面/translation 里给出的是计算机翻译出的该蛋白的序列。 在 Features 里还有两个新的条目之前没有见到过。sig_peptide 和 mat_peptide。sig_peptide，也就是 signal peptide，指出了编码信号肽的碱基的位置。信号肽决定了蛋白质的亚细胞定位，也就是蛋白质工作的地方。mat_peptide，也就是 mature peptide，指出了编码成熟肽链的碱基的位置。 基因组里的 DNA 序列，是非线性分布的基因序列。我们仍然浏览编码人的 dUTPase 的dut 基因序列。输入检索号 AF018430，搜索！

（三）基因组数据库 Ensemble

人的基因组有 33 亿个碱基分布在23 个染色体上。 从 Ensembl 数据库（http://www.ensembl.org）查看人的基因组。Ensembl 是由欧洲生物信息学研究所 EBI 和英国桑格研究院合作开发的。它收入了各种动物的基因组，特别是那些离我们人类近的脊椎动物的基因组。这些基因组的注释都是通过配套开发的软件自动添加的。Ensembl 主页左下角有人，老鼠，斑马鱼这三个点击率最高的基因组的快速链接。 点击Human进入之后，我们点这个查看染色体，就可以看到人的所有染色体的图例 前面一直研究的那个编码 dUPTase 的 dut 基因就在 15 号染色体上。点一下 15 号染色体，在弹出窗口中选择染色体概要（chromosome summary）。这时我们会得到 15 号染色体的一个一览图。里面包括编码蛋白的基因、非编码基因、假基因分别在染色体上不同区段内的含量，以及 GC百分比（红线），和卫星 DNA 百分比（黑线）。染色体统计表给出了 15 号染色体的长度，以及各种类型的基因的个数。 从 Genbank 我们了解到，dut 基因的第三号外显子位于 15 号染色体的长臂条带 21.1 附近。所以我们进一步进入这个条带看一下。点击条带 21.1，选择区间链接。这时，这个区间内所有的基因就都被显示在一张图上。 在以 dut 基因为中心显示的放大图谱中，点击 dut 或者对应的区域，在弹出的概况窗口中选择 Ensemble 数据库的检索号。之后就会出现 dut 基因在 Ensemble 数据库中的详细记录。

（四）微生物宏基因组数据库 JCVI

微生物宏基因组数据库是非常有用的一级核酸数据库资源。说到微生物宏基因组学，不得不介绍的是美国基因组研究所 TIGR 和克莱格反特研究所 JCVI。美国基因组研究所致力于微生物基因组的研究，也有部分植物基因组项目。TIGR 是NCBI 基因组资源的有力补充，因为它不仅拥有已完成测序的基因组，还有那些测序中的基因组信息。在植物基因组项目中可以找到拟南芥、玉米、苜蓿和柳树的基因组信息。在微生物与环境基因组目中，特别值得关注的是“人类微生物组计划”，HMP。 从 JCVI 主页（http://www.jcvi.org/）的统计表中我们可以看到不同器官中有多少微生物基因组已被测序并被注释。点击下方的统计链接。可以得到 HMP 中已研究的所有微生物基因组。这些微生物在人体中存在的位置，测序及注释是已完成还是在进行中。已完成的基因组后面会有三个链接： WGS 是全基因组鸟枪法测序项目数据库记录的链接。SRA 是高通量测序数据库记录的链接。这两个链接里记录的是测序的信息。相比之下对大家更为有意义的是ANNOTATION 链接里的内容，他列出了某个基因组在 Genbank 中所有注释的链接。比如微生物 Acinetobacter radioresistens SK82 的基因组共分成 82 条序列记录在 Genbank 数据库中。  

三、 二级核酸数据库

二级核酸数据库包括的内容非常多。其中 NCBI 下属的三个数据库经常会用到。他们是RefSeq 数据库，dbEST 数据库和 Gene 数据库。RefSeq 数据库，也叫参考序列数据库，是通过自动及人工精选出的非冗余数据库，包括基因组序列、转录序列和蛋白质序列。凡是叫ref 什么的数据库都是非冗余数据库，就是已经帮你把重复的内容去除掉了。dbEST 数据库，也就是表达序列标签数据库，存储的是不同物种的表达序列标签。Gene 数据库以基因为记录对象为用户提供基因序列注释和检索服务，收录了来自 5300 多个物种的 430 万条基因记录。 此外，非编码 RNA 数据库，提供非编码 RNA 的序列和功能信息。非编码 RNA 不编码蛋白质但在细胞中起调节作用。目前该数据库包含来源于 99 种细菌，古细菌和真核生物的3 万多条序列。microRNA 数据库主要存放已发表的 microRNA 序列和注释。这个数据库可以分析 microRNA 在基因组中的定位和挖掘 microRNA 序列间的关系。

RefSeq	https://www.ncbi.nlm.nih.gov/refseq/
dbEST	https://www.ncbi.nlm.nih.gov/dbEST/
Gene	https://www.ncbi.nlm.nih.gov/gene
ncRNA	http://biobases.ibch.poznan.pl/ncRNA
microRNA	http://www.mirbase.org/

四、一级蛋白质序列数据库

蛋白质数据库的种类比核酸数据库要多，但它的注释要比核酸数据库直白得多。像核酸数据库一样，蛋白质数据库也分为一级和二级。一级蛋白质数据库又分为蛋白质序列数据库和蛋白质结构数据库。这两种数据库里存放的都是通过实验方法直接获得的基础数据。而二级蛋白质数据库都是在一级数据库的基础上分析加工出来的。一级蛋白质序列数据库包含三大蛋白质序列数据库，Swiss-Prot、TrEMBL 和 PIR，这三个数据库共同构成UniProt数据库。 Swiss-Prot 是一个人工注释的蛋白质序列数据库。它拥有注释可信度高，冗余度小的优点。它是由欧洲生物信息学研究所 EBI 与瑞士生物信息学研究所 SIB 共同管理的。 TrEMBL也是 EBI 和 SIB 共同管理的一个数据库，他与 Swiss-Prot 的区别是：TrEMBL 里的蛋白质序列注释是由计算机完成的，它包含了 EMBL 核酸序列数据库中为蛋白质编码的核酸序列的所有翻译产物。换言之，TrEMBL 是通过计算机，把核酸序列数据库里能编码蛋白的核酸序列都翻译成了蛋白质序列，然后把这些计算机翻译出来的蛋白质序列存入其中。TrEMBL 把已经包含在 Swiss-Prot 数据库中的序列剔除掉了。也就是在 Swiss-Prot 里已经有人工注释的蛋白质序列在 TrEMBL 里就不再出现了。 PIR 数据库是蛋白质信息资源数据库，他设在美国 Georgetown 大学医学中心。是一个支持基因组学，蛋白质组学和系统生物学研究的综合公共生物信息学资源。

（一）UniProt 数据库介绍

UniProt 数据库有三个层次。 第一层叫 UniParc，收录了所有 UniProt 数据库子库中的蛋白质序列，量大，粗糙。 第二层是 UniRef，他归纳了 UniProt 几个主要数据库并且是将重复序列去除后的数据库。 第三层是 UniProtKB，他有详细注释并与其他数据库有链接，分为 UniProtKB 下的 Swiss-Prot和 UniProtKB 下的 TrEMBL 数据库。 关系稍有点复杂，但实际上我们最常用的就是 UniProtKB下的 Swiss-Prot 数据库。   从 UniProt 数据库查看一条蛋白质序列（http://www.uniprot.org/）。在UniProt数据库的首页上有一个关于 UniProtKB 数据库的统计表。可以看到，TrEMBL 数据库里存储的序列数量远远大于 Swiss-Prot 中的。统计表里清楚的写着：TrEMBL 是自动注释的，没有经过检查，而 Swiss-Prot 是人工注释的，并且经过检查。 UniProt 数据库的首页上也有一个搜索条，选择UniprotKB 数据库，然后输入“human dutpase”，第一条就是我们要的。Entry 这一列是蛋白质序列在 UniProtKB 数据库中的检索号，Entry_Name 是检索名，检索号与检索名平行运行，都是一条序列在数据库中的唯一标识，两者作用相同，只是写法不同。从检索名可以更直观的知道是哪个物种的什么蛋白质。从加星文档图标（Entry_Name后一列）我们可以获知序列是被人工检查过的还是没有。也就是说，有加星文档图标的是 Swiss-Prot 中的数据，没有的是 TrEMBL 里的。后面这几列，依次是蛋白质的名字，编码这一蛋白质的基因的名字，所属物种以及序列长度。点击第一条序列的检索号，打开这条数据库记录。 UniProtKB 中的数据库记录分成几个部分，左侧是注释标签，点击其中某一个标签可以直接跳转到该部分注释。上方是工具标签，可以用于和其他序列进行比较，格式转换，存储等。工具标签下方是这条蛋白质序列的基本信息，蛋白质的名字，基因的名字，所属物种，以及状态。这里有加星文档图标，是被人工检查过的，应该属于 Swiss-Prot 数据库。注释打分 5 星，说明注释得很全面，并且这些注释在蛋白质水平上有实验依据。 Function：功能这部分注释很详细的说明了这个蛋白质的功能。从这里可以得知dUTPase 是一种在核酸代谢过程中的酶、它的催化反应方程式、它的辅助因子、它参与的代谢途径等。每条注释信息都提供出处来源，让你有据可查。 Names & Taxomomy： 给出了蛋白质的各种名字，包括全称、缩写以及别名。还列出了所属物种以及该物种的分类学谱系等。 Subcellular location：提供蛋白质亚细胞定位（subcellular localization）的信息。目前，研究亚细胞定位的数据来源基本都是Swiss-Prot 数据库。 Pathology & Biotechnology：提供蛋白质突变或缺失导致的疾病及表型信息。 PTM/Processing：提供蛋白质翻译后修饰或翻译后加工的相关信息。 Expression：提供了基因在 mRNA 水平上的表达信息，或者在细胞中蛋白质水平上的表达信息，或者在不同器官组织中的表达信息。 Interaction：提供了蛋白质之间相互作用的信息。包括 UniProtKB 中直接与这个蛋白质有两两相互作用的蛋白质序列的链接，以及这个蛋白质在各种蛋白质相互作用数据库或蛋白质网络数据库中涉及的数据库记录链接。 Structure：提供蛋白质二级结构和三级结构信息。只有那些已通过实验方法测定三级结构并且已提交到蛋白质结构数据库 PDB 的蛋白质才有结构注释。二级结构以图形拓扑的形式呈现。三级结构列出了该蛋白质在蛋白质结构数据库 PDB 中涉及的数据库记录链接。这些结构经常只对应蛋白质的部分序列。 Family & Domains：提供蛋白质家族及结构域信息。 Sequence：提供蛋白质氨基酸序列信息。含有多个异构体的蛋白质会显示多条序列。 Cross-references：列出了所有通往其他含有该蛋白质信息的数据库的链接。 Publications：列出了有关这个蛋白质已发表的所有文献的信息。 Entry information：提供有关这条数据库记录的录入信息，外加一个免责声明。 Miscellaneous：杂项，包含任何无法归入前几项的内容。 Similar Proteins：在 UniRef 数据库里找到与该蛋白质在序列水平上相似的其他蛋白质，并按相似度高低分组。

五、一级蛋白质结构数据库

蛋白质结构数据库 PDB（http://www.rcsb.org ）是全世界唯一存储生物大分子3D 结构的数据库，这些生物大分子除了蛋白质以外还包括核酸以及核酸和蛋白质的复合物，只有通过实验方法获得的 3D 结构才会被收入其中。 在 PDB 网站的搜索条中输入 “3H6X”，也就是把作者的名字和蛋白质的名字同时输入搜索条，然后点 go。 PDB 数据库的检索号，俗称 PDB ID，是由字母和数字组成的四位编号。一个结构对应一个 PDB ID，而不是一个蛋白质对应一个 PDB ID，因为同一个蛋白质在 PDB 数据库中可以有很多个结构。他们可以是不同作者提交的，也可以是一个蛋白的不同结构形态。网页上的信息都是关于这个结构的基本描述以及解析结构所用的实验参数。真正的结构信息要从Download files 里面下载。结构信息存储在 PDB 格式的一个纯文本文件里，这种文件叫做 PDB文件。PDB 文件都是以 PDB ID 命名，以“.PDB”为后缀，可以用记事本打开。 PDB 文件和我们GenBank 以及 UniProtKB 的纯文本数据库记录差不多。也是每行有条目索引词，后面是具体内容。我们通过浏览 3H6X 的 PDB 文件，看看这样的文本记录如何呈现 3D 结构。 第一部分：头信息 HEADER：蛋白质结构的基本信息描述，包括分子类别，存储日期，PDB ID TITLE：结构的标题 COMPND：对结构中各个分子的描述。从这里可以看出 3H6X 这个结构是由三条链形成的三聚体结构。 SOURCE：结构中所包括的每一个分子的实验来源。 KEYWDS：用于数据库搜索的关键词 EXPDTA：测定结构所采用的实验方法。PDB 中绝大部分结构都是通过 X 射线衍射法测定的，少数是核磁共振法，极少数是使用包括电子显微镜在内的其他方法测定的。 AUTHOR：作者信息 REVDAT：历史上曾经对该数据库记录进行过的修改。 JRNL：发表结构的文献信息。 REMARK：无法归入其他部分的注释。 第二部分：一级结构信息（也就是氨基酸序列） DBREF：该蛋白质在蛋白质序列数据库里的检索号等信息。 SEQRES：氨基酸序列。 MODRES：对标准残基上的修饰，比如第 56 号位置的蛋氨酸被硒代蛋氨酸所取代。 第三部分：非标准残基信息 HET：非标准残基及位置。 HETNAM：非标准残基的化学名称。 FORMUL：非标准残基的化学式。 第四部分：二级结构信息 HELIX：位于螺旋结构上的氨基酸所在位置及所属链。 SHEET：位于折片结构上的氨基酸所在位置及所属链。 TURN：位于转角结构上的氨基酸所在位置及所属链。 Link：残基间的化学键。比如 106 号氨基酸上的 C 与 107 号氨基酸上的 N 之间的化学 键是肽键！键长 1.32 埃。除了肽键还可能有氢键，二硫键等等。 第五部分：实验参数信息 CRYST1：晶胞参数。 ORIGXn：直角-PDB 坐标。 SCALE：直角部分结晶学坐标。 第六部分：3D 坐标信息 ATOM：PDB 文件中最重要的，也是篇幅最长的就是 3D 坐标部分。每一行是一个原子。包括原子号，原子名，这个原子所在氨基酸的名字，属于哪条分子链以及所在氨基酸的编号。后面这三个数就是这个原子在三维空间里的坐标，X 轴 Y 轴和 Z 轴。通过这个 3D 坐标，蛋白质的每一个氨基酸上的每一个原子都能找到自己的空间位置。所有原子按照各自的空间位置站好，就构成了整个蛋白质的空间结构。至此，我们终于知道了 PDB 是如何存储 3D 结构了。它存储的实际上是原子的 3D 坐标。 CONECT：原子间化学键连接信息。 MASTER：版权拥有者信息。 END*：结束符。 PDB 数据库就提供这样一个在线的可视化软件，叫 JSmol。JSmol 基于 JAVA 开发，所以需要先安装 java 运行环境。JAVA 可以到 JAVA官网下载。安装好 JAVA 后，重启浏览器，打开这个网页，点 JSmol 链接，之后，接受 java，信任 java，运行 java。（如果 IE 浏览器打开 JSmol 有问题，可以尝试 360 等其他浏览器）。网页加载完成之后，页面上会出现一个图片。这个图片貌似是当前这个蛋白质的结构图。   按照 JSMOL 的操作规则，按住鼠标左键拖拽是旋转结构，鼠标中键可以放大缩小，右键可以打开 JSmol 菜单，进行更多操作。通过设置右边的参数，可以改变 3D 结构的显示方式。在线版本的 JSmol 可以方便快捷的查看结构，但是功能不够全，特别是缺少分析功能。在后续的文章中，我将为你详细讲解功能更为强大的 3D 可视化软件的使用，并且在这一章里还将涉及 PDB 数据库的更多内容。

六、二级蛋白质结构数据库

（一）结构域家族数据库 Pfam

Pfam 数据库（http://pfam.xfam.org/）是一个蛋白质结构域家族的集合，包括了一万六千多个蛋白质家族。以Toll 样受体蛋白为例解读搜索结果。Pfam 主页上的搜索工具可以帮助我们查找某条序列上有哪些结构域。  

这是一条 Toll 样受体蛋白的序列，序列如下：

MMSASRLAGTLIPAMAFLSCVRPESWEPCVEVVPNITYQCMELNFYKIPDNLPFSTKNLDLSFNPLRHLGSYSFFSFPELQVLDLSRCEIQTIEDGAYQSLSHLSTLILTGNPIQSLALGAFSGLSSLQKLVAVETNLASLENFPIGHLKTLKELNVAHNLIQSFKLPEYFSNLTNLEHLDLSSNKIQSIYCTDLRVLHQMPLLNLSLDLSLNPMNFIQPGAFKEIRLHKLTLRNNFDSLNVMKTCIQGLAGLEVHRLVLGEFRNEGNLEKFDKSALEGLCNLTIEEFRLAYLDYYLDDIIDLFNCLTNVSSFSLVSVTIERVKDFSYNFGWQHLELVNCKFGQFPTLKLKSLKRLTFTSNKGGNAFSEVDLPSLEFLDLSRNGLSFKGCCSQSDFGTTSLKYLDLSFNGVITMSSNFLGLEQLEHLDFQHSNLKQMSEFSVFLSLRNLIYLDISHTHTRVAFNGIFNGLSSLEVLKMAGNSFQENFLPDIFTELRNLTFLDLSQCQLEQLSPTAFNSLSSLQVLNMSHNNFFSLDTFPYKCLNSLQVLDYSLNHIMTSKKQELQHFPSSLAFLNLTQNDFACTCEHQSFLQWIKDQRQLLVEVERMECATPSDKQGMPVLSLNITCQMNKTIIGVSVLSVLVVSVVAVLVYKFYFHLMLLAGCIKYGRGENIYDAFVIYSSQDEDWVRNELVKNLEEGVPPFQLCLHYRDFIPGVAIAANIIHEGFHKSRKVIVVVSQHFIQSRWCIFEYEIAQTWQFLSSRAGIIFIVLQKVEKTLLRQQVELYRLLSRNTYLEWEDSVLGRHIFWRRLRKALLDGKSWNPEGTVGTGCNWQEATSI

搜索结果显示，一共找到 4 个区域匹配 Pfam 数据库中已记录的结构域。前三个是 Toll 样受体蛋白胞外域典型的重复序列片段。最后一个是 TIR 结构域，以TIR 结构域为例查看详细信息如下： Summary 里可以获得这个结构域的功能注释以及结构信息。 Domain Organization 里可以看到目前有多少蛋白质拥有 TIR 结构域，以及 TIR 结构域和其他结构域之间的组合搭配关系。 Structure 会列出目前所有包含 TIR 结构域的蛋白质结构，以及他们在序列数据库UniProt 和结构数据库 PDB 中的链接。同时，也提供 JSmol 在线结构查看工具。

（二）结构分类数据库 CATH

根据结构域的空间特征可以对结构域进行分类。CATH 和 SCOP 是两个重要的蛋白质结构分类数据库。CATH 数据库（http://www.cathdb.info/）由伦敦大学 1993 年创建。CATH这个数据库的名字 C、A、T、H 是数据库中四种结构分类层次的首字母。也就是，所有蛋白质结构域在 CATH 中被首先分成 4 种 CLASS，这就是 C。四种 CLASS 分别是全α型，全β型， α +β型，低二级结构型。每一个 Class 中的结构域又被具体分为不同的 architecture，也就是 A。A 这一层是按照螺旋和折叠所形成的超二级结构排列方式分类的。比如α +β这个 class 下的结构可以进一步分为桶状的，三明治状的，还有滚轴状等 Architecture。每种 Architecture 里的结构域，又可以根据二级结构的形状和二级结构间的联系更进一步分为不同的 topology，也就是T。最后再通过序列比较以及结构比较确定同源性分类，划分出不同的 homologous superfamily，也就是 H。这样每个结构从粗到细，即从 A 到 H，会有四个层次的分类。注意结构分类是以结构域为单位进行的，而不是针对整个蛋白。所以 PDB 中的一个蛋白质结构可能对应 CATH中多个结构域分类。CATH 在分类时既使用计算机程序，也进行人工检查。 此外，CATH-Gene3D 还为超过 500 万条来自公共数据库的蛋白质序列进行了结构分类预测。Gene3D 里的信息为绝大多数还未解析 3D 结构的蛋白质提供了重要的功能研究依据。 搜索条输入 3H6X，这是之前在 PDB 数据库里查看过的 dUTPase的结构。结果显示 dUTPase 蛋白的结构分类代码是 2.70.40.10。CATH 为每一层的每一种结构分类命名，并用数字代号代表这一分类。因此每个结构域会具有一个分类代码。第一个数字是 C 这一层的分类代码，第 2 个数字是 A 这一层的分类代码，第 3 个数字是 T 这一层的分类代码，第 4 个数字是 H 这一层的分类代码。    这里，CATH 把所有拥有 2.70.40.10结构分类的结构域，根据他们的序列相似度不同，进行了聚类。不同深浅的圈代表不同的序列相似度。通过这张图，我们可以了解到具有相同结构分类的蛋白质他们在序列水平上的亲缘关系远近（下图左）。   此外，CATH 还从 2.70.40.10 这个结构分类里挑出了 19 个有代表性的结构域，并且把他们的 3D 结构叠加在了一起。从这个图上（上图右），我们可以看到这个结构分类的总体特征以及差异产生的位置。

（三）结构分类数据库 SCOP2

SCOP 数据库与 CATH 类似，也属于蛋白质结构分类数据库，但 SCOP 的分类原则更多考虑蛋白质间的进化关系，而且分类主要依赖于人工验证。和 CATH 一样，SCOP 的结构分类也基于四个层次。第一层也叫 Class，也是基于二级结构成分分类。Class 之下是 Fold，主要考虑结构的空间几何关系。再往下是 Superfamily，基于远源的蛋白质进化关系分类。最后是 Family，基于近源的蛋白质进化关系分类。注意 SCOP 和 CATH 里面都有提到 Superfamily这个词，但两者的含义并不相同。CATH 里 Superfamily 是指的从 C 到 A 到 T 再到 H 这样四层的一个精细结构分类。而 Scop 中，Superfamily 是结构分类的第三个层次的名称。目前，SCOP 已升级为 SCOP2（http://scop2.mrc-lmb.cam.ac.uk） SCOP2 的主页上也有搜索条，可以查看某一个 PDB 结构的结构分类。   搜索结果中的第 2 到第 5 条，就是该蛋白质结构的四层分类。第一层 Class，第二层 Fold，第三层 Superfamily，第四层 Family。第一层 Class 之上是 SCOP 数据库的根。第 4 层 family 之下是这个蛋白质的名字，再往下是所属物种。虽然从这个谱系上看有 7 个层次，但实际上真正的结构分类只有中间四层。

七、专项数据库

（一）京东基因与基因组百科全书 KEGG

KEGG，全称京东基因与基因组百科全书（http://www.genome.jp/kegg）。它是关于基因、蛋白质、生化反应以及通路的综合生物信息数据库。由多个子库构成。 这些子库中，KEGGPATHWAY 数据库包含了大量物种的代谢与生物信号传导通路信息。Pathway 数据库下又分为 7 个部分：1）Metabolism，2）Genetic Information Processing，3）Environmental Information Processing，4）Cellular Processes，5）Organismal Systems，6）Human Diseases，7）Drug Development。其中 Metabolism 代谢通路这部分，又具体分为几个专题：1）Global/overview，2）Carbohydrate，3）Energy，4）Lipid，5）Nucleotide，6）Amino acid，7）Other amino，8）Glycan，9）Cofactor/vitamin，10）Terpenoid/PK，11）Other secondary metabolite，12）Xenobiotics，13）Chemical structure。 选择KEGG PATHEAY，然后选择Metabolism下的Global/overview，最后选择01100这个选项   这是会看到一个像电路板的图，其中有一个圈这个圈的名字是 TCA 循环，也就是三羧酸循环。这个图上的每一个圆点儿代表一个化合物，把鼠标放在某一个点上，会出现化合物的分子式，点击可进入相应数据库查看详细。   图上的每一条线代表一个生化反应。把鼠标放到三羧酸循环的名字上，可以看到一个更加详细的通路图，我们点击这个名字。得到三羧酸循环详细的通路图，其中圆圈是化合物，箭头代表反应及反映方向。方块中的是酶。虚箭头指向其他途径，中间过程省略没有列出。当点击某一个酶的时候会直接进入 KO 数据库。KO 是 KEGG 中的一个“专有名词”，表示蛋白质或者说酶的一个分类体系。序列高度相似，并且在同一条通路上有相似功能的蛋白质被归为一组，然后打上一个 KO 标签。从对应的 KO 数据库记录中可以查看当前这个酶。 接下来我们看一下 Toll 样受体的信号传导通路。它位于 KEGG Pathway 数据库里的Organismal Systems 部分的 Immune system 专题里。人的各种 Toll 样受体信号传导通路图如下。 可以看到在细胞膜和内质网上，有很多种 Toll 样受体。它们识别了不同的入侵物后，激活下游蛋白，一个接一个的传递信号，直至产生各种细胞因子，激发炎性反应。如果点击其中的 Toll 样受体 4，可以看到这个蛋白质的详细信息，包括它参与的各种 Pathway。此外，还有该蛋白可引发的疾病信息，并且可以链接到 KEGG 人类疾病数据库。此外，数据库记录里还提供了两个相关的 Drug，其中一个 Drug 叫 Eritoran。它是 Toll 样受体 4 的拮抗剂。因为它长得和 Toll 样受体 4 的激动剂 LPS 脂多糖很像，所以可以被 Toll 样受体 4 捕获。但是因为它又比激动剂 LPS 少了两条链，所以 Toll 样受体 4 捕获它之后不能激活下游的信号传导，从而使 Toll 样受体 4 丧失免疫功能。这种药可以用于 Toll 样受体 4 引发的自身免疫疾病的治疗。

（二）人类孟德尔遗传在线 OMIM

一个有关人类遗传病的数据库，人类孟德尔遗传。它是一个将遗传病分类并链接到相关人类基因组中的数据库。它的在线版本是人类孟德尔遗传在线 OMIM。OMIM为临床医生和科研人员提供了权威可信的关于遗传疾病及相关疾病基因位点的详细信息。从NCBI 的OMIM子库页面点击Getting Started进入数据库（http://www.ncbi.nlm.nih.gov/omim），或者直接从OMIM主页进入（http://www.omim.org/）。 以阿尔茨海默症（AD）为例，在搜索条中输入：alzheimer disease。搜索结果里排在第一位的就是我们想要的。点击进入后数据库给出了与 AD 相关的致病基因，如下图。 包括他们在染色体中的位置，所引发表型的数据库编号，以及基因的数据库编号等。此外，页面上还提供大量的文字信息。如果我们点击某一个染色体定位的话。会出现这个位置附近基因的列表，以及引发的各种疾病。点击某一基因的数据库编号，可以查看这个基因的详细信息。

第二章 序列比较

生物序列有自己的书写格式，而且格式很多。不同的处理软件会用到不同的格式，但是最常用的，大多数软件都识别的格式是 FASTA 格式。我们用一致度（identity）和相似度（similarity）这两个指标来定量描述序列有多相似。

一、 替换记分矩阵

替换记分矩阵是反映残基之间相互替换率的矩阵。也就是说，它描述了残基两两相似的量化关系。比如图2.1 就是一个替换记分矩阵。矩阵中行和列分别是 20 种氨基酸，且两两之间有一个分值。根据这个分值就可以知道谁和谁相似，谁和谁不相似。替换记分矩阵有很多种。DNA 序列有 DNA 序列的替换记分矩阵，蛋白质序列有蛋白质序列的替换记分矩阵，两者不可混用。 图 2.1. BLOSUM-62 替换记分矩阵

（一）核酸序列的替换记分矩阵

DNA 序列的替换记分矩阵主要有三种。一个是等价矩阵（图2.2-1）。这个矩阵最简单.其中，相同核苷酸之间的匹配得分为 1，不同核苷酸间的替换得分为 0。由于不含有碱基的理化信息和不区别对待不同的替换，在实际的序列比较中很少使用，一般只用于理论计算。 另一种是转换-颠换矩阵（图2.2-2）。转换即嘌呤变嘌呤，嘧啶变嘧啶；颠换即嘌呤变嘧啶，或者嘧啶变嘌呤。大自然更倾向于接受嘌呤和嘌呤之间的替换，以及嘧啶和嘧啶之间的替换，而嘌呤和嘧啶之间的替换会导致不好的事情发生，这种替换大多在进化过程中已经被淘汰了。为了反映这一情况，转换-颠换矩阵中，转换的得分比颠换要高为-1 分，而颠换的得分为-5分。 还有一种叫 BLAST 矩阵。经过大量实际比对发现，如果令被比对的两个核苷酸相同时得分为+5 分，不相同为-4分，这时比对效果最好。这个矩阵广泛地被 DNA 序列比较所采用。 图2.2 DNA序列的三种替换记分矩阵

（二）蛋白质序列的替换记分矩阵

蛋白质最常用的两种矩阵是PAM矩阵和BLOSUM矩阵。PAM矩阵基于进化原理。如果两种氨基酸替换频繁，说明自然界容易接受这种替换，那么这一对氨基酸替换的得分就应该高。PAM 矩阵是目前蛋白质序列比较中最广泛使用的记分方法之一。基础的 PAM-1矩阵反应的是进化产生的每一百个氨基酸平均发生一个突变的量值，由统计方法得到。PAM-1自乘n次，可以得到PAM-n ，表示发生了更多次突变。我们需要根据要比较的序列之间的亲缘关系远近，来选择适合的 PAM 矩阵。如果序列亲缘关系远，也就是说序列间会有很多突变，那就选 PAM 后面跟一个大数字的矩阵。如果亲缘关系近，也就是突变比较少，序列间大多数地方都是一样的，那就选 PAM 后面跟一个小数字的矩阵。 图2.3是 PAM250矩阵。对角线上的数值为匹配氨基酸的得分。其他位置上≥0 的得分代表对应的一对氨基酸为相似氨基酸，<0的是不相似的氨基酸。 图2.3 PAM250 矩阵 BLOSUM 矩阵有和 PAM 矩阵相同的地方，也有不同的地方。相同的是，BLOSUM 矩阵后面也带有一个编号，有很多种 BLOSUM 矩阵。不同的是，BLOSUM 矩阵都是通过对大量符合特定要求的序列计算而来的。这点和 PAM 矩阵不同的。PAM-1 矩阵是基于相似度大于85%的序列计算产生的，也就是通过关系较近的序列计算出来的。那些进化距离较远的矩阵，如 PAM-250，是通过 PAM-1 自乘得到的。也就是说，BLOSUM 矩阵的相似性是根据真实数据产生的，而 PAM 矩阵是通过矩阵自乘外推而来的。和 PAM 矩阵的另一个不同之处是BLOSUM 矩阵的编号。这些编号，比如 BLOSUM80 中的 80，代表这个矩阵是由一致度≥80%的序列计算而来的。同理，BLOSUM62 是指这个矩阵是由一致度≥62%的序列计算而来的。因此，BLOSUM 后面跟一个小数字的矩阵适合用于比较相似度低的序列，也就是亲缘关系远的序列；而 BLOSUM 后面跟一个大数字的矩阵适合比较相似度高的序列，也就是亲缘关系近的序列。 图2.4是BLOSUM 62矩阵.样子和 PAM 矩阵差不多，但是里面的数值是不一样的。同样，≥0 的得分代表对应的一对氨基酸为相似氨基酸，<0 的是不相似的氨基酸。 图2.4 BLOSUM62 矩阵 总之，PAM1对应的氨基酸差异是1%，这是基础矩阵，由实际数据计算得出。而 PAM11 是由 PAM1 自乘 11 次得到的，他对应的氨基酸差异可不是 11%，而是大约在 10%左右（图2.5）。同样，PAM80 对应的差异也不是 80%，而是在 50%左右。如果你要比对的序列亲缘关系远，比如氨基酸差异在 80%左右，那就得选 PAM 自乘次数非常多的矩阵，适合的是 PAM246。但是现成的 PAM 矩阵也不是什么号的都有，只有几个关键号的。比如这个 PAM246 就没有，有的是 PAM250。 图2.5氨基酸差异与矩阵编号对照图   再来看BLOSUME。BLOSUME 后面的号和 PAM 刚好相反，因为它对应的是序列的相似度。差异在 80%左右意味着相似度在 20%左右，所以这个档次上的序列适合用的 BLOSUM 矩阵就是BLOSUM20。概括的说，PAM 后面的数体现的是序列的差异度，但不直接等于差异度，只是成对应关系而已；BLOSUM 后面的数体现是的序列的相似度并且直接等于相似度。所以我们看到，随着差异度的增大，适用的 PAM 矩阵后面的编号是增大的，而 BLOSUM 矩阵后面的编号是减小的（图2.5）。 进一步总结一下（图2.6），亲缘关系较近的序列之间的比较，用 PAM 数小的矩阵或BLOSUM 数大的矩阵；而亲缘关系较远的序列之间的比较，用 PAM 数大的矩阵或 BLOSUM数小的矩阵。对于关系较远的序列之间的比较，由于 PAM250 是通过矩阵自乘推算而来的，所以其准确度受到一定限制。相比之下BLOSUM 矩阵更具优势。对于关系较近的序列之间的比较，用 PAM 或 BLOSUM 矩阵做出的比对结果，差别不大。如果关于要比较的序列你不知道亲缘关系远近，那么就闭着眼睛用BLOSUM62 吧！ 图2.6 序列亲缘关系远近与矩阵的选择 除了 PAM 和 BLOSUM 矩阵，还有两个蛋白质的替换记分矩阵。一个是遗传密码矩阵（图2.7），它是通过计算一个氨基酸转换成另一个氨基酸所需的密码子变化的数目而得到的。矩阵的值对应为据此付出的代价。如果变化一个碱基就可以使一个氨基酸的密码子转换为另一个氨基酸的密码子，则这两个氨基酸的替换代价为 1；如果需要 2 个碱基的改变，则替换代价为 2；再比如从蛋氨酸（Met）到酪氨酸（Tyr）三个密码子都要变，则代价为 3。遗传密码矩阵常用于进化距离的计算，它的优点是计算结果可以直接用于绘制进化树，但是它在蛋白质序列比对，尤其是相似程度很低的蛋白质序列比对中，很少被使用。 另一个疏水矩阵（图2.8），它是根据氨基酸残基替换前后疏水性的变化而得到的矩阵。若一次氨基酸替换导致疏水特性不发生太大的变化，则这种替换得分高，否则替换得分低。疏水矩阵物理意义明确，有一定的理化性质依据，适用于偏重蛋白质功能方面的序列比对。在这个矩阵里，氨基酸按照亲疏水性排列。前边是亲水的，后面是疏水的。 有了替换记分矩阵，就可以知道哪些氨基酸是相似的。比如下面这个例子里： 序列1：C LH K 序列2：C I H L 我们可以从替换记分矩阵中读出I和L相似，K和L不相似。因此，它们的相似度就是 2个相同的加上1个相似的，除以长度4，等于75%。 这是两个序列具有相同长度，如果两个序列的长度不相同，怎么计算它们的一致度和相似度呢？后面再做介绍。 图2.7遗传密码矩阵 图2.8疏水矩阵  

二、 序列两两比较之打点法

（一）打点法的用途

比较两个序列的方法有打点法和序列比对法。打点法是最简单的比较两个序列的方法，理论上可以用纸和笔来完成。如果要比较下面这两条序列： Seq 1:THEFASTCAT Seq 2:THEFATCAT 我们需要把序列 1 整齐的水平书写，然后把序列 2 整齐的竖直书写，然后依次横横竖竖的比较每一个位置上的残基。相同的话就在这个位置上打个点，不同话，什么也不干（图2.9） 图2.9 序列 1 和序列 2 的打点图 从这样一个打点矩阵上可以总结出一些东西来。首先，说是打点你也可以打叉。这就是做个标记，标记下这个位置上对应的两个残基是相同的。第二点发现，这个矩阵中绝大部分地方是没有点的，只有少数位置上有点。第三个发现，这个矩阵中打点打出了一条较为明显的对角线。在打点矩阵中，连续的对角线及对角线的平行线代表两条序列中相同的区域。这个矩阵中在主对角线位置上连续的红色的对角线说明这个位置对应的序列 1 的部分和序列 2 的部分是完全相同的，都是 THEFA。此外，跟红对角线平行的蓝色平行线和绿色平行线，同样指出了序列1和序列2中两条相同的序列。也就是序列1和序列2中对应位置的 TCAT，以及序列 1 和序列 2 中对应位置的 AT。由这三条线，我们找到了序列 1 和序列 2 中三条相同的子序列。最后，我们放眼全局，红色的线和蓝色的线加起来基本上构成了一条主对角线。由此我们可以得出结论：序列 1 和序列 2 是比较相似的两条序列。事实上，如果直接看一下这两条序列，确实是挺相似的。如果是风马牛不相及的两条序列，做出的打点矩阵里是不会出现对角线的，哪怕是模糊的对角线，也不会出现。比如，让序列 1 和序列 3 打点做出的打点图（图2.10）中，完全是散点，根本就没有连续的线，更别提主对角线了。 图2.10 序列1和序列3的打点图 除了可以用打点法给两条不同的序列打点，还可以用一条序列自己跟自己打点。这样可以发现序列中重复的片段。比如我们让下面这条序列自己和自己打点： Seq4：THEFASTHESTHE 这样的打点矩阵必然是对称的，并且一定有一条主对角线（图2.11）。此外，在横向或纵向上，与主对角线平行的短平行线所对应的序列片段就是重复的部分。其中，红色短平行线对应的THE在序列中重复出现了3次。包括主对角线在内，平行线出现的次数就是重复的次数。 图2.11 序列自己和自己打点寻找重复序列 用这种方法我们还可以快捷的发现序列中的串联重复序列以及重复的次数。我们只要数数在半个矩阵中包括主对角线在内的所有等距的平行线的个数，就可以知道重复的次数，而且最短的平行线对应的序列就是重复单元（图2.12）。短的串联复序列具有高度多态性，也就是说不同的个体间重复次数存在差异，而且这种差异在基因遗传过程中一般遵循孟德尔共显性遗传规律，所以快速查找某些特定的短的串联复序列的重复次数可以用于法医学的个体识别或亲子鉴定等领域。 图2.12 打点法寻找串联重复序列

（二）Dotlet介绍

打点用笔和纸有点不现实，我们得用计算机软件来执行，可以自动打点的软件有很多（表2.1）。我们挑其中最常用的 Dotlet 软件做为演示（http://myhits.isb-sib.ch/cgi-bin/dotlet）。Dotlet 基于Java开发，所以电脑需要安装Java，如果没有安装，可自行去JAVA官网下载安装。 表2.1 打点法在线软件 网页正常打开后如图2.13所示。我们看到上边一行是参数设置选项，包括选择水平位置的序列，选择垂直位置的序列，选择替换记分矩阵，选择以多长的序列为单元打一个点，还有窗口显示的比例，以及计算按钮。左边是打点图的显示窗口，右边是调控打点图显示效果的区域，下面还有一块是序列信息显示区。网页打开时除了 input 按钮，其他位置都是灰的。 图2.13 Dotlet 在线打点工具界面 点击 input 按钮输入要打点的序列。在新打开的窗口里把要打点的序列拷贝进来。拷贝，黏贴。注意，序列输入窗口里只能输入纯序列部分。fasta格式里带大于号的第一行不能拷贝进来，否则Dotlet 会把这部分内容也当成序列进行打点。由于Dotlet不能自动识别FASTA格式里哪是名字，哪是序列，所以我们还需要手动的给这条序列取个名字，比如叫seq1，然后点ok，一条序列就输入进来了（图2.14）。 图2.14 打点序列的输入窗口   序列输入之后，参数设置行里水平序列的名字和垂直序列的名字都变成了seq1（图2.15）。 图2.15  Dotlet 参数选择 在这个例子里面我们让输入的这条 seq1 序列自己和自己打点。也就是水平位置的序列和垂直位置的序列都是 seq1。替换记分矩阵可以选讲过的各种 BLOSUM 矩阵或者PAM矩阵。如果选了这两种矩阵，表明并不是完全一致才打点，而是相似就打点。如果想要严格要求只有完全一致才打点的话，选Identity。再来看这个打点单元长度下拉条。这里如果选择1的话，就是一次只比较一个字母，这跟我们前面手动打点时的情况是一样的。如果选其他的，比如选择15，那就是一次比较15个字母，也就是看15个字母长度的序列整体的相似度如何来确定打不打点。注意这里不是比较完前15个字母，然后再从第16个字母开始比较后面的15个字母，而是第1次比较第1到第15个字母，然后再比较第2到第16个字母，再是第3到第17个字母，依次类推。所以选15和选1所进行的比较的次数是几乎一样多的，只是每次比较的序列单元长度不一样而已。显示比例窗口可以选择打点结果图是1:1原图显示，还是缩小为原图的 50%，或者25%显示。因为对于比较长的序列，这么小的窗口很难把原图显示完整。总结一下，第一个例子里水平序列seq1，垂直序列也是seq1，替换积分矩阵选Blosum62，比较单元长度15个字母，窗口显示比例1:1，最后点compute。 注意默认的颜色方案是在越相似的地方打的点的颜色越浅，越不相似的地方颜色越深。所以整体感觉像是在一张黑纸上打白点。因为是 seq1 自己和自己打点，所以应该有一条明显的主对角线，这是所有自己和自己打点的序列都会出现的情况。除此之外，还有一条与主对角线平行的次对角线，涉及 seq1大约1/2的长度。这说明 seq1 的前半部分和后半部分非常的相似（图2.16）。把鼠标点在这条次对角线上的任意位置，之后就可以从下面的序列显示区域看到，这一点，对应横着这条 seq1 的位置 153 开始的这一段，同时对应竖着这条seq1的位置11开始的这一段。为什么这一点对应的是一段呢？因为我们这里选的单元长度是15，也就是一次比较 15 个字母。从这里也可以看出，如果把 seq1 的前半段和后半段重合起来的话，它们是完全一样的序列。 图2.16 Seq1 自己和自己打点计算出的打点图 可以通过调整灰度条，来屏蔽大多数低分值的点，让他们统统变成黑色背景，并且强化高分值的点，让他们以纯白色突出显示出来（图2.17） 图2.17 调整灰度条改变打点图显示效果 上一个例子只是让大家熟悉一下 Dotlet 的使用界面及结果的查看方法，并没有实际的生物学意义。下面我们比较两条不同的序列，看看它们是否相似。将seq2 和 seq3 两条序列 input 进来。水平序列选 seq2，垂直序列选seq3，其他参数不变，点计算按钮。虽然这次的主对角线不如上一个例子里的明显，但是我们还是可以一眼就看出它来（图2.18）。说明这两条序列整体上十分相似。通过调整显示方案可以让主对角线清晰呈现。 图2.18  Seq2 和 seq3 打点计算出的打点图

三、双序列全局比对

双序列全局比对，经典的全局比对算法是 Needleman-Wunsch 算法。用 Needleman-Wunsch 算法为序列 p 和序列 q 创建全局比对。输入值除了两条序列之外，还要有替换积分矩阵以确定不同字母间的相似度得分，以及空位罚分（图2.19），空位罚分就是当字母对空位的时候应该得几分。我们还是希望一致或相似的字母尽可能的对在一起，字母对空位的情况和不相似的字母对在一起的情况一样，都不是我们希望的，还是少出现为好，所以通常字母对空位会得到一个负分，这个负分就叫做空位罚分。这里我们让空位罚分，也就是 gap 分值为-5 分。在比对中没有空位对空位的情况。输入值就是这些。 图2.19 全局比对输入值：序列 p 和序列 q，替换记分矩阵，空位罚分 接下来创建一个得分矩阵（图 2.20-A），并根据公式（图2.20-B）把得分矩阵填满。填满后全局比对就会跃然于纸上。得分矩阵的第一行是序列 p，第一列是序列 q，这一步和打点法很像。不过要注意，p和q的前面各留一个空列和一个空行，也就是第 0 列和第 0 行。 图 2.20 全局比对计算公式及得分矩阵 现在开始给得分矩阵赋值（图2.21）。根据公式：s(0,0)是初始值0。第0行：s(0,j) = gap * j，j 从 1 到 m，m 是序列 p 的长度。也就是 s(0,1)=gap1=-5，s(0,2)=gap2=-10，依次类推。第 0 行实际是一种极端情况的假设。也就是当序列 p 全部对空位时的得分。A 对空位是-5 分，AC 都对空位就累计到了-10 分，ACG 都对空位就累积到了-15 分，如果序列p 全部对空位，最终的累积得分就是-25 分。 第 0 列：s(i,0) = gap * i 第 0 列和第 0 行一样，也是反映了序列 q 如果全部对空位的累计得分。对一个空位累积gap1=-5 分，对两个空位累积 gap2=-10 分，对三个空位累积 gap3=-15 分，对四个空位累积 gap4=-20 分。 图2.21 得分矩阵中的第0行和第0列 第 0 行和第 0 列相对简单，其他的格就稍微复杂一点儿了。接下来填 s(1,1)（图2.22）。这个格里的值来源于三个值中的最大值。哪那三个值呢，一个是上面格 s(0,1)里的值加 gap，一个是左面格 s(1,0)里的值加 gap，还有一个是斜上格 s(0,0)里的值加当前这个位置字母对字母在替换记分矩阵里的分值 w(i,j)。什么意思呢？就是累积到这个位置时，是字母对字母得分高，还是序列 p 的字母对空位得分高，还是序列 q 的字母对空位得分高？有且只有这三种情况，我们要的是得分最高的那种情况。逐个看一下，上面格 s(0,1)+gap=-5+-5=-10。左面格 s(1,0) +gap=-5+-5=-10。斜上格 s(0,0)+w(1,1)=0+10=10。max(-10,-10,10)=10。所以当前这个格 s(1,1)的分值就是 10。此外，我们还需要用箭头记录一下这个 10 是从哪里来的。它是从斜上这个格来的，所以我们画一个指向斜上的箭头。 图2.22 得分矩阵中的 s(1,1) 接下这个格 s(1,2)值的计算（图2.23），仍然是找三个值中的最大值。上面格 s(0,2)+gap=-10+-5=-15。左面格 s(1,1)+gap=10+-5=5。斜上格 s(0,1)+w(1,2)=-5+-3=-8。max(-15,5,-8)=5。大值是5，来源于左面格 s(1,1)，画上向左的箭头。 图2.23 得分矩阵中的 s(1,2) 按照上面的公式，将整个得分矩阵填满。这时，我们再回过头来看一下第一行和第一列（图2.24）。其实，第一行的每一个值都是从左边的格加 gap 来的。所以我们给它们补上向左的箭头。第一列的每一个值都是从上边的格加 gap 来的。所以我们给它们补上向上的箭头。至此，所有的箭头和数值就都填好了。填满之后，右下角的分数就是整个全局比对最终的得分。然后从这个位置开始追溯箭头一直到左上角的零，并且把这些箭头标记出来。 图2.24 填满分值和箭头的得分矩阵 图2.24中标出的红色箭头是写出全局比对的唯一依据。追溯箭头是从右下角到左上角，但是写全局比对是从左上角开始，如果是斜箭头则是字符对字符，如果是水平箭头或垂直箭头则是字符对空位，箭头指着的序列为空位。我们看第一个是斜箭头，字母对字母，就是 A对 A，第二个是水平箭头，字母对空位，箭头指着的序列是空位，也就是C对空位。然后斜箭头G对A，斜箭头T对 T，斜箭头 C 对 C，一直写到右下角，全局比对就出现了（图2.25）。唯一的一个空位插在序列 q 的 A 与 A 之间，这样最终的比对得分最高。不信的话可以试试，其他任何一种插入空位的比对结果，得分都不会超过 21 分。因为我们在得分矩阵的创建过程中，每一步都是在上一步最优的情况下得出的当前最优结果。 图2.25 序列 p 和序列 q 的全局比对 四、双序列局部比对 一长一短的两条序列，比较局部比比较全长更有意义。局部比对的算法和全局比对很相似，只是在选最大值时通过增加了第四个元素“0”，来达到比对局部的效果。序列p和序列q，一长一短，其他输入值跟全局比对的一样（图2.26） 图2.26 局部比对输入值：序列p和序列 q替换记分矩阵，空位罚分 局部比对的计算公式在全局比对的基础上增加了第四个元素“0”。得分矩阵初始值仍是0，但第一行和第一列与全局比对不同，全是0（图2.27）。 图2.27 局部比对计算公式及得分矩阵 从 s(1,1)开始要选择四个值中的最大值。除了上面格 s(0,1)+gap=0+-5=-5，左边格 s(1,0)+gap=0+-5=-5，斜上格 s(0,0)+w(1,1)=0+-3=-3，还有一个 0。max(-5,-5,-3,0)=0。并且这个 0 既不是从上面格，也不是从左边格，以及斜上格三个方向来的，而是来自于公式里增加的“0”，所以不用画箭头（图2.28）。 s(1,4)的计算：上面格 s(0,4)+gap=0+-5=-5，左边格 s(1,3)+gap=4+-5=-1，斜上格 s(0,3)+w(1,4)=0+0=0，还有一个0。max(-5,-1,0,0)=0。这个0和s(1,1)的0是不一样的。这个0应该画上斜上的箭头（图2.29），因为它可以来自公式中的 0，也可以来自斜上格。而s(1,1)的0没有箭头因为它只来自公式中的0。两种情况虽然都是0，但来源不同，一定要通过箭头标识清楚。 图2.28 得分矩阵中的 s(1,1) 图2.29 得分矩阵中的 s(1,4) 按照公式，填充满整个得分矩阵（图2.30）。与全局比对不同，局部比对的得分不是在右下角，而是在整个矩阵中找最大值。这个最大值才是局部比对的最终得分，他可能出现在任何一个位置。这次箭头追溯也不是从右下角到左上角，而是从刚刚找到的最大值开始追溯到没有箭头为止。追溯箭头终止的位置也可以是得分矩阵中的任何一个位置。 图2.30填满分值和箭头的得分矩阵 最后根据标记好的箭头写出比对结果（图2.31）。从左上到右下标记的红色箭头依次是：斜箭头字母对字母，C 对 C；斜箭头字母对字母，G 对 G。相比这两条序列的全局比对结果，两端的空位在局部比对中就全部被忽略掉了。 图2.31 序列 p 和序列 q 的局部比对与全局比对的比较

五、一致度和相似度

两条长度不同的序列做全局比对，然后计算全局比对中一致字符的个数和相似字符的个数，再除以全局比对的长度，就可以得到它们的一致度和相似度了。比如下面这两条序列： 首先做出它们的全局比对，比对中一致字符的个数是 4 个，全局比对长度 6，一致度=67%。相似字符个数 1，相似度就是(4+1)/6=83%。 把长度相同的两个序列计算一致度和相似度的方法重新规范一下。尽管长度相同，但是做出的全局比对的长度并不一定等于序列的长度，比如下面这两条序列： 上下各加入一个空位，全局比对的长度就不等于序列的长度了。所以不管两条序列长度是否相同，都要先对它们做全局比对。让两条序列先以最优的方式比对起来，再从全局比对中数出一致字符和相似字符的个数，除以全局比对的长度，来得到它们的一致度和相似度。

六、在线双序列比对

（一）EMBL全局双序列比对工具

目前，使用率最高的是 EMBL 网站的双序列比对工具（http://www.ebi.ac.uk/Tools/psa）。打开页面，上面有全局比对工具、局部比对工具、还有基因组比对工具。首先看全局比对中的蛋白质序列比对工具（图2.31）。输入值非常简单，把要比较的两条蛋白质序列贴在输入框里或者上传。如果想要进一步设置比对的参数，可以点 More options。从这里可以选择使用哪种替换记分矩阵。按照之前讲过的原则，选择 PAM 矩阵或 BLOSUM 矩阵。如果实在不知道选哪个矩阵，就闭着眼睛选 BLOSUME62，下拉菜单里默认选的就是BLOSUM62。除了选择替换记分矩阵，这里还可以设置空位罚分，也就是gap的分值。这里实际上是让你选空位对字母的情况罚几分，所以显示的是正数，但在计算的过程中还是按照负数处理。这里的 gap 分好几种，一种叫“gap 开头（GAP OPEN）”，另一种叫“gap延长（GAP EXTEND）”，gap 开头就是连续的一串 gap 里面打头的那一个。gap 延长就是剩下的那些gap，这一串里，第一个gap 是gap开头，后面的都是 gap 延长。单独的一个gap按gap开头算。gap开头和gap延长可以分别定义不同的罚分。默认情况下，gap 开头罚分多，gap 延长罚分少（图2.32）。当gap开头小，gap延长大的时候，做出来的比对里面，gap很分散，极少有连续长串的gap出现（图2.33-A）。开头的一串gap是个例外，因为 seq2 太短，seq1的这一段只能跟gap相对。其他部分的gap都是分散出现的。这和我们默认参数（gap开头大，gap延长小）做出来的比对结果是截然不同的（图2.33-B）。在实际应用中，需要根据不同的情况选取不同的 gap 罚分，以满足不同的生物学意义。如果你对结果没有什么预期，那就请保持默认的参数。 除此之外，结尾的 gap 也可以划分出不同的种类并赋予不同的罚分，如果把 END GAPPENALTY 选成 true，就可以设置结尾的 gap 罚分了。结尾 gap 不太常用，特别是在做亲缘关系较近的序列比对时，是否设置结尾 gap，比对结果差别不大。 图2.31  EMBL全局双序列比对输入页面 图2.32  EMBL全局双序列比对输入页面 图2.33 设置不同的 GAP OPEN 和 GAP EXTEND

（二）EMBL局部双序列比对工具

EMBL的局部双序列比对工具可以选择经典的 Smith-Waterman 算法。More options 里面的参数设置和全局比对是一样的。在这个例子里，我们保持所有参数都为默认值，点提交（图2.34）。 图2.34  EMBL 局部双序列比对输入页面 从比对结果可以看出（图2.35），只有中间黑色的相似的部分出现在比对结果中了，两头红色的不相似的部分被忽略掉了。也就是只返回了局部最相似，得分最高的片段的比对结果。 如果给这两条序列做全局比对的话，会发现，绝大部分位置对得都很差，只有中间这一段对的还不错（图2.36）。所以，有时候两条序列并不同源，它们只是有一个功能相似的区域，这时用局部比对我们就能很快找到这一区域在两条序列中的位置。但是如果做全局比对的话，结果就不如局部比对明显了。 图2.35 局部双序列比对结果 图2.36 全局比对与局部比对的比较

（三）其他在线双序列比对工具

可以做双序列比对的工具很多（图2.37）。不同网站都有自己的比对工具，所使用的算法也不尽相同，但是它们的核心算法都是讲过的 Nidelmann-Wunsch 和 Smith-Waterman 算法，只是在他们的基础上有所变化，有所升级。 Biotools 的双序列比对工具无论是核酸序列还是蛋白质序列都能做。全局比对，局部比对，还包括蛋白质二级结构辅助的序列比对。功能比 EMBL 的只多不少。当然这都不是最关键的。关键的大招是除了能给出序列比对，还能给出得分矩阵。   图2.37 其他在线双序列比对工具  

七、BLAST比对

之前用EMBL的双序列比对工具做全局比对，虽然很快就出结果了，但至少也要经历一两秒钟的时间。而数据库中有几百万条序列，全部比对一遍，耗时太长。因此，我们需要快速的数据库相似性搜索工具。目前世界上广泛使用的就是 BLAST。它可以在尽可能准确的前提下，快速的从数据库中找到跟某一条序列相似的序列。BLAST 是 Basic Local Alignment Search Tool 的首字母缩写，直译过来就是基本局部比对搜索工具。BLAST 的基本原理很简单，要点是片段对的概念。所谓片段对是指两个给定序列中的一对子序列，它们的长度相等，且可以形成无空位的完全匹配。 BLAST 实际上是综合在一起的一组工具的统称，它不仅可用于直接对蛋白质序列数据库和核酸序列数据库进行搜索，而且可以将待搜索的核酸序列翻译成蛋白质序列后再进行搜索，或者反之，以提高搜索效率。因此BLAST可以分为BLASTp，BLASTn，BLASTx，tBLASTn和 tBLASTx（图2.38）。 BLASTp 也就是用蛋白质序列搜索蛋白质序列数据库， BLASTn 是用核酸序列搜索核酸序列数据库。 BLASTx 是将核酸序列按6条链翻译成蛋白质序列后搜索蛋白质序列数据库。 tBLASTn 是用蛋白质序列搜核酸序列数据库，核酸数据库中的核酸序列要按6条链翻译成蛋白质序列后再被搜索。 tBLASTx是将核酸序列按 6 条链翻译成蛋白质序列后搜索核酸序列数据库，核酸数据库中的所有核酸序列也要按 6 条链翻译成的蛋白质序列后再被搜索。 图2.38 各种 BLAST 示意图

（一）NCBI BLASTp

在 BLASTp 输入界面里（图2.39）：1）输入待搜索的蛋白质序列，这条序列可以在示例文件 blast.fasta 里面找到。2）指定搜索跟输入序列哪部分相似的序列，如果空着就是全长搜索。3）给搜索任务起一个名字，如果输入的是 FASTA 格式的序列，那么在输入框里面点一下，序列的名字就会被自动识别出来。4）如果在 Align two or more sequences 前面打勾的话，可以同时提交多个 BLAST 任务。 图2.39 NCBI BLASTp 输入界面上半部分 在输入界面的下部（图2.40）选择：1）被搜索的数据库。我们看到，虽然是 NCBI 的 BLAST工具，可以选择的数据库却不只 NCBI 下属的数据库，还包括其他组织的数据库，比如 PDB，Swissprot。事实上，各大数据库网站的 BLAST 工具都可以实现跨平台搜索。我们这次用 NCBI的 BLAST 工具搜索 SwissProt 数据库。2）Organism 可以把搜索范围限定在某一特定物种内，或者排除某一物种。3）在算法选择这一栏里，有之前提到的三种不同的 BLAST 算法，标准BLAST，PSI-BLAST 和 PHI-BLAST。这一次我们先尝试标准 BLAST。所有参数设置完毕之后，点 BLAST。 图2.40 NCBI BLASTp 输入界面下半部分 图2.41是搜索结果。最上面是第一部分搜索任务描述部分。输入界面里设置的各种参数都会在这里列出。第二部分（Graphic Summary）是图形化搜索结果部分。在图形化搜索结果里，BLAST工具识别出输入序列的第 25 到第 170 个氨基酸这一段属于 TIR 蛋白质家族。这部分里彩色线条构成的图告诉我们，一共从数据库中找到 50 个 hits，也就是高分匹配片段。注意这些线代表的是 50 个高分匹配片段而不是 50 条序列。一个高分匹配片段有可能是一条全长的序列，也就是全长匹配，也有可能只是某条序列的一部分，也就是局部匹配。代表这些高分匹配片段的线拥有不同颜色和不同的长短。如果把鼠标放到某一条线上，可以看到这条匹配片段的具体信息，包括他所在序列的数据库编号，序列的名字，匹配得分，期望值 E 值。匹配得分在 200 以上的用红线表示，80 到 200 之间的用粉线，50 到 80 的绿线，40-50 的蓝线，40 以下的黑线，所以颜色反映的是匹配的好坏程度。如果某一个高分匹配片段和输入序列是从头到尾匹配，就是全长的线，比如最上面的三条红线。如果只匹配输入序列的一部分，则是一条短线，短线所在的位置就是与输入序列匹配的位置。第三部分（Descriptions）是这 50 个高分匹配片段所在序列的详细信息列表。每条序列都有一个匹配得分和覆盖度。这两项决定了第二部分彩图中每条线的颜色和长短。除了匹配得分和覆盖度，表中还列出了其他指标。尤为重要的是 E-value。E-value 也叫做期望值或 E值。E 值越接近零，说 图2.41 NCBI BLASTp 搜索结果 明输入序列与当前这条序列为同一条序列的可能性越大。第三部分的表就是根据 E 值由低到高排序的。随着 E 值增大，匹配得分是成反比逐渐降低的。但是一致度与 E 值并非完全成反比。BLAST 没有做双序列比对，为了提高速度，它牺牲了一定的准确度。表中的一致度，是 BLAST 搜索完成后，针对搜索到的这 50 条序列专门做双序列比对而得到的。BLAST 牺牲掉的准确度对高度相似的序列，也就是亲缘关系近的序列构成不了威胁，不会把它们落掉，但是对于那些只有一点点相似，也就是远源的序列，就有点麻烦了，它们很有可能被丢掉而没有被BLAST发现。

（二）NCBI PSI-BLAST

为了提高速度，标准 BLAST 牺牲了一定的准确度，牺牲掉的准确度对高度相似的序列，也就是亲缘关系近的序列构成不了威胁，不会把它们落掉，但是对于那些只有一点点相似，也就是远源的序列，就有点麻烦了。它们很可能被落掉而没有被BLAST 发现。 要解决这个问题，可以用 PSI-BLAST。PSI-BLAST 可以搜罗出一个庞大的蛋白质家族，当然也包括标准 BLAST 不小心漏掉的那些远房亲戚。换言之，标准 BLAST 找到了直接认识的朋友，但朋友的朋友都丢掉了。PSI 是 Position-Specific Iterated 首字母缩写，中文是位点特异性迭代。PSI-BLAST 的特色是搜完一遍再搜一遍，且从第二次搜索开始，每次搜索前都利用上一次搜索到的结果创建一个位置特异权重矩阵以扩大本次搜索的范围。如此反复直至没有新的结果产生为止。位置特异权重矩阵（Position-Specific Scoring Matrix，简称 PSSM）是以矩阵的形式，统计一个多序列比对中，每个位置上不同残基出现的百分比。假设A的朋友只有B，B 朋友除A外还有C。如果输入序列的第一个位置是 A，那么在第一轮没有 PSSM辅助的情况下，只有第一个位置是 A 或 B 的序列被找到了。它们是图2.42-A 中所示的四条序列。根据这四条序列创建的 PSSM（图2.42-B）得知，第一个位置可以是 A，也可以是 B，那么在第二轮搜索中，除了 A 的朋友 B 之外，B 的朋友 C 也可以出现在第一个位置了。这样如此反复，我们就可以找到朋友的朋友了。 图2.42 位置特异权重矩阵 我们用 NCBI 的PSI-BLAST工具再次搜索图2.39中的序列，仍然在 SwissProt 数据库中搜索，算法选 PSI-BLAST（图2.43）。 图2.43 PSI-BLAST搜索 第一轮搜索结果（图2.44）和标准 BLAST 是一样的。略有不同的是，在第三部分（Descriptions）列表的最右侧多了两列。绿色标记的一列是将要为下一轮搜索创建 PSSM 所选择的序列，默认第一轮搜索到的序列都将用来创建 PSSM。红色标记的一列是已在本轮搜索中用来创建PSSM 的序列。因为是第一轮搜索，之前还没有搜索到任何序列，也就是第一轮的搜索过程中没有使用 PSSM，所以这一列都为空。接下来点“go”进行第二轮搜索。“go”左侧的输入框里可以指定列出搜索到的前多少条序列。因为 PSI-BLAST 无休止的一直做下去，也许会把数据库中的全部序列都找出来。 图2.44 NCBI PSI-BLAST 第一轮搜索结果 第二轮的搜索结果里（图2.45），红色标记的一列，也就是已在第二轮搜索中被用来创建PSSM 的序列，大部分都打勾了。但是，再往后面会看到有些标黄的序列没有打勾。这些没有打勾的标黄的序列就是在第二轮搜索中新找到的序列。它们将用于创建下一轮搜索使用的PSSM，但是在本轮搜索中，它们没有被用到，所以没有打勾。而没有标黄的这些打勾的序列，是在第一轮中就已经找到的序列。 仔细看一下这些新找到的序列，它们和输入序列之间的一致度其实并不比第一轮中找到的序列低。但是因为算法上的原因，它们被标准 BLAST 漏掉了。好在通过 PSI-BLAST，它们又被找回来了。 如果想要快速找到本轮搜索新找到的序列，也就是第一条标黄的序列，可以点击结果页面上方“skip to the first new sequence”链接（图2.46）。 图2.45 NCBI PSI-BLAST 第二轮搜索结果 图2.46 “skip to the first new sequence”链接

（三）NCBI PHI-BLAST

PHI-BLAST 和 PSI-BLAST 不同，PSI-BLAST 是撒大网搜索，而 PHI-BLAST 则是精准搜索。PHI 是 Pattern-Hit Initiated 首字母缩写，中文是模式识别。PHI-BLAST 能找到与输入序列相似的并符合某种特征模式的蛋白质序列。模式 Pattern 是对特征的描述。比如发生 N 糖基化位点的序列都符合这样一个特定的模式：发生糖基化的天冬酰胺，后面一定紧跟一个脯氨酸以外的氨基酸，再紧跟丝氨酸或者苏氨酸，再紧跟一个脯氨酸以外的氨基酸。 特定模式可通过正则表达式来表述。所谓正则表达式就是这句话的一个简约的规范性字符书写法。发生 N 糖基化位点的序列符合的特定模式翻译成正则表达式为 N{P}[ST]{P}。其中，N 是天冬酰胺，P 是脯氨酸，S 是丝氨酸，T 是苏氨酸。{}代表除大括号里的氨基酸以外的任意氨基酸，[ ]代表中括号中的任意一个氨基酸。PHI-BLAST 可以根据给入的正则表达式对搜索到的相似序列进行模式匹配，符合正则表达式的才会被作为结果输出。 正则表达式，{ }代表除什么以外，[ ]代表其中之一，x 代表任意字母，(3,7)代表 3 到 7 个某字符。那么正则表达式{L}GEx[GAS][LIVM]x(3,7)的意思是，除 L 以外的任意一个字母，紧接 G，再紧接 E，再接一个任意字符，之后是 GAS 中的任意一个，再接 LIVM 中的任意一个，最后再接 3 到 7 个任意字符。这种序列特征模式可能代表某个翻译后修饰的发生位点，也可以代表一个酶的活性位点，或者一个蛋白质家族的结构域、功能域等有重要意义的地方。 在 NCBI BLAST 工具的输入页面，当算法选择了 PHI-BLAST 之后，会自动出现模式输入框（图 1）。输入正则表达式 S[IVFL]TPS(2)（含义为：一个 S 后面紧接 IVFL 中的任意一个字母，再接 T，再接 P，再接两个 S）。这次搜索找到的相似序列中，只有符合该模式的才会被作为结果返回。 图2.46  NCBI PHI-BLAST 输入界面 结果显示只有几条严格符合模式的序列被找到了（图2.47）。此外，PHI-BLAST 可以和PSI-BLAST 联合使用，以找到更多符合模式的远房亲戚们。 图2.47  NCBI PHI-BLAST 结果界面

（四）其他 BLAST

除了前面三种BLAST，NCBI还开发了一个SMART-BLAST（聪明的 BLAST）。SMART-BLAST 可谓是标准 BLASTp 的简约强化版。操作极其简单，简单到只需要输入序列（图2.48）。 图2.48 NCBI SMART-BLAST 输入界面 SMART-BLAST 虽然操作简单，但返回的结果却并不简单（图2.49），它包括数据库中与输入序列最相似的三条序列，以及研究的最透彻的物种中可以展现一定进化关系的两条相似序列。图2.49中黄色的是你输入的序列，绿色的是研究的最好的模式物种中与你输入序列相似的序列。旁边还直接给出了这三条序列的系统发生树。总之，SMART-BLAST 可以帮你从一大堆结果中挑选出你最想要的。如果你很懒，或者你很茫然，可以试试这个聪明的 BLAST。 图2.49 NCBI SMART-BLAST 搜索结果 除了BLAST，还有其他的数据库相似性搜索工具，比如 WU-BLAST。WU 代表Washington University，它比 BLAST 更灵敏，在插入空位的算法上更灵活。再有 SSEARCH，它有点儿慢，但是比 BLAST 更准确。FASTA 也是一个搜索工具，它也是有点儿慢，但是对于 DNA 序列的比较比 BLAST 更准确，尤其适合短的序列。最早被 FASTA 使用的序列格式就叫 FASTA 格式，并沿用至今。

八、在线多序列比对工具

多序列比对的定义很简单，两条以上的生物序列进行的全局比对就是多序列比对。为了看清楚每一列的保守情况和理化性质，通常会给多序列比对根据不同的原则赋予丰富的色彩。目前世界上最流行的多序列比对工具是 CLUSTAL 系列、TCOFFEE和MUSCLE。其中CLUSTAL 系列使用率最高；TCOFFEE 最新，而且还有很多变形；MUSCLE 最快，而且胃口大，能接受的序列数量是其他工具比不了的。各大生物网站，比如 EBI、Expasy 等都可以在线使用这些工具。 多序列比对的主要用途： 1）我们可以通过多序列比对确定某一个未知序列是否属于某一个家族。 2）可以用多序列比对构建系统发生树，查看物种间或者序列间的进化关系。事实上，做多序列比对是构建系统发生树的必要步骤之一。 3）模式识别。一些特别保守的序列片段往往对应着重要的功能区。通过多序列比对，可以找到这些保守片段，并由此推测出潜在功能区。 4）可以把已知的有特殊功能的序列片段通过多序列比对做出匹配模型。然后根据这个模型推测未知的序列片段是否也具有这个功能。除此之外，多序列比对在生物信息学分析的很多方面都有应用，比如用来预测蛋白质的二级结构和三级结构，预测RNA 的二级结构等等。 在开始做多序列比对之前，需要注意以下几点： 1）做多序列比对的序列个数不能太多，一般10到15条序列刚好，最好不要超过 50 条。序列太多，任何软件都受不了。 2）关系太远的序列不适合做多序列比对。两两之间序列相似度低于30%的一组序列，做多序列比对要么做不出来，要么即使勉强做出来了，做得也是零七八碎，没有任何意义。 3）关系太近的序列不适合做多序列比对。两两之间序列相似度大于 90%的序列，有再多条都只等于一条。做出来的多序列比对无非就是把各条序列抄写了一遍，没有任何意义。 4）短序列受不了。多序列比对支持一组差不多长的序列，个别很短的序列纯属捣乱分子，受不了。第五条，有重复域的序列受不了。如果序列里包含重复片段，大多数多序列比对的程序都会出错，甚至崩溃。 此外，给序列起名字也要谨慎。 1）序列名字里不要有“空格”，用下划线代替“空格”是个好习惯。 2）不要使用特殊字符，比如中文，@，#之类的。 3）序列名字不要太长，最好在 15 个字符以内。 4）一组序列里，不要有重名的序列。

（一）EMBL——Clustal Omega

EMBL的多序列比对工具（http://www.ebi.ac.uk/Tools/msa）。EMBL 的多序列比对工具很多，包括CLUSTAL 系列、TCOFFEE、MUSCLE。EMBL比对工具中 CLUSTAL 系列的最新版本是Clustal Omega。 要比较的多条序列都用FASTA 格式书写。这样程序会根据“>”自动识别出每条序列以及它们的名字。这就是 FASTA格式的好处。进入 Clustal Omega 参数输入界面（图2.50），黏贴或者上传序列。 图2.50  EMBL-Clustal Omega 参数设置界面 在这个界面里可以点通过点击 More options，设置各种设置。在这个例子里，所有参数都使用默认值。参数里有输出格式（OUTPUT FORMAT）和输出顺序（ORDER）这两个参数。输出格式里可以选择常用的多序列比对格式。我们选标准的 Clustal 格式。这是最常见的多序列比对格式。输出顺序参数可以设定多序列比对中各个序列的排列顺序。“aligned”是按照比对过程中自动创建的计算顺序排列；“input”是按照输入序列的原始顺序排列。输入序列是按照 TLR10、9、8、7…这样的顺序排列的。输出顺序参数设定为 aligned，看看比对结果里序列的排列顺序是否发生了变化。 做多序列比对的时间要比双序列略长。序列越多，序列越长，则时间越长。Clustal 格式的输出结果如图2.51 所示。可以看到，比对中序列的排列顺序跟输入的时候不一样了，这是按照比对创建时的计算顺序排列的。请点击 Download Alignment File 保存将当前结果，以便后面章节进一步加工。保存的文件后缀名是“.clustal”。它是一个纯文本文件，用写字板或者记事本都可以打开。 如果想要添加色彩，点击“Show Colors”。之后，不同的氨基酸根据它们的理化性质 不同会显示出不同的颜色（图2.52） 除了颜色之外，多序列比对每一段的最后一行都有些星星点点的标记。这些标记和双序列比对中的竖线、双点、单点的意思类似，但并不完全相同。如果某一列是完全保守的一列，也就是说这一列里的字母完全相同，那么这一列的下面就打一个“*”。如果这一列的残基有大致相似的分子大小及相同的亲疏水性，也就是这一列的字母要么相同要么相似，没有不相似的，那么就打一个“：”。如果这一列残基的分子大小及亲疏水性被一定程度上保留了，但是有替换发生在不相似的残基间，也就是这一列的字母有相似的也有不相似的，那么就打一个“.”。 图2.51 EMBL-Clustal Omega 多序列比对结果 图2.52 EMBL-Clustal Omega 多序列比对结果 Result Summary 标签里，给出了全部结果信息的下载列表和一个Jalview 的按钮，Jalview 是多序列比对编辑软件（图2.53）。在下载列表里，如果打开“Percent Identity Matrix”链接，可以得到所有序列两两之间的一致度矩阵。一致度矩阵的第一行省略掉了。它和第一列完全相同，都是序列的名字并且按照相同的顺序排列。所以这个矩阵是以对角线为轴对称的，并且对角线上是某条序列自己和自己的一致度，都是 100%。这个矩阵可以帮助我们更好的了解这些序列之间的关系。比如我们可以从中发现，一致度最高的一对序列是 TLR1 和 TLR6。 图2.53  Result Summary：Percent Identity Marix 除了通过一致度矩阵了解序列间的关系，还可以通过 Phylogenetic Tree 标签下的 GuideTree 清楚的看出哪条序列和哪条序列更相似（图2.54）。Phylogenetic Tree 翻译成中文是系统发生树。但是这里要特别注意，这不是真正意义上的系统发生树！它只是在创建多序列比对的过程中用到的树（Guide Tree），没有经过距离校正，所以不能当作系统发生树来使用。如果想要根据多序列比对结果构建系统发生树，可以在 Alignments 标签下，点击“Send toClustalW2_Phylogeny”链接，把做好的多序列比对发送给专门做系统发生树的工具。 图2.54 Phylogenetic Tree 标签

（二）TCOFFEE——Expresso

TCOFFEE 是一个非常流行的多序列比对工具。TCOFFEE 与 CLUSTAL 系列在所使用的算法上类似，准确度上比 CLUSTAL 系列略高，但计算耗时也比 CLUSTAL 系列略高。最关键的是 TCOFFEE 有很多种变形，也就是说它有更多的功能。许多网站都提供 TCOFFEE的在线使用，比如 EMBL 的多序列比对工具里就有TCOFFEE。TCOFFEE也有自己的网站（http://tcoffee.crg.cat）。TCOFFEE 本身是一个标准的多序列比对工具，跟CLUSTAL没有什么区别。我们来看它的变形，也就是根据比对序列种类的不同，TCOFFEE网站下特有的比对工具（图2.55）。 针对蛋白质序列的比对工具，除了 TCOFFEE 以外，还有 Expresso，M-Coffee, TM-Coffee以及 PSI-Coffee。其中，Expresso 最有特色，它是为序列加入结构信息后再做多序列比对的工具。因为有结构信息的辅助，它可以大大提高比对的准确度。M-Coffee 可以把多个比对的结果整合成一个。TM-Coffee 专为穿膜蛋白打造，PSI-Coffee 专为远源序列打造。同样的还有针对RNA和DNA序列的Coffee。 做 Expresso 的序列我们选用网站提供的示例序列（图2.56）。Show more options 下，可以通过各种方式给入输入序列的结构信息。如果你有这些序列现成的结构文件，也就是 PDB文件，可以直接把它们上传上来。三条序列对应三个上传链接。可以上传的结构文件不只限于 PDB 数据库下载的，也包括还未正式发表的解析结构或者计算机预测的结构，只要是用PDB 文件格式保存的，都可以。 如果没有现成的结构文件，但是这些序列在 PDB 数据库里有对应结构的话，你可以从接下来的输入框里，按照规定的写法，指定哪条序列对应 PDB 数据库中的哪个结构（图2.57），如果这里输入了信息，Expresso 会自动从 PDB 网站下载指定的 PDB 文件。那些已经本地上传的结构，Expresso 也会根据序列信息自动匹配出它是哪条序列的结构，不需要再在这里列出了。如果对结构信息一无所知，只需要将“MODE_PDB”钩选。 图2.55 TCOFFEE 网站下特有的比对工具 图2.56  TCOFFEE Expresso 序列和结构信息输入界面 图2.57  Expresso 其他给入结构信息的方式 之后，Expresso 就会自己在网络上为所有没有指定结构信息的序列搜索相应的结构。你提供给 Expresso 的结构信息越多，计算时间就会越短；你提供结构信息越少，计算时间就会越长。如果只勾选了“MODE_PDB”，那么需要等待的时间会很长，因为 Expresso 首先要搜索，然后要下载，最后要计算。因此，留下你的Email信息是很有必要的。 比对结果的页面链接会发送至邮箱。打开链接后会得到Expresso做出的比对结果（图2.58）。TCOFFEE 系列各比对工具做出的多序列比对的颜色代表比对质量的好坏。越红质量越好，越蓝质量越坏。这次的比对结果非常令人满意。如果配合上这些序列二级结构信息看一下的话，你会发现，螺旋和螺旋很好的对在了一起，折叠和折叠很好的对在了一起。 图2.58 Expresso 多序列比对结果 同样的序列做普通的 TCOFFEE，质量远不如 Expresso（图2.59）。可以看到二级结构全部错位。所以，如果你有序列的结构信息的话，用 Expresso 相比用普通的比对工具会大大提高比对质量。 图2.59  TCOFFEE 多序列比对结果

（三）多序列比对的保存格式

无论用哪种多序列比对工具，CLUSTAL 系列也好，TCOFFEE 系列也罢，都提供多种格式的多序列比对结果。比如有漂亮的网页格式的，标准的 Clustal 格式的，还有写完一条序列再写下一条的 FASTA 格式的，以及方便下一步建树使用的 Phylip 格式。 要保存哪种格式主要看你下一步要干什么。在选择保存格式之前，需要问自己几个问题。 1）你选的这个格式大多数软件都支持吗？ 2）你的同事能用吗？ 3）你需要的信息这个格式都提供了吗？ 4）这个格式适合进一步加工吗？ 通过回答这几个问题来确定你最终需要的格式。如果比对工具输出的格式里没有你想要的哪种，可以通过第三方软件进行格式转换。比如 fmtseq 工具（http://www.bioinformatics.org/JaMBW/1/2），它可以实现 20 多种格式间的转换。

（四）多序列比对的编辑和发布：Jalview 的介绍和操作

在 EMBL Clustal Omega 比对结果的 Result Summary 标签下有Jalview按钮。这个按钮可以快速启动 Jalview，但这里启动的在线版本功能不完整。完全版的 jalview 可以从 Jalview 官网（http://www.jalview.org）在线启动，或者下载安装到本地。 通过 Jalview 除了可以加工多序列比对，还可以针对比对中的序列做各种各样的分析、比如构建系统发生树、预测蛋白质二级结构、查看结构域家族、从 PDB 数据库中查询三级结构等。 点击 File 菜单- Input Alignment- From File - 打开我们用Clustal Omega做出并保存的多序列比对结果（clustal格式文件）。因为“.clustal”不是 Jalview 熟悉的后缀名，所以需要把文件类型改成“所有文件”才能看到它（图2.60）。 图2.60  打开多序列比对 在打开的多序列比对窗口的下方有三行柱状图（图2.61）。它们体现了比对中每个位置的保守度高低（Conservation）、比对质量高低（Score）、以及共有序列（Consensus）。从保守度行，可以很清楚的找到保守区大致的位置。共有序列指的是某一列出现频率最高的那个字母，比如第 58 列中 W 出现的频率最高，是 100%。如果某一列拥有的最 图2.61 多序列比对窗口 高出现频率的字母是两个或两个以上的话，会以“+”显示。把鼠标放在“+”上就可以看到是哪些字母出现的频率一样高。共有序列可以一定程度上体现出某个保守区域所具有的序列特征。以后如果看到和这段序列长相极其相似的序列，它很可能能跟这个保守区的功能相似。 在Colour菜单下有很多种颜色方案。能够和保守度这一行柱状图配合的颜色方案是 Percentage Identity。选了这个颜色方案之后，每一列会根据这一列的保守度用深浅不同的蓝色表示。蓝色越深说明这一列越保守，反之越不保守。再配合 Colour 菜单下的“By Conservation”参数，可以从弹出的参数设定窗口中设定保守程度达到百分之多少以上的才给赋予不同的蓝色，阈值以下的都是白色。 另一个较常用的颜色方案是 Clustal 系列配色方案。这个配色方案和 EMBL 多序列比对工具做出的结果页面里“Show Colors”之后的颜色方案是基本相同的。具体哪个氨基酸选用哪个颜色可以参见图2.62。我们从文献里看到的彩色多序列比对，大多是用的这种颜色方案。 图2.62  Clustal 系列颜色方案 除了给多序列比对上彩妆，有时还需要给它修理一下局部瑕疵，也就是对局部位置进行手动调整。比如，从前期实验我们得知，图2.63中方框所示区域的 TLR2、10、6、1 这四条序列的第 53 列应该往右挪一列，跟 TLR9、8、7 这三条序列的第 54 列对在一起。TLR2、10、6、1 这四条序列的第 53 列补空位。其他位置不动。 多序列比对的外观也很重要。默认情况下，多序列比对是单行显示的。由于序列长，需要拖动窗口拉条才能浏览全部。这样不利于查看分析，也不利于将导出的比对图片插入文献。如果想要让多序列比对根据 Jalview 窗口的宽度自动换行，可以在 Format 菜单下勾选“Wrap”。此外，还可以通过“Font…”窗口对字体格式、大小等进行调整。如果你只需要多序列比对，而不需要有关保守度等的注释行。可以关闭 View 标签下的“Show annotations”选项，以达到去掉注释行的目的。 Jalview 除了有编辑多序列比对的功能还有很多分析功能。比如，可以按照序列的名字、两两一致度或其他规则给比对中的序列重新排序以及为选中的两条序列做双序列全局比对（图2.64）、为选中的一组序列计算各种系统发生树（图2.65）、或者用在线软件为某一条序列预测二级结构（图2.66）。Web service 菜单下的所有功能都需要网络支持才能运行。 图2.63 多序列比对局部位置调整 图2.64 序列排序和双序列全局比对 2.65  计算系统发生树 图2.66  预测蛋白质二级结构 除了 Jalview，还有很多比对美化工具（图2.67）。Boxshade 擅长黑白制图。因很多学术期刊只收取彩图的编辑费，所以黑白图可以节省科研经费。ESPript 的功能十分强大。MView擅长把彩色多序列比对转换成 HTML 源代码。这样就可以将它直接插入网页，并方便以文本形式选取。 图2.67 多序列比对编辑工具列表

九、寻找保守区域

（一）序列标识图 WebLogo

如果用一句话来描述你究竟想从多序列比对中得到什么，答案是你想要找到序列中重要的位置。说得更专业一点，就是要找到保守区域。可以借助软件来更好的寻来保守区域。 序列标识图（sequence logo）就是序列的 logo，它是以图形的方式依次绘出序列比对中各个位置上出现的残基，每个位置上残基的累积可以反应出该位置上残基的一致性。每个残基对应图形字符的大小与残基在该位置上出现的频率成正比。 但图形字符的大小并不等于频率百分比，而是经过简单统计计算后转化的结果。图2.68 是用一款流行的软件 WebLogo 创建的序列标识图。 图2.68 序列标识图 要创建序列标识图，首先需要一个多序列比对。多序列比对中的一列对应序列标识图中的一个位置。然后分别计算每一列中不同残基出现的频率，再根据以下公式（图2.69）把频率转换成高度值，最后根据高度值写出不同残基的彩色字母图形。 图2.69 频率转换成高度值 如果某一列非常保守，字母高度就高。反之，如果某一列没有什么特征，各种残基都有出现，杂乱无章，那么就会看到一堆比较矮的字母摞在一起。这里再次强调，字母的高度和它在某一列中出现的频率成正比，但是并不等于频率。试想一下，如果字母高度就是频率的话，那么序列标识图中每个位置上字母摞起来的总高度应该是一样的，都是 100%。但是从图2.69 中可以看到，序列标识图上每个位置字母摞起来的总高度是不一样的，这是因为在字母高度的计算过程中涉及了熵值。某一列中字母出现的情况越混乱，熵值越大，字母越矮。字母出现的情况越有规律，熵值越小，字母越高。所以序列标识图可以很好的展现多序列比对中每一列的保守程度，即，它们是杂乱无章的，还有有规律可循的。并且把可循的规律图形化的展现出来。这就是我们为什么要给序列打上 logo 的原因。 WebLogo 是一款在线创建序列标识图的软件（http://weblogo.threeplusone.com/） 主页面上点“Create your own logos”，然后输入多序列比对（图2.70）。WebLogo 可以接受大多数常见的多序列比对格式。示例文件 promoter.fasta 是一组启动子序列的多序列比对，以 FASTA 格式存储。FASTA 格式的多序列比对要求把多序列比对中的每一条序列连同插入的空位一起按 FASTA 格式书写，写完一条序列再写下一条。这和之前讲过的 Clustal 格式不太一样。在序列输入框的下方可以设置不同参数，以定义序列标识图的样式，比如设置序列标识图的创建范围、定义字母的颜色方案等。保持所有参数默认，点“Create WebLogo”。 图2.70  WebLogo 输入页面 图2.71为创建出的序列标识图。从图中可以清晰的看到：输入的这些启动子序列上TATA-Box 的共有特征序列，以及它们出现的位置。 图2.71  WebLogo 结果页面    

（二）序列基序 MEME

MEME 是一款寻找序列基序（motif）的软件。在核酸或蛋白质序列中存在一些有特定模式的序列片段，这些片段称为序列的基序（motif）。序列的基序与生物功能密切相关。比如，发生 N 糖基化位点的基序：发生糖基化的天冬酰胺后面一定紧跟一个脯氨酸以外的氨基酸，再紧跟丝氨酸或者苏氨酸，再紧跟一个脯氨酸以外的氨基酸。这个特定模式可通过正则表达式来规范描述，也可以通过序列标识图来直观描述。基序的发现要通过大量相关序列的分析。MEME 就是一款可以自动从一组相关的核酸或蛋白质序列中发现序列基序的软件。 MEME 是 The MEME Suite 在线软件套装中的一员（http://meme-suite.org/）。MEME 的使用非常简单，只需要将待分析的序列上传即可（图2.72）。而且，上传的序列为原始序列，不需要提前为它们做多序列比对。你也可以指定返回排名前几的基序。MEME 的等待时间稍长，大约 10 分钟以上，所以最好留下邮箱。 图2.72  MEME 输入页面 Meme 的返回结果被保存成各种格式：HTML、XML、test 等。便于在线查看的是“MEMEHTML output”，即网页格式。 网页格式的 MEME 结果页面中，给出了找到的排名前三的基序（图2.72）。它们以序列标识图的形式展现出来。同时还提供这三个基序在每条序列中的大体位置。如果要进一步了解某个基序，可以点击序列标识图右侧的“More”下面的“”箭头，以查看详细（图2.73）。点击后，会得到大比例序列标识图，以及该基序在每条序列中对应的序列片段和它们出现的具体位置。此外，还可以点击序列标识图右侧的“Submit/Download”下面的“”箭头（图2.74），将某个基序提交至各种数据库，并进行针对该基序的序列相似性搜索，已找到数据库中 含有该基序的序列，进而推测该基序的功能。这步操作是通过 The MEME Suite 软件套装下的另一个软件 FIMO来实现的。 图2.72  MEME HTML 结果页面 图2.73  More 链接查看基序详情 图2.74 提交基序给FIMO进行数据库相似性搜索

（三）PRINTS 指纹图谱数据库

目前，科学家已经对现有的蛋白质序列进行了充分的研究，而且早已发现并总结了这些序列上的重要基序。相关研究成果汇入了 PRINTS 蛋白质序列指纹图谱数据库（http://www.bioinf.manchester.ac.uk/dbbrowser/PRINTS/）。所谓蛋白质的指纹是指一组保守的序列基序，用于刻画蛋白质家族的特征。这些基序由多序列比对结果获得，且它们在氨基酸序列水平上是不相邻的，但是在三维结构中可能紧密地结合在一起。PRINTS 数据库存储了目前已发现的绝大多数蛋白质家族的指纹图谱。对于一个陌生的蛋白质，只要看看它的序列是否符合某个蛋白质家族的图谱就可以对它进行分类并预测它的功能。 要浏览 PRINTS 数据库，可以输入数据库编号、关键词、或标题等以查找某一个指纹图谱。比如点击“By text”通过关键词搜索（图2.75）。输入条中输入“TRANSFERRIN”，也就是搜索转铁蛋白家族的图谱。搜索返回转铁蛋白家族的指纹图谱链接。 点击结果页面中的“TRANSFERRIN”链接后，会显示包括指纹图谱的基本信息、与其他数据库之间的交叉链接、构建指纹图谱所使用的蛋白质序列、以及指纹图谱中每个基序等具体信息（图2.76）。 点击“View alignment”链接后，可以看到创建指纹图谱所使用的多序列比对（图2.77）。 点击“View structure”链接后，网页会打开一个三维视图插件，并以该家族中某一特征蛋白质具有的三维结构为例，在线显示指纹图谱中各个基序在三维结构中的位置（图2.78）。从该三维结构图中可以看出，紫色的基序在氨基酸序列水平上并不相邻，但是在三维空间结构中是紧密联系在一起的，并形成蛋白质的重要功能区。 图2.75  关键词搜索转铁蛋白家族图谱 图2.76 TRANSFERRIN 结果页面 图2.77  View alignment 结果页面 图2.78  View structure 结果页面 除了浏览某一指纹图谱，PRINTS 还提供指纹匹配服务。也就是搜索某一序列所匹配的指纹图谱。此功能通过 PRINTS 主页也上的“FPScan”链接实现（图2.79）。注意输入的待搜索序列只能是“a raw sequence”，也就是纯序列。换言之，FASTA 格式中带大于号的第一行不能拷贝进输入框。 图2.79 FPScan 输入页面 提交后返回的结果页面中，跟输入序列匹配的指纹图谱，根据匹配得分的高低被排列出来（只列出前十名）（图2.80）。此外，还单独列出了排名前三的指纹图谱。由此可知，得分最高的是视紫红质家族的指纹图谱。 图2.80 FPScan结果页面 点击排名第一的视紫红质家族的“Graphic”链接，可以得到该家族指纹图谱中各个基序在输入序列中所匹配的位置（图2.81）。结果页面的下部还提供了视紫红质家族的 6 个基序在输入序列中所对应的具体序列片段。由此，可以推测，输入序列属于视紫红质家族，并具有该家族蛋白质的主要功能。 图2.81 视紫红质家族指纹图谱中各个基序在输入序列中所匹配的位置

第三章 分子进化与系统发生

拉马提出进废退理论。他说生物经常使用的器官会逐渐发达，不经常使用的器官会逐渐退化。而且这种后天获得的性状是可以遗传的，因此生物可以把后天锻炼的成果遗传给下一代。达尔文认为，所有的生物物种都存在趋利的适应性变化。并且这些适应性变化会通过一种他称之为“自然选择”的过程遗传下去。而大自然是数百万年间推动演变进化的唯一力量。这个理念让达尔文成为了现代生物学之父。达尔文认为“物竞天择”意味着，存在一种最初的生物，之后通过某种方式得到了改良。 拉马克和达尔文告诉我们，进化很重要！ 传统上有两种研究进化的方法。一种是看死的，一种是看活的。看死的，也就是看化石，这是进化最直接，最确凿的证据。第二种研究方法是利用比较形态学、比较解剖学和生理学等手段，确定大致的进化框架。这种方法比第一种方法更容易实现。但是这种方法仅局限于大致的框架，很多细节是存在争议的。这两种传统方法都有局限性，如今研究的是分子进化。也就是利用软件，从分子水平上构建物种的进化树。这里说的分子水平是指 DNA、RNA、以及蛋白质序列。 分子进化理论是 1962 年美国科学家 Linus Pauling提出的。这种理论与传统研究方法的最大区别是，它研究的是 DNA、RNA 以及蛋白质序列这些分子水平上的信息，而不是物种的外在特征。并且基于某一个特定的分子在不同物种中的序列差异来构建系统发生树。此外，分子进化有两个基本的假设条件，只有接受这两个假设，分子进化理论才能得以实施。第一、DNA、RNA 或蛋白质序列包含了物种的所有进化史信息。第二是分子钟理论。这个理论说的是，一个特定基因或蛋白质的进化变异速度在不同物种中是基本恒定的。所谓变异速度是指一定时间内不同碱基或氨基酸突变的个数。这个进化变异速度被认为是恒定的，跟物种没有关系。所以，拿蛋白质来说，两个蛋白质在序列上越相似，他们距离共同祖先就越近。分子钟理论是进化研究领域被普遍认可的理论，但是至今也没有直接的证据证实。 同源（Homologs）：来源于共同祖先的相似序列为同源序列。也就是说，相似序列有两种，一种是来源于共同祖先的，那么他们可以叫同源，另一种不是来源于共同祖先的，那么他们尽管相似也不能叫同源。第二种情况出现的概率虽然低，但还是存在的，所以相似序列并不一定是同源序列。同源又分为三种，直系同源，旁系同源和异同源。同源只是对性质的一种判定，只能定性描述，不能定量描述。 直系同源（Orthologs）是指来自于不同物种的由垂直家系，也就是物种形成，进化而来的基因，并且典型的保留与原始基因相同的功能。也就是说，随着进化分支，一个基因进入了不同的物种，并保留了原有功能。这时，不同物种中的这个基因就属于直系同源。 旁系同源（Paralogs）是指在同一物种中的来源于基因复制的基因，可能会进化出新的但与原功能相关的功能来。 异同源（Xenologs）是指通过水平基因转移，来源于共生或病毒侵染所产生的相似基因。异同源的产生不是垂直进化而来的，也不是平行复制产生的，而是由于原核生物与真核生物的接触，比如病毒感染，在跨度巨大的物种间跳跃转移产生的。

一、系统进化树

研究分子进化所要构建的系统发生树（Phylogenetic tree），也叫分子树。 对于一个未知的基因或蛋白质序列，可以利用系统发生树确定与其亲缘关系最近的物种。比如你得到了一个新发现的细菌的核糖体 RNA，你可以把它跟所有已知的核糖体 RNA 放在一起，然后用他们构建一棵系统发生树。这样就可以从树上推测出谁和这个新细菌的关系最近。系统发生树还可以预测一个新发现的基因或蛋白质的功能。 系统发生树还分为有根树和无根树（图3.1），顾名思义，有根树就是有根，无根树就是无根。其实两者是可以互换的。如果我们按住无根树上某一个点，然后用把梳子将树上所有的枝条都以这个点为中心向右梳理，就能把它梳成有根树的样子。按住的这个点就是根。所以对于一棵树来说，根的位置是主观的，你想让他在哪它就在哪里。但是你不能随意指定哪个内节点当根，毕竟根有其自身的生物学意义，它应该是所有叶子的共同祖先。那么我们如何确定根的位置呢？可以通过外类群（outgroup）来确定，从而把无根树变成有根树。 图3.1  无根树和有根树 有根树反映了树上基因或蛋白质进化的时间顺序,通过分析有根树的树枝的长度，可以了解不同的基因或蛋白质以什么方式和速率进化。而无根树只反映分类单元之间的距离,而不涉及谁是谁的祖先问题。做有根树需要指定外类群。所谓外类群，就是你所研究的内容之外的一个群。比如你要分析某一个基因在不同人种间的进化关系，那就可以额外选择黑猩猩加入进来，作为外类群一同参与建树。 构建系统发生树的方法很多。最常用的有基于距离的构建方法，包括非加权分组平均法（Unweighted Pair Group Method with Arithmetic mean，UPGMA），最近邻居法（NeighborJoining method，NJ），最小二乘法（Generalized Least Squares，GLS）等。还有最大简约法（Maximum Parsimony，MP），最大似然法（Maximum Likelihood，ML），贝叶斯推断法（Bayesian Inference，BI）等。 从计算速度来看，最快的是基于距离的方法，几十条序列几秒钟即可完成。其次是最大简约法。最大似然法就要慢得多。最慢的是贝叶斯法。但是从计算准确度来看，算得最慢的贝叶斯法确是最准确，而算得最快的基于距离法结果确是最粗糙。从实用的角度，建议使用最大似然法。因为这种方法无论从速度还是准确度都比较适中。 目前流行的建树软件（图3.2），比如 PHILIP 和 MEGA，基本能够包括上述所有算法。如果想要构建 ML 树，也可以尝试专门构建 ML 树的 PHYML。贝叶斯的算法以 MrBayes 为代表，只是计算速度比较慢。如果构建的系统发生树要用于发表生物信息学领域的文章，需要两种以上的构建方法锁定同一个结果才能审稿通过。如果是用于发表以生物实验为主的文章用一种构建方法就可以了。 虽然软件可以快速自动地完成系统发生树的构建，但是对于基本算法的了解还是必不可少的。以非加权分组平均法（UPGMA法）为例，介绍如何通过计算所有序列两两间的距离，再根据距离远近构建系统发生树。序列两两间的距离可以用双序列比对得出的一致度/相似度代表，或用其他简化值代替。 图3.2 构建系统发生树的软件 比如，有如下 A、B、C、D 四条序列： A：TAGG B：TACG C：AAGC D：AGCC 在接下来的例子里，我们用序列间不同的碱基数目作为序列间遗传距离的度量。首先，计算出每两条序列间有几个碱基不同，并以用矩阵的形式记录下这些距离，找出距离最小的一对序列。A和B之间的距离最小，d[AB]=1。然后将 A 与 B 合并聚集，其分支点为 d[AB]/2=1/2=0.5。即，A、B之间的距离等于1，从中间折叠后每边各 0.5。现在，把(AB)看成一个整体，分别计算它们与C和D的距离。(AB)和C的距离等于A和 C 的距离加上B和C的距离除以2，即，d[(AB)C]=(d[AC]+d[BC])/2=(2+3)/2=2.5。同样，(AB)和D的距离等于A和D的距离加上B和D的距离除以2，即，d[(AB)D]=(d[AD]+d[BD])/2=(4+3)/2=3.5。据此，计算出新的距离矩阵，并找出新矩阵中最小的距离。C 和 D 之间的距离最小，d[CD]=2。将C和D进行合并聚集，其分支点为 d[CD]/2=2/2=1（图3.3）。 图3.3  非加权分组平均法 接下来，把(CD)看成一个整体，计算它们与(AB)之间的距离。(CD)与(AB)之间的距离等于C和(AB)的距离加上D和(AB)的距离除以2，即，d[(CD)(AB)]= (d[C(AB)]+d[D(AB)])/2=(2.5+3.5)/2=3。最后，将(AB)与(CD)进行合并聚集，归为一类，分支点为 d[(CD)(AB)]/2=3/2=1.5。这样，A、B、C、D 四条序列的系统发生树就构建好了。树上，枝的长短直接反应了它们与共同祖先的距离。 分子进化的研究对象是核酸和蛋白质序列。如果要研究某个基因的进化，是应该选用它的 DNA 序列，还是翻译后的蛋白质序列呢？序列的选取要遵循以下原则：1）如果 DNA 序列两两间的一致度≥70%，选用 DNA 序列。因为，如果 DNA 序列都如此相似，它们对应的蛋白质序列会相似到几乎看不出区别。这对于构建系统发生树是不利的。所以这种情况选用 DNA 序列，而不选蛋白质序列。2）如果 DNA 序列两两间的一致度<70%，DNA 序列和蛋白质序列都可以选用。

二、利用MEGA软件构建系统发育树

MEGA7 是完全的图形化界面操作（http://www.megasoftware.net/），以10 条人的不同 Toll 样受体胞内域的氨基酸序列为例。 点击 MEGA 主窗口上的 File —Open A File，找到自己保存好的序列文件（fasta格式）。这时，MEGA 会询问以何种方式打开（想要做系统发生树先要做多序列比对，然后把多序列比对结果提交给建树软件去建树。所以，此时可以输入一个现有的多序列比对，也可以输入原始序列，让MEGA先来做多序列比对，再建树。）在这里选择“Align”，之后，在弹出的 Alignment Explorer窗口上点击 Alignment —Align by ClustalW。MEGA 提供 ClustalW 和 Muscle 两种多序列比对方法。这里选择熟悉的 ClustalW 方法。弹出窗口询问“Nothing selected for alignment. Select all? （是否要选择所有序列来做多序 列比对) ”，选择OK。之后，弹出多序列比对参数设置窗口（图3.4）。这个窗口和EMBL 的多序列比对在线工具一样，可以设置替换记分矩阵、不同的空位罚分等参数。MEGA 的所有默认参数都不是随便设置的，这些经过反复考量默认设置好的参数保证了 MEGA 傻瓜机全自动档的品质。所以，当你无从下手的时候，直接点 OK，接受这些默认参数，开始计算多序列比对。 图3.4  多序列比对参数设置窗口 将这个多序列比对作为中间结果保存下来。在Alignment Explorer 窗口上点 Data  Export Alignment MEGA Format。注意这里一定选 MEGA format 以方便MEGA 继续加工。其他格式适用于其他软件。选择 MEGA 格式文件，文件自动赋予后缀名“.meg”。保存后，弹出窗口，要求给这组数据命名。生成的.meg格式文件属纯文本文件，可以用记事本打开查看，里面是用 MEGA 固有格式书写的多序列比对结果。可以双击“TIR.meg”文件将其导入 MEGA。也可直接将其拖入 MEGA 主窗口中。拖入后，主窗口里增加了一个TA图标。点击该图标，弹出新窗口Sequence Data Explorer，其中显示出 MEGA 格式的多序列比对。再点击 Sequence Data Explorer 窗口上的TA按钮。点击后，多序列比对中的 10 条序列之上增加了一行。这一行是根据多序列比对结果分析得出共有序列（consensus sequence），也就是每一列里出现次数最多的那个字母。多序列比对中，每一列里的字母如果和共有序列的字母相同则打点，不同则标出不同的字母，空位还是空位（图3.5）。 图3.5 如果你想进一步了解这些序列的保守程度，可以点击C按钮，以黄色标记保守列，或者点击V按钮，以标黄标记不保守列。通过对比对结果的进一步分析，可以淘汰掉一些序列。比如海选的序列里有有一些明显不合群的序列，就可以把它们前面的钩去掉，不让它们参与建树，以免影响建树质量。此外，还可以对这些序列进行分组标记。点击分组按钮（图3.5，TA左边第二个）会弹出一个对话框（图3.6）。 图3.6 分组对话框 我们让TLR3、7、8、9这四个内质网上的蛋白质序列在一组，点击＋号图标，更改组名为“ER”。然后按住 Ctrl 键同时选中 Ungrouped Taxa 列表里的 TLR3、7、8、9 四条序列的名字。选中后点击＋号图标，将选中的序列归入ER组。再次点击＋号图标，创建第二组，更名为CM，即细胞膜上的蛋白质。按住 Ctrl 选中 Ungrouped Taxa 列表里剩下的所有序列，点击＋号图标，将它们归入CM组。当序列数量较多的时候，人为分组，可以从树上更加清楚的看出组内哪些成员叛逃了去了别组。此外，输入序列的名字较长。这样的名字如果作为构建的系统发生树上叶子的名字，会破坏树的外观也不利于信息的解读。因为，需要人为修改一下序列的名字，可点击－号下的图标进行更改，把名字改为能区分彼此的关键词，比如只保留 TLR*。全部更改好之后，点Save，准备工作全部完成。 点击 MEGA 主窗口上的 Phylogeny 下拉菜单，之前提到的各种建树方法这里基本都有。选 Neighbor Joining（最近邻居法）。弹出窗口询问是否使用当前 TIR.meg 里面的数据，选Yes。接下来，弹出参数设置窗口（Analysis Preferences）（图3.7）。参数设置对构建的系统发生树的准确程度来说非常重要。在树构建好之后，还经常需要根据树的具体情况，重新设置参数，并重新建树，如此反复，直至结果令人满意为止。同样的，如果对于参数依然摸不着头脑，不知道应该怎样设置，那么继续使用 MEGA 的自动挡，接受默认参数，也能做出基本满意的系统发生树。如果想要做出更好的树，至少应该掌握其中三个参数的设置。第一个参数是Test of phylogeny（建树的检验方法）。默认值为不进行检验。这显然比较偷工减料，不检验怎么知道建出来的树质量如何，是否可信！因此，检验方法选常用的Bootstrp method（步长检验）。步长检验需要设定检验次数，通常为100 的倍数，比如设置为500。 图3.7  参数设置窗口 步长检验是根据所选的建树方法，计算并绘制指定次数株系统发生树。因为大多数建树方法的核心算法都是统计概率模型，所以每次计算出的树都会有所差别。而建好的系统发生树上每个节点上都会标记一个数字，它代表了指定次数次计算所得出的系统发生树中有百分之多少棵树都含有这一节点。一般来说，绝大多数节点上的数值都大于70%的树才可信。个别低于70%的节点可以暂且容忍，或通过添加，删减序列来改善质量。 第二个参数是，Substitution Model。它是选择计算遗传距离时使用的计算模型。理论上应该尝试各种模型，根据检验结果选择最合适的模型进行计算。但在实际操作中，可先尝试选用较简单的距离模型，比如 p-distance。第三个参数是 Gap/Missing Data Treatment。大多数建树方法会要求删除多序列比对中含有空位的列。但是根据遗传距离度量方法的不同，删除原则也不同。如果是以序列间不同残基的个数来度量遗传距离的话，这里需要选择 Complete deletion（全部删除）。如果是其他方法，比如这里选用的 NJ 方法，可以选择 Partial deletion（部分删除）。删除程度定在 50%，即，保留一半含有空位的列。 按照以上方案参数设置后，点击Compute开始构建系统发生树。经过一番计算之后，新窗口Tree Explorer 里展示的就是创建好的系统发生树（图3.8）。 图3.8  系统发生树 Original Tree 是步长检验构建的 500 株树中的一株，未经过多棵树合并，所以树枝的长短可以精确代表遗传距离。此外，从这株树也可以看出之前的人为分组情况是不是发生了意想不到的变化。比如，TLR5 似乎脱离了CM 组，成为了外类群，从而确定了树根。树构建好之后，外形也许还不太令人满意。比如也许你想要将树的外形改成圆形或三角型。 可以通过树状按钮选择。或者你想要调整树枝的粗细或字体的大小，可以从 View 下拉菜单下的 Option 选项卡中调整。调整好之后，就可以把这棵树保存成图片了。保存图片可以点Image 下拉菜单，选择保存格式。或者将窗口放大截图。

第四章 蛋白质结构分析和预测

蛋白质的结构可分为四级。一级结构也就是氨基酸序列，二级结构是周期性的结构构象，比如α螺旋β折叠等。三级结构是整条多肽链的三维空间结构。四级结构是几个蛋白质分子形成的复合体结构，比如三聚体，四聚体等。 蛋白质经过折叠后会形成规则的片段，这些规则的片段构成了蛋白质的二级结构单元。三种常见的二级结构单元包括螺旋、β折叠、和转角。螺旋中最常见的就是α螺旋，但不只有α螺旋，还有其他的螺旋，比如 3 转角螺旋，5 转角螺旋等。β折叠由平行排列的β折片组成。这些折片在序列上可能相隔很远，但是在空间结构上并排在一起，彼此间形成氢键。除了螺旋和折叠外，蛋白质结构中还存在大量的无规律松散结构 coil。如果这些无规律的肽链突然发生了急转弯，这个转弯结构就叫做β转角。 蛋白质的二级结构经常用图形来形象的描述。比如图4.1中黄色的箭头代表对应的氨基酸具有β折片结构。波浪线代表螺旋结构，小鼓包是转角。此外，以字母形式书写的二级结构序列能够更加精准的描述。其中，E 代表β折叠，H 代表α螺旋，T 代表转角。没有写任何字母的地方是松散的 coil 结构。   图4.1  二级结构的图形及字母形式表述

一、 蛋白质的二级结构

（一）DSSP指认

截至目前，并非所有已发现的蛋白质都有明确的二级结构信息。只有通过实验方法已经解析出三维空间结构的蛋白质才有二级结构信息。这些二级结构信息是研究人员根据DSSP（蛋白质二级结构定义词典）将三级结构里的二级结构单元指认出来的。然后再按照规定的格式，记录下蛋白质中每个氨基酸处于哪种二级结构单元。这样一个记录蛋白质二级结构信息的文件叫做DSSP文件。蛋白质结构数据库PDB中的每一个蛋白质三级结构都有自己对应的 DSSP 文件。DSSP 文件里不同字母所代表的不同二级结构单元和PDB里面的记录方式是统一的。 DSSP 的主页上，Introduction 部分有一个 Web server 链接。这个链接很容易让人误以为可以通过它预测某条氨基酸序列的二级结构。这是不对的。DSSP 网站的 Web Server 可以指认蛋白质结构文件，也就是 PDB 文件中的二级结构，并创建出相应的 DSSP 文件。提交的PDB 文件可以是用实验方法刚刚解析出来，还没有提交 PDB 数据库的蛋白质三级结构，也可以是用计算方法预测出来的蛋白质三级结构模型。总之，输入值必须是三级结构，而不是一级的氨基酸序列。至于那些已经提交到 PDB 数据库中的蛋白质结构对应的 DSSP 文件可以从DSSP网站（http://swift.cmbi.ru.nl/gv/dssp/）提供的fpt网址（ftp://ftp.cmbi.ru.nl/pub/software/dssp/）直接下载。 此外，PDB数据库也提供直接下载DSSP文件的网址（图4.2）。DSSP 文件虽然记录的信息全，但查看起来并不直观。如果想要更加方便直观的查看蛋白质的二级结构信息，可以去 PDB 数据库网站。 图4.2  PDB 数据库图形化二级结构和DSSP文件下载链接

（二）PDB获取

PDB 数据库中，一个蛋白质结构记录中的二级结构信息在 Sequence 标签下（图4.2）。从序列图形化部分可以看到二级结构对应在一级结构上的图形化表示。点击左侧的“ViewSequence & DSSP Image”可以获得直观的一级结构对二级结构的序列表示（图4.3）。图4.3中的序列有两行，上面的一行是一级结构，下面的是二级结构。这个页面看上去很不错，序列10个字母一间隔，50个字母一行，而且不同的二级结构还对应不同的字母颜色。但是在接下来的分析研究工作中，我们往往需要的是像氨基酸序列那样的 FASTA 格式的二级结构序列。想要从这个网页上单独保存下二级结构序列是很麻烦的事儿。需要一行一行的拷贝黏贴，还需要删除行号。有位困难的是去除其中的空格，因为很难区分是格式里的空格还是代表松散结构的空格。所以，这种形式的二级结构信息便于浏览，但是不便于保存。非常遗憾的是，PDB里没有现成的针对某一个蛋白质的 FASTA 格式二级结构序列下载链接。“Download FASTA File”链接只能下载 FASTA 格式的一级结构序列，也就是氨基酸序列。 图4.3  PDB 中的一级结构序列和二级结构序列对应图 此外，PDB数据库中有一个叫做“ss.txt”的文件：http://www.rcsb.org/pdb/files/ss.txt.gz(。这个文件里面有PDB所有蛋白质结构的一级和二级结构的 FASTA 格式序列。但是这个文件非常大！仅仅打开文件就要耗费许久时间，使用起来相当的不方便。可以通过一个小程序输出fasta格式文件（http://1.51.215.28/~gongj/biotools/）（图4.4）。只需要输入PDB ID，程序就会自动下载相应的DSSP文件，并从中抽取出一级和二级结构的序列信息，最后以FASTA格式输出。 图4.4  BioTools二级结构自动获取工具

（三）蛋白质的二级结构软件预测

用计算机软件来预测蛋白质的二级结构。预测的结果和真实情况会有一定出入，究竟差多少，取决于预测软件的准确度。可以预测蛋白质二级结构的软件很多，而且都可以在线使用。图4.5列出了一些常用的在线预测软件。 图4.5  常用蛋白质二级结构预测软件 以PSIPRED为例，看看如何从氨基酸序列预测二级结构（图4.6）。PSIPRED（http://bioinf.cs.ucl.ac.uk/psipred）是一个蛋白质序列分析平台，它不仅可以预测二级结构，还有很多其他分析功能，比如预测三级结构。选择第一个“PSIPRED”工具，即，预测蛋白质二级结构的工具。输入的氨基酸序列（见附件 psipred.fasta）在 PDB 数据库中已有对应的空间结构（PDB ID：3CIG），因此二级结构也是已知的。之所以输入这样一条序列，是为了将预测结构和真实结构进行比较，从而评估一下软件的预测准确度如何。预测一般需要30 分钟。可在线等待结果，也可以查收结果邮件。需要注意的是，像大多数在线软件一样，PSIPRED 不支持免费的商业邮箱，比如hotmail或者 QQ 邮箱等。此外，最好给预测任务起个名字。最后，点击“Predict”。 图4.6  PSIPRED预测工具输入界面 预测结果页面 Summary 标签下的内容如图4.7所示，，粉红色的位置是α螺旋出现的位置，黄色的是β折片，没有底色的是松散的 coil 结构。如果预测出有错乱的结构也会被标出。目前的二级结构预测软件，都只能预测出α螺旋和β折片，对其他不常见的二级结构单元并不进行预测 图4.7  PSIPRED预测工具结果界面Summary标签 结果页面中，PSIPRED 标签下有全面的结果描述（图4.8）。描述一共有四行，最下面是一级结构序列，往上是二级结构序列，再往上是二级结构图形化的描述，最上面的柱状图反应的是每个位置预测的可信度。 图4.8  PSIPRED 预测工具结果界面 PSIPRED 标签 从结果页面的Download标签下可以下载纯文本格式的结果文件（图4.9）。纯文本格式的结果更利于下一步分析加工。 现在，我们把 PSIPRED 预测结果和 DSSP 里的真实二级结构放在一起（图4.10）。绿色的是预测结果，粉色的是真实结构。可以看到图 6 中显示出的这部分结果里，绝大多数α螺旋和β折片，也就是H和E，都被正确的预测出来了，只有少数几个短的β折片没有被预测出来，准确度超过90%。如果某一个软件的预测结果心存疑虑，可以尝试用多个软件进行预测，最后把各个软件的结果综合在一起。如果大多数软件都预测出来某个位置有螺旋，那么这个位置就应该有螺旋。如果只有一个软件预测出某个位置有折片，而其他软件都没有预测出来，那么这个位置可能就是没有。至此，无论是已有空间结构的蛋白质，还是未知空间结构的蛋白质，我们都已经掌握了获得二级结构的方法了。 图4.9 PSIPRED 预测工具结果界面 Downloads 标签 4.10 PSIPRED预测结果和DSSP真实二级结构的比较

二、蛋白质的三级结构

（一）三级结构搜索

蛋白质的三级结构是指整条多肽链的三维空间结构，也就是包括碳骨架和侧链在内的所有原子的空间排列。第一个蛋白质的三维空间结构于 1958 年用 X-射线衍射法（X-ray Crystallography）测定。这种方法目前仍然是获取蛋白质三级结构的主要方法。PDB 数据库中绝大多数蛋白质结构都是用这种方法测定的。另一个测定蛋白质三维空间结构的方法是核磁共振法（Nuclear Magnetic Resonance, NMR）。无法结晶的蛋白质，可以利用核磁共振法在液体环境中进行结构测定。但是核磁共振法只能用于质量小于 70 千道尔顿的分子，大约对应 200 个氨基酸的长度。除此之外，还有一些不太常用的方法也可以测定分子的三维空间结构，比如冷冻电子显微镜技术（Cyro-Electron Microscopy）。无论用什么方法测定的空间结构，都要提交到 PDB 数据库。所以我们获取蛋白质三级结构最直接的办法就是去PDB 搜索（图4.11，http://www.rcsb.org/）。 图4.11  PDB 首页搜索结构 如图4.11所示，从PDB首页的搜索条里，可以通过搜索PDB ID、分子名称、作者姓名等关键词来查找蛋白质三级结构。此外，利用高级搜索工具，可以通过序列相似性搜索获得与输入序列在序列水平上相似的蛋白质的三级结构（图4.12）。搜索方法选 BLAST，输入序列，点击“Result Count”。此次搜索一共找到108个在序列水平上和输入序列相似的蛋白质。点击链接“108 PDB Entities” 图4.12 高级搜索工具BLAST搜索 搜索结果中，排在第一位结构是人的 dUTPase 蛋白的三维结构，PDB ID 为 2HQU（图4.13）。这个结构所对应的序列与输入序列中黄色片段之间的一致度是 100%。输入的序列中蓝色区域是信号肽。信号肽在蛋白质到达亚细胞定位之后就被切掉了，所以解析的成熟蛋白质结构里不会有这一段。此外，成熟肽段 N 端的一小部分，由于实验技术等原因，也没有被解析出来，这在 PDB 结构中是很常见的。有时，在序列中间也会有未解析出的断口。甚至有时，为了得到稳定的晶体状态，需要突变个别的氨基酸或者删除一截肽段。这些技术手段都会使得结构中的序列和蛋白质本身的序列有所差别。 图4.13  BLAST搜索结果

（二）三级结构可视化软件VMD

能够实现蛋白质三维结构可视化的软件非常多。比专业级的PyMOL（https://pymol.org/2/）。这个软件已经被世界上著名的生物医药软件公司“薛定谔公司（Schrödinger）”收购。这种专业级的可视化软件不仅能够做出非常漂亮的图片，它还有强大的插件支持各种各样的蛋白质结构分析，这款软件需要购买，如果你发表的文章里提到某些内容是使用PyMOL制作的，而文章中所有作者和作者单位都没有PyMOL的购买记录的话，你可能会面临薛定谔公司的追责。 下面给大家介绍一个功能同样强大的免费蛋白质三维结构可视化软件，VMD（http://www.ks.uiuc.edu/Research/vmd）。VMD由伊利诺伊大学研发。下载 VMD 需要先注册获得一个账户，之后就可以根据你的操作系统和机器配置选择合适的版本下载了。当然，如前所述，注册和下载对于非商业用途的用户都是免费的。VMD 的安装也极其简单。不需要预装任何语言环境，完全图形化安装过程，绝对可以轻松搞定。 对于 Windows 用户，VMD 默认安装在 C :\Program Files(x86)\University of Illinois\VMD 文件夹下。打开 VMD。VMD打开后，会弹出三个窗口：VMD Main（主窗口），VMD Display（显示窗口），和命令窗口。 图4.14 VMD 主窗口、显示窗口和命令窗口 打开一个蛋白质结构（图4.15）：主窗口中点 File\New Molecule，弹出新窗口Molecule File Browser（文件读取窗口），文件读取窗口中点 Browse，自动打开 VMD安装目录，进入proteins文件夹选择VMD自带的演示结构bpti.pdb，文件读取窗口中点 Load，显示窗口出现蛋白质结构。 图4.15  打开一个蛋白质结构 默认的蛋白质结构显示方案（图4.16）：把显示窗口拖大后可以看到，VMD 读取了 PDB文件中的原子坐标，把每一个原子以细线的形式展示在 3D 空间中。不同的原子应的细线颜色不同。碳原子是青色的，氮原子是蓝色的，氧原子是红色的，氢原子是白色的，还有少量黄色的硫原子。 图4.16  默认的蛋白质结构显示方案 1、VMD 中鼠标的使用：把鼠标移到显示窗口里，按住左键，随意拖动，蛋白质就会在 3D 空间内任意旋转；按住鼠标右键拖动，蛋白质会在当前平面内 360 度旋转；前后滚动鼠标中键，可以将蛋白质放大或缩小。此外，鼠标还有更多的使用方法：主窗口上的 Mouse 菜单里可以切换鼠标模式。默认的鼠标模式是 Rotate Mode（旋转模式，R）。R模式下，鼠标在显示窗口内为单箭头。操作即为上述三种。将鼠标模式改为 Translate Mode（移动模式，T）后，鼠标在显示窗口内变为十字箭头 。T模式下，按住鼠标左键拖动为移动结构；右键和中键都为放缩功能。此外，还有Scale Mode（缩放模式，S）。S模式下，鼠标在显示窗口内变为左右箭头，按住鼠标左键或右键后，左右移动，即可完成连续的缩放。这种连续的缩放，不同于滚动中键实现的缩放。在默认的情况下，所有旋转是围绕整个结构的中心点进行的。我们也可以自己指定旋转的中心。这需要把鼠标模式切换为 Center（中心模式，C）。C模式下，鼠标在显示窗口内为十字。将十字放在想要作为旋转中心的原子上点一下，再按住鼠标左键旋转，就会以新定义的中心进行旋转。 2、恢复结构初始位置：主窗口中点Display\Reset View。或者在显示窗口内单击鼠标左键以激活窗口后，点击“=”键。 3、改变蛋白质结构的外观：有关外观的设置在主窗口中的 Graphic 菜单下的 Representation 窗口里（图4.17）。一个 Representation（显示状态）由三个元素构成。第一个元素是用什么样式（Style）显示，当前使用的样式是以细线显示原子（Lines）。第二个元素是用什么颜色（Color）显示，当前使用的颜色是按原子名定义的不同颜色（Name）。最后一个元素是要显示什么内容（Selection），当前显示的内容是所有原子（all）。这三个元素分别在 Representation窗口里的 Draw style 标签下的 Drawing Method 下拉条、Coloring Method 下拉条和Selected Atoms 输入框里进行设置。 图4.17 Representation 窗口 4、Drawing Method 下拉条（图4.18）：Lines 以细线显示原子。CPK 以不同大小的球来显示原子，原子间的连线是相应的化学键，比如碳与碳之间的共价键，硫与硫之间的二硫键等。NewCartoon 只显示蛋白质的碳骨架（backbone），并形象的展示出不同的二级结构。每一种 Drawing Method 都可以再进一步设置显示效果。比如对于 CPK，可以调整原子球的大小（Sphere Scale）、改变化学键的粗细（Bond Radius）、以及设置更高或更低的分辨率（Sphere/Bond Resolution）。 4.18  不同的几种Drawing Method 5、Coloring Method 下拉条（图4.19）：Name 颜色方案是一种原子一种颜色，常见的比如碳原子青色、氧原子红色、氮原子蓝色、硫原子黄色。在 NewCartoon 显示样式下，只有碳骨架被现实出来，再配以 Name 颜色方案，整个结构都是青色的。Secondary Structure 颜色方案可以为不同的二级结构赋予不同的颜色，常见的比如α螺旋紫色、β折片黄色、转角青色、松散 coil 结构白色。颜色方案 Res Type 根据氨基酸类型的不同赋予不同的颜色，比如非极性氨基酸白色、极性带正电荷的氨基酸（碱性的氨基酸）蓝色、极性带负电荷的氨基酸（酸性的氨基酸）红色、极性不带电荷的氨基酸绿色。这种颜色方案适合在显示氨基酸侧链的 Drawing Method 下观看，比如 CPK 样式。颜色方案 ResName 为20 种氨基酸设置了 20 种不同的颜色。 4.19  不同 Drawing Method 和 Coloring Method 组合 6、Selected Atoms 输入框：输入框里输入需要显示的内容，比如，写“all”代表显示所有原子，也就是整个蛋白质、写“backbone”代表显示碳骨架。输入框里允许输入的关键词和语法在 Selections 标签下有详细定义。可以写 Singlewords 里面单个的单词，也可以写 Keyword和 Value 组成的词组，还可以利用逻辑词“and/or/not”把多个单词和词组串成句子。举例： （1）点击 Reset 清空输入框里的内容→Keyword 里双击 ResName→Value 里双击 ALA（输入框里出现“ResName ALA”） →点击 Apply。显示名字为 Alanine 的氨基酸上的原子，即显示所有丙氨酸。配合 Drawing Method 设置为 CPK，Coloring Method 设置为 ResName。 （2）点击 Reset 清空输入框里的内容→Keyword 里双击 Resid→Value 里双击1→点击Apply。显示第一个氨基酸。利用这个 Keyword 和 Value 组合可以根据残基的编号选择某个或某一段氨基酸，比如，想要显示第1到第10个氨基酸，可以直接在输入框里输入“resid 1 to 10”，回车。此时显示的就是前10个氨基酸。 （3）Draw style 标签下，Drawing Method 设置为NewCartoon，Coloring Method 设置为 Secondary Structure→Selections 标签下，点击 Reset 清空输入框里的内容→Singlewords 里面双击 alpha_helix→点击“or” → Singlewords 里面双击 beta_sheet（输入框里出现“alpha_helix or beta_sheet”） →点击 Apply。通过逻辑词显示出所有的α helix 和β sheet。 VMD 可以设置多个representations（简称 rep），也就是将多个显示状态的试试效果叠加在一起（图4.20）。 7、representations（简称 rep）：设置第一个 rep：当前处于编辑状态下的 rep 背景色为浅绿色。设置 Drawing Method为 NewCartoon，Coloring Method 为 Secondary Structure，Selected Atoms为“all”，回车（或点击右下角 Apply 按钮）。创建第二个 rep：点击 Create Rep。点击后，复制产生了和第一个一摸一样的第二个 rep。浅绿色背景自动跳转到第二个 rep，即目前第二个 rep 处于可编辑状态。设置Drawing Method 为 cpk，Coloring Method 为 colorid，并选择选“1 red”（红色），Selected Atoms 为“resname PRO”，回车（或点击右下角 Apply 按钮）。设置两个 representations 后，蛋白质中的所有脯氨酸就以红色的原子球的形式叠加显示在以二级结构着色的碳骨架上了。 图4.20 设置两个 representations 创建第三个rep：点击 Create Rep。设置 Drawing Method 为 Surface，Coloring Method 为 colorid，并选择选“8 white”（白色），Selected Atoms 为“all”，回车（或点击右下角 Apply 按钮）。此时，可以看到填满肉的蛋白质，也就是蛋白质的外表面。继续设置 Coloring Method 右侧的 Material（材质）为 Transparent（透明材质）。如果显卡支持，可以打开 GLSL 显示模式：主窗口中点 DisplayRendermode选中 GLSL。GLSL 打开后的 Transparent Surface 变得更加柔美了（图4.21）。 图4.21  设置三个 representations 及 GLSL 显示模式 如果要删除某一个rep，比如删除第二个rep，需要先选中它（选中后背景色变为浅绿色），再点击 Delete Rep，就删掉了。或者你可以双击某一个 rep，比如双击第三个rep，将其暂时隐藏，等需要的时候再双击它取消隐藏（图4.22）。 图4.22  删除和隐藏 representation 8、保存与打开：当把蛋白质结构的显示效果调整到比较满意的状态后，可以保存当前所有的representations（注意保存的是显示状态，而不是结构）：主窗口中点击 File→Save Visualization State→保存在桌面上，起名叫 mystate.vmd。接下来关闭 VMD 再重新打开。这次我们不需要 Load 蛋白质结构，分别设置三个 rep，我们只需要直接载入刚刚保存的mystate.vmd ， 即可恢复刚刚的显示状态：主窗口中点击File→ LoadVisualization State→找到并打开mystate.vmd。这时，之前保存的显示状态就自动显示出来了（图4.23）。 图4.23 保存和载入显示状态 9、调换背景颜色（图4.24）：主窗口中点击 Graphics→Colors→弹出 Color Controls颜色控制窗口→Categories选Display→Names 选 Background→Colors选8 white。 图4.24  Color Controls 窗口调换颜色方案 Color Controls 窗口里可以调整各种VMD 默认的颜色方案。比如Coloring Method选Name的时候，默认的颜色方案是：氢原子白色，氧原子红色，碳原子青色等。这里，你可以根据需要把它们设置成其他颜色。 10、隐藏坐标轴：主窗口中点击 Display→Axes→off。 11、显示Lable：有时我们需要标记出某些氨基酸的序号或原子的名称（图4.25）。主窗口中点击 Mouse→Lable→Atom。把鼠标模式设置为Lable模式，让它标记原子。此时，鼠标变为十字。接下来在显示窗口中，要标记哪个原子就在哪里点一下。点击后出现文本显示的氨基酸名字、氨基酸序号以及被点击的原子的名字。 图4.25 显示Lable 12、调整 Lable： 1）改变字体颜色（图4.26）：主窗口中点击 Graphics→Colors→弹出 Color Controls颜色控制窗口→Categories 选 Lable→Names 选 Atoms→Colors选16 black。 图4.26 Color Controls 窗口改变 Lable 字体颜色 2）改变字体大小/粗细（图4.27）：主窗口中点击Graphics→Lables→弹出 Lable 窗口→Global Properties 标签下设置 Text Size 改变字体大小，Thickness 改变粗细。 3）改变显示位置和内容（图4.28）：主窗口中点击 Graphics→Lables→弹出 Lable 窗口→Properties 标签下选中要调整的 Lable→按住鼠标左键在 Offset 坐标系内移动来改变 Lable 的位置写入 Format 来改变 Lable 的内容。默认 Format 为：%R（氨基酸名字）%d（氨基酸序号）:%a（原子名字）。比如，只想显示氨基酸名字和编号，删掉“:%a”这部分代表原子信息的代码即可。   图4.27  Lable 窗口改变 Lable 字体大小/粗细 图4.28  Lable 窗口改变 Lable 位置和内容 13、保存结构图： 以刚刚bpti.pdb 这个蛋白质结构文件为例，设置了三个 representations，添加并调整了Lable，得到了最终的蛋白质3D 结构（图4.29）。最后要做的就是保存结构图片。屏幕截图是较为方便快捷的方法，但是所得图片分辨率不高。如果想要高质量的图片，可以用主窗口下的FileRender弹出 File Render Controls 窗口。File Render Controls 窗口里可以选择多种图片导出方式。 图4.29  bpti.pdb 的最终结构图片

（三）蛋白质三级结构的预测

蛋白质结构预测的方法有从头计算法，同源建模法，穿线法和综合法。

1、同源建模法（SWISS-MODEL）

原理：相似的氨基酸对于着相似的蛋白质结构，其步骤流程为：（1）找到与目标序列同源的已知结构作为模板（目标序列与模板序列间的一致度不能小于30%）；（2）为目标序列与模板序列（可以多条）创建序列比对。通常比对软件自动创建的序列比对还需要进一步人工校正；（3）根据第二步创建的序列比对，用同源建模软件预测结构模型；（4）评估模型质量并根据评估结果重复以上过程，直至模型质量合格。 SWISS-MODEL（www.swissmodel.expasy.org）是一款用同源建模法预测蛋白质三级结构的全自动软件（图4.30） 图4.30  SWISS-MODEL首页 点击Start Modeling 进入，输入序列。这里选择全自动，点击Build Model（图4.31） 图4.31 SWISS-MODEL输入页面 结果页面如图4.32， 图4.32 结果页面

2、穿线法

原理：不相似的氨基酸序列也可以对应着相似的蛋白质结构 穿线法首选I-TISSER（https://zhanglab.ccmb.med.umich.edu/I-TASSER/） 图4.32 I-TISSER首页   使用前必须先注册  

3、从头计算法（QUARK）

QUARK（https://zhanglab.ccmb.med.umich.edu/QUARK/）是一款用从头计算法预测蛋白质三级结构的在线软件，适用于没有同源模板的蛋白质，并且氨基酸序列长度200以内。

4、蛋白质三级结构预测方法的选择

**  **

5、蛋白质三级结构预测结果模型质量评估

  I-TASSER、Swiss-Model以及QUARK都自带评估体系，一定程度上能反应模型的质量，但不能完全反应，至少有3个评估体系完全认为模型可靠，模型才可用。这需要借助第三方评估软件，比如：SAVES（http://services.mbi.ucla.edu/SAVES/）平台提供的6种评估软件。无论用哪一种，只要有3种认为模型合格，就可以使用。 用ProQ（http://proq.bioinfo.se/cgi-bin/ProQ/ProQ.cgi）评估该模型质量 ModFold评估 ** ** ** ** ** ** ** ** ** **

（四）蛋白质四级结构分析

蛋白质四级结构的获取数据库如下图： ** **

1、蛋白质与蛋白质分子对接

** ** **

2、蛋白质相互作用分析软件

3、小分子化合物-蛋白质分子对接

** **

4、分子动力学模拟

** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** ** **