基因组组装

1.k-mer

那么我们首先要看一下k-mer是什么。
它的定义是：是指将一条序列分成包含k个碱基的子字符串，如果reads长度为L，k-mer长度设为k，则产生的k-mers数目为：L-k+1，例如序列AACTGACT，设置k为3，则可以将其分割为AAC ACT CTG TGA GAC ACT共6个k-mers。其中k一定是奇数，如果是偶数遇到回文序列可能会产生完全相同的k-mers。
这是从网络上抄下来的定义，参考：
https://www.jianshu.com/p/d434abfa7eaa

2.与k-mer相关的公式及意义

其中
k-mer（number）= Read_length - k-mer_size + 1
如果假设基因组大小为G，由于测序是随机的，read长度为L，n（base）和n（kmer）分别为测序的碱基和kmer个数，C（base）和C（kmer）分别为碱基和kmer的覆盖深度
那么有：

一般来说，k-mer的深度比碱基深度要低一些

3.k-mer杂合度

对于杂合基因组，来说k-mer可分为杂合k-mer和纯合k-mer，其中k为k-mer的个数，那么对有些位点来说，有2k个杂合k-mer覆盖，如图红色部分：

以此，我们定义α1/2为杂合k-mer种类数的百分比，n（kspecies）为所有k-mer种类数，那么杂合度计算如下：

4.k-mer分析常用软件

k-mer分析的常用软件有Jellyfish，GCE，GenomeScope
这样我们就可以看到一些kmer的分布，基因组大小和重复序列等一些指标

5.基因组组装技术

从测序的代数看，
一代测序双脱氧终止法：成本高，读长长，但是通量低
二代测序焦磷酸（illumina）：成本低，读长短，通量高且拼装的原理是基于德布鲁因图来拼装的
三代PacBio，Nanopore：成本高，读长超长，通量高，基于overlap拼接
通常来说拼接顺序如下：

reads：为我们测序的读段reads
contig：片段重叠群，指拼接软件基于短序列之间的重叠区（overlap），拼接获得的较长序列。
片段框架，由先后顺序已知的Contigs组成的序列，中间有Gap。基因组de novo测序中，通过reads拼接获得Contigs后，往往还需要构建454 Paired end库或Illumina Matepair库，以获得一定大小片段（如3Kb、6Kb、10Kb、20Kb）两端的序列。基于这些序列，可以确定一些Contig之间的顺序关系，这些先后顺序已知的Contigs组成Scaffold。
N50：Reads拼接后会获得一些不同长度的Contigs。将所有的Contigs长度相加，获得Contig总长度。然后将所有的Contigs按照从长到短进行排序，如获得Contig1，Contig2，Contig3…Contig25。将Contigs按照这个顺序依次相加，当相加的长度达到Contig总长度的一半时，最后一个加上的Contig长度即为Contig N50，可以作为基因组拼接的结果好坏的一个判断标准。此概念很容易被误认为所有Contigs长度排名第50的序列长度，与之类似的有N90，N50与N90同样适用于Scaffolds。

6.二代数据组装基因组常用算法和软件

目前主流算法有：

1.OLC

Over lap Layout Consensus

基于OLC算法的拼接，基本原理如下：
首先找到reads与reads之间的overlap区域，然后根据这些overlap区域拼接出scaffolds，最后基于这些scaffold的关系确定先后顺序，完成拼接

2.DBG

这个方法主要基于图论里面的德布鲁因图来完成拼接的

比方说以读段AGATACT为例，那么3-mer（k=3）构成的德布鲁因图如上，每个3-mer前后错位1个碱基，依次排列开。
那么对于两条读段来说，首先找到两条序列重叠的部分，发现它们之间仅一个碱基不同，根据德布鲁因图，发现它们其中三个k-mer会有所不同，并形成一个小包

那么小包外的部分就确定为我们的基因组的一部分了
对于10个读段来说：

对于不同的k-mer，会隆起一个小包，那么小包外的部分即可以确定为基因组的一部分，对于小包来说，根据更多的读段，在进行打分，从而选择最优的那条路径

目前二代组装常用的软件有：
1.OLC：Canu，Celera Assembler，Falcon等
2.DBG：SOAPdenovo2，SPAdes，ALLPath-LG，Velvet等

7.三代数据组装基因组

三代数据读长较长，需要30X以上的深度
它的拼接基本想法和二代的OLC比较类似，先寻找测序reads之间的overlap关系，把它们拼接成Contig，然后根据图论的方法，解出最优路径，然后对Contig进行矫正，从而拿到一致性序列

只不过较二代数据来说，由于三代读长较长，所以基于长reads拼接较为准确直观
目前三代拼接的软件有：
1.CANU
2.WTDBG2

8.遗传图谱挂载染色体

一般来说，无论是二代还是三代数据，按照上面两节介绍的技术来组装，只能组装到scaffolds的水平，接下来拼接到染色体水平，需要利用HI-C来做染色体挂载

Hi-C建库

Hi-C建库原理是，在同一条染色体上的DNA片段因为空间距离更近，所以更容易被蛋白交联在一起。所以Hi-C测的是被蛋白交联的DNA片段，所以这样组装就可以将这些Contig依据染色体进行分组

像热图红色部分即为染色体分组
目前染色体挂载常用的软件有：ALLHIC，3d-DNA，LACHESIS，SALSA2

基因组注释

1.基因组序列的结构组成注释

说到基因组组装，我们就应该先了解下基因组序列的组成。
对于单倍体生物来说，每个细胞都应该都有一套完整的基因组，那么该基因组包括了编码序列和非编码序列，那么对于有性生殖生物，基因组包括核基因组，叶绿体基因组和线粒体基因组，不过这里我们重点讨论核基因组

对于我们的基因组来说，大部分都是重复序列和非编码序列，只有5%是蛋白质编码序列，这些蛋白质编码区域散落在基因组的各个部位
那么对于蛋白质编码区域，其结构为：

目前对基因组结构预测（注释）的方法主要有两种：
1.同源比对预测
2.从头预测

（1）同源比对方法

这个方法的实现要基于同源物种的基因组是被注释过的。该方法是利用近缘物种已知的基因组进行比对，从而发现同源序列，对于近缘物种来说，我们认为它们的同源序列具有相似的结构特性和功能特性

那么发现同源序列后，根据一些分子（碱基）特征，包括翻译起始位点，终止位点，内含子和外显子部位等进行注释
或者根据转录本序列，结合表达序列标签（EST序列）来进行注释

（2）从头预测法

从头预测主要是根据编码区的统计信号以及序列信号来进行预测，主要有HMM模型和贝叶斯理论

一.HMM模型

对于HMM模型来说，最早是用于预测元和生物基因组结构的
它的基本原理是对于一段序列来说，构造出几个状态。以最简单的模型为例：
外显子态，内含子态和它们之间的5'端剪切位点态

那么对于该序列的每一个碱基，我们都可以看做是在这三个状态中的其中一个状态，那么在往下一个碱基转移的时候，可以看做是三个状态相互转移的过程
那么给出转移概率表：

这个表的意思是，确定第一个碱基状态是外显子的概率为1，在此状态基础上为C碱基的概率为1/4，在此基础上，第二个状态还是外显子态的概率为9/10，在此状态基础上为T碱基的概率为1/4，以此类推，计算出所有有可能的路径得分
例如：
第一条路径：
第一条路径得分

最后一条路径：

最后一条路径得分

可参考：https://www.jianshu.com/p/866bfb75586a

至于这样的转移概率和发射概率是如何得到的，通常是基于先验知识统计而得到的，一般来说，A或G为剪切位点的概率较大，所以出现A或G位点均有可能是可变剪切位点

二.贝叶斯模型

贝叶斯模型是HMM模型的升级版，它专门用于解决路径得分相近的情况，基于上述例子，19位的G得分为-41.22，23位的G的得分是-41.71，由于得分相近，那我们要如何判断哪一个剪切位点才是最优的情况？
首先，根据贝叶斯理论：
已知p(x | π)，反推p( π | x )；即已知在某种状态的条件下，出现某个碱基的概率为p(x | π)，那么反推在已知某个碱基的条件下，为某个状态的概率

通常我们要求的p(πi)是为剪切位点态的概率，因为我们的目的是判断剪切位点，所以p(πi)为确定的事件，故概率为1。
并以此计算这两个G位点的后验概率，比较大小即可得到最优可变剪切位点

2. 基因结构注释软件

刚才我们介绍了基因结构预测的原理——HMM模型
那么基因结构预测层次分为DNA水平，RNA水平，蛋白质水平

一.DNA水平

正如我们上面介绍的HMM模型，根据已经研究好的物种基因区域的特征，那么我们抓住这些特征构建HMM模型，那么用这个模型去预测未知物种的基因结构，但是这样做最好是用近缘物种训练模型最好
对于找不到近缘物种的物种怎么做鉴定呢？那我们可以找做实验克隆出来的该物种的基因序列做训练集训练模型

二.RNA水平

利用三代全长转录本比对到参考基因组上，那么全长转录本比对上的部分就很可能是外显子区域（blastn），或者用二代测序技术进行拼接（Trinity进行转录本拼接）然后再比对到参考基因组

三.蛋白质水平

我们寻找的蛋白质的数据源，大部分都是来自于近缘物种，利用近缘物种来进行预测，利用blast比对到参考基因组上（或者利用exonerate比对）

四.可信度检验

基于RNA和蛋白水平的预测后，我们还要整合所有证据进行可信度分析
整合证据常用Maker，EVM等，输出文件是GFF/GTF，cds序列，protein序列

五.参数训练

我们利用HMM模型预测新物种基因结构，我们需要利用近缘物种训练模型参数，那么我们利用RNA水平的证据训练模型参数，通常Braker就可以训练，当然也有已经训练好的模型参数可以直接用

3.基因功能注释

一.利用序列相似性进行功能注释

基因功能注释可以通过NR，Uniprot等数据库进行注释，它的基本原理是利用blast进行搜索，获得最佳匹配的基因片段，根据数据库里面的对应片段的功能，来推测未知基因的功能
或者是用Interpro数据库进行注释，输入DNA序列进行远程比对，最后返回注释结果

二.利用GO，KEGG信息进行注释

对于GO注释信息，我们可以利用WEGO（http://wego.genomics.org.cn/
），AgriGO（http://systemsbiology.cau.edu.cn/agriGOv2/
），GOEAST（http://omicslab.genetics.ac.cn/GOEAST/index.php
）等在线分析，但这些大部分都是根据gene ID来做的注释，如果是一段未知序列，我们通常采用blast2go来获得注释（比方说NR库）
对于KEGG，可以利用官网提供的KAAS来做分析，这个功能可以利用序列相似性来进行功能注释

三.功能注释

（1）同源注释

同源注释的基本原理是基于序列相似性（软件：blast/Diamond，数据库：NR，Uniprot）或者基于结构域（软件：HMMER/Interproscan，数据库：Pfam/Interpro）
如果是少量蛋白质注释直接在NCBI或者Uniprot上直接在线注释
如果是大量的蛋白质注释，可以在本地批量注释

4.重复序列注释

基因组上有很多重复序列，可分为：
1.简单重复序列
2.散在重复序列（转座子），左右为重复序列，中间有一定的功能域；转座子又分为DNA转座子和RNA转座子（需逆转录）
3.WGD，全基因组复制，后期丢失一些序列
4.基因的多拷贝，基因家族
5.染色体大片段复制

那么重复序列注释的原理是什么呢？一般重复序列我们可以分为已知重复序列和未知重复序列
那么对于已知重复序列，我们通常在Dfam和Reobase数据库上做blast来做注释
那么对于未知重复序列，我们一般分为简单重复序列（用trf做鉴定），重复次数较高的序列（RECON，RepeatScout做鉴定，这两款软件原理是基于k-mer）然后与近缘物种数据库进行比对，做家族分类，LTR（LtrHarvest，Ltr_retriever做鉴定）

注意在注释重复序列，首先构建RMBLAST数据库，然后再生成某物种自身的重复序列数据库，再鉴定重复序列并可视化

参考：部分参考《生物信息学》樊龙江