序列拼接简介

丁香园论坛：/bbs/thread/1247063#1247063
问：从测序仪的光密度采样与分析、碱基读出、载体标识与去除、拼接、填补序列间隙、到重复序列标识、读框预测和基因标注的每一步都是紧密依赖基因组信息学的软件和数据库。

答：
一、这应该是对DNA及mRNA的基本分析，有很多免费的软件可以利用，但是想做流程就需要用perl一样的胶水语言进行组合。

从测序仪结果开始：
phred 进行base calling，即碱基判读
cross_match 去除载体及引物序列
repeatmask屏蔽重复序列
longorf.pl 进行最长读码框预测
blast2/blat定位样本序列到已知基因或者基因组
用emboss软件包中各种软件可以进行进一步分析，如DNA/RNA/PRotein二级结构，跨膜区，信号肽分析等。

GO数据库对基因分类研究
还可以进行分支研究，如利用测序结果进行SNP/Mutation研究，用polyphred/consed，或者mutation surveyor软件（有demo版和商业版）
再以后的研究思路就非常细化了，可以结合具体分析目的进行。

二、1、基因组序列拼接——phred/phrap/consed
Phred 简介
Phred是一个采用快速傅利叶变换分析技术以及动态规划算法从DNA测序所得到的图形数据中提取DNA序列排列顺序信息（Base Calling）得到DNA序列的软件。

Phred 对序列中的每一个数据产生一个被广泛接受的带有质量控制标准（quality scores）的“Base Call”。

Phred质量指标x就相当于约10-x/10的误差概率。

因此，PHRED质量指标20就相当于在原始数据中一个Base Call的精确度为99%。

Phred可以读取DNA测序仪生成的色谱图文件（二进制格式），通过分析每个碱基的“质量”信息而输出每个测序序列的碱基序列和质量信息文件（文本格式）。

它自动的判断并读取ABI 373、377、3700和MegaBase等大多数DNA测序仪产生的色谱图文件，而且还可以自动识别经过gzip或Unix compress压缩的数据文件。

在完成Base calling以后，Phred可以输出FASTA格式或者是XBAP、PHD和SCF格式的文件。

序列的质量信息被写入FASTA或是PHD的输出文件中。

而这些文件可以在下一步的组装中被Phrap等其他程序所使用。

Phrap简介
Phrap是一个用于将鸟枪法测序的原始序列拼接成连接群（Contig）的软件。

这个软件的核
心算法是Smith-Waterman动态规划算法。

Phrap结合相应的质量控制标准值对Phred Base Calling所得到的DNA序列进行拼接。

在拼接过程中Phrap可以自动地在Phred质量指标的基础上计算拼接后得到的连接群序列中的每个位置上数据的误差的概率，并给出其质量控制标准值（quality scores）。

Phrap质量指标x就相当于约10-x/10的误差概率。

因此，PHRAP 质量指标20就相当于在装配起的序列中的一个碱基的精确度为99%。

Consed简介
Consed是一个基于Unix/Linux系统的专门用来查看和编辑phrap的拼接结果的可视化工具，目前的最高版本是13.0，支持的操作系统包括Solaris、Redhat Linux、Irix、DEC-Alpha和HP-UX等。

由于Consed软件也是华盛顿大学（University of Washington）开发的，所以它从开发阶段起就紧密的与phrap结合——Phil Green（Phrap的作者）直接指导了它的设计。

因此，只有使用Consed软件才能充分发挥phrap拼接的所有功能（比如有效的读取phrap 所产生的一系列标记等）。

而且Consed的图形界面使得phrap的拼接结果变得很直观，用户可以直接查看测序色谱图，并可以很容易的检查一些低质量区域、错拼或由于Base-Calling 的误差而造成的不匹配等等，甚至还可以直接根据叠连群的序列设计PCR引物等等。

另外，在后期的Consed版本中还加入了一个新的称为Autofinish的功能，可以很大程度的简化补洞的工作量并降低Finishing阶段的成本。

关于重复序列标识
RepeatMasker是专门用于查找并遮蔽重复序列区域的程序。

其核心程序其实是phrap包中的cross_match程序，而它本身只是一个perl脚本，另外加上一些自带的通用重复序列库。

因此，对于具体项目而言，应该自己添加针对自身项目的一些特有的重复序列库（其实就是一个标准的Fasta格式的文件）。

另外，consed包中自带有一个tagRepeat.perl程序，做法与RepeatMasker类似，也是用cross_match来标记重复序列区域。

但是它可以将其转为consed能够识别的标签，在consed 的图形窗口中直接显示，因此在finishing过程中会非常有用。

【注】上面的这一系列程序都是只能在Unix/Linux系统中运行的！
2、基因组序列基本注释
ORF识别
以前曾经写过，请参见过去的帖子：对于细菌等原核生物的ORF的识别，目前理论手段已经发展的比较成熟，而且特异性（specificity）和灵敏度（sensitivity）都相对较高。

其中使用最为广泛的软件是由TIGR中心开发的Glimmer，这个名称是“基因定位和内插马尔可夫模型”（Gene Locator and Interpolated Markov Modeler）的缩写，它是专门用于在微生物基因组中定位基因的软件，对于细菌（bacteria）和古细菌（archaea）特别有效。

Glimmer 是对学术研究机构免费的，详细的信息和如何获取可以从TIGR的网站上的相关网页中了解
到（/software/glimmer/）。

目前可以得到的是它的改进版本2.10的源代码包，解开编译（make）以后将产生四个独立的可执行程序：①long-orfs，用于根据启动子和终止子的位置，给出可能的范围较长且重叠较少的ORF；②extract，用于根据指定的信息从序列文件（FASTA格式）中提取特定的子序列；③build-icm，用于建立并输出内插马尔可夫模型（IMMs）；④glimmer2，用于根据给定的序列文件和内插马尔可夫模型（由build-icm程序产生的）列出所有可能的ORF。

具体的参数和用法可以参阅其软件包中提供的针对每个程序的名为“*.readme”的文档。

在实际使用时可以直接运行其提供的简单的shell脚本run-glimmer2，也可以根据各自的需要调节参数或单独运行。

另外，对于真核生物的ORF预测，虽然总体的效果不如原核生物，但是目前也已经开发出了不少现成的软件工具。

比较常用的是由斯坦福大学开发（现在由麻省理工学院生物系的Burge实验室维护）的Genscan软件，它是一个“普适”的基因预测程序，它可以用于分析许多物种的基因组序列，包括人、脊椎动物、无脊椎动物和植物等。

它主要是基于给定物种的基因结构和组成特性的一个统计模型来预测出最有可能的基因，并可以给出文本和图形（PostScript格式）两种结果形式。

Genscan所使用的这个统计模型中考虑到了许多主要的基因序列的结构特性，例如典型的基因密度、一般每个基因中的外显子数目、不同类型外显子的大小分布等。

并且也考虑到了许多重要的基因组成特性，例如翻译起始位置的特殊组成（Kozak框）以及终止信号、TATA框、加帽位点和poly A信号等。

麻省理工学院提供了使用Genscan的WEB服务器（/GENSCAN.html）和Email服务器（/GENSCANM.html）。

同时Genesan对学术机构也是可以免费获得的，其主页上提供了针对多种Unix/Linux系统的版本，而且本地使用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”后缀），其中显示了预测出的基因和外显子的分布情况。

或者/gorf/gorf.html
功能预测
最常用的手段是使用由NCBI开发的BLAST软件包中所提供的一系列用于序列比对的程序，将在新基因组中预测出的所有的基因在公共数据库，比如GenBank所提供的非冗余蛋白数据库（简称NR），进行检索。

对于这样大批量（通常在103数量级以上）的检索，完全通过互联网上所提供的Web服务显然是烦琐而且低效的，所以建立一个本地的（local）蛋白数据库和BLAST软件工具是非常必要的，它无疑将大大提高注释过程的效率。

NR蛋白质数据库可以从NCBI的Ftp服务器上直接下载（ftp:///blast/db/nr.Z），而BLAST软件包则需要根据所要使用操作系统而选择的相应的编译过的可执行程序下载（ftp:///blast/executables/）。

关于BLAST软件包的使用，可能很多人已经比较熟悉了，而且过去的很多专门介绍生物信息学的书中都有所涉及，因此这里不再做更多的讨论了。

如果不清楚，可以参阅过去的帖子：Blast 简介
当BLAST程序无法找到已知功能的同源基因时，还可以进行更进一步的预测。

这里最为常用的手段是使用华盛顿大学（Washington University）开发的HMMER软件包所提供的工具对Pfam数据库进行检索。

Pfam数据库（protein families database of alignments and HMMs）是由英国的Sanger中心在1998年建立并维护至今的，其中收集了大量通过多序列比对产生的在一定范围具有隐马尔可夫模型（Hidden Markov Models，简称HMMs）的蛋白质结构域和蛋白质家族。

我们知道，在进行多序列比对时，由于有多个亲源关系不等的序列包括在内，因此需要插入一些空位来使比对的序列形成正确的匹配，随着空位的插入，一些具有保守性的区域便形成了。

这些保守区域通常由10-20个氨基酸的长度，并对应着蛋白质的核心结构或功能区域，因此它们的特征就可以用来对蛋白质家族的成员进行鉴别。

HMMER软件包中包含了一系列基于隐马尔可夫模型的应用程序（也包含一些其他的有用的但非HMMs的程序），它是一个完全开发源码的软件，目前最新的版本是2.3.2，可以从它的主页直接下载其源代码包或者选择相应系统的可执行包。

其中的hmmbuild程序可用于处理多序列比对来产生具隐马尔可夫模型的蛋白质结构域，而hmmpfam程序则可以检索基于隐马尔可夫模型（HMMs）的数据库（如Pfam等），因为它采用了位置特异的计分方法（position specific score），来处理多序列比对中不同位置上不同保守程度的氨基酸以及不同程度的插入或缺失，而不象传统的序列比对算法,如FASTA和Smith-Waterman算法等，用的都是位置无关（position-independent）的计分方法。

如果使用上述的两种方法都无法确定出一个预测出基因的功能，那么它有可能就是一个新的未知功能的基因，通常可以将它们的产物（product）描述成“hypothetical protein”，这样可以方便以后的功能基因组或蛋白质组的研究。

RNA识别
对于rRNA，由于它的序列极端的保守，而且过去由于研究进化关系的需要，许多物种的rRNA 序列已经被测定了。

所以一个比较简单的方法就是找到这个物种的rRNA的序列（如果它是
已经测定的）或者选取一个已知的与其亲缘关系最相近的其他物种的rRNA的序列，用BLASTN程序将它们在基因组序列中进行检索即可。

而对于tRNA，由于其特殊的结构和序列特征，目前已经有许多软件工具可以对它进行识别。

最常用的软件是由华盛顿大学（Washington University）开发的专门用于识别基因组DNA 或RNA序列上的tRNA基因的tRNAscan-SE。

它结合了Cove概率法RNA预测软件包的特异性和tRNAscan1.3的速度和灵敏度，同时又新加入了一种在此程序中称为EufindtRNA的算法（此算法最初是由Pavesi等人在1994年用于查找真核生物体中的Pol III tRNA启动子的）。

tRNAscan-SE软件在处理过程中首先使用tRNAscan和EufindtRNA算法来进行“预查找”，从而识别出一部分候选的tRNA区域序列，然后将这些子序列送给Cove算法再做进一步的分析，最终会将被Cove确认为正确的预测结果输出。

用这种方法，tRNAscan-SE程序取得了三个方面的优势：①假阳性（false positive）的比率很低——每预测150亿个碱基不超过一个；②很好的结合了tRNAscan和EufindtRNA两种算法的灵敏度（可以识别出99％的真实的tRNA）；③运行速度比Cove程序快1000到3000倍，比原来的tRNAscan1.3程序也快30到90倍（因为tRNAscan-SE对原来的tRNAscan1.3和Pavesi的算法都进行了代码上的优化）。

3、更进一步的生物信息学分析
上面只是最最基本的简单注释的内容，接下来的分析还有很多，例如基因组组分分析（composition analysis）、基因功能分类（biological role category）、比较基因组研究（comparative genomics）和物种进化分析（phylogenetic analysis）等等。

这里涉及到的数据库和软件就太多了，如EMBOSS、BioEdit、KEGG、WIT、BioCyc、BioCarta、GO、COG、ACT、MUMmer、Alfresco、VISTA、Pipmaker、GenomeComp、phylip、PAUP……有些我也不太熟悉，不敢妄加评论。

另外，对于基因组数据的处理，通常由于其大规模的分析要求，很多软件工具未必能解决问题。

而且，各种软件之间的输入和输出的格式过滤、处理等等的问题都是必须自己设法解决的。

所以我认为，分析人员最好具备一定的编程能力，这里最适合的语言首推Perl。

强大的文本处理能力和在Unix/Linux系统上的通用性都是其它语言无法匹敌的。

（呵呵，个人意见，仅供参考）。