基因组组装及有参基因组间隙填补

   目前的基因组组装软件甚多，像常用的SPAdes、megahit、SOAPdenovo2、IDBA等。基因组的组装分为基因有参基因组组装和无参基因组组装（de novo assembl）在组装的计算量上无参基因组装计算量大、耗时久。而且对于不同的无参基因组组装每个软件的效果不同。例如对于细菌病毒的组装用SPAdes结果可能是最好的，而对于动物、植物的组装可能SPAdes又不如其他组装软件。

        以无参基因组组装为例，要获得良好的组装效果首先在于测序数据的质量。一般来说最好情况是测序深度足够深、样本量足够大保证建库的随机性。覆盖率足够高，这样才是获得良好组装结果的先决条件。组装基因组目前是基于contigs、scaffold、Chromosomes三个层次。组装过程是先将reads根据重叠区组装成contigs再将contigs连接组成scaffold，最后根据遗传图谱将scaffold组成Chromosomes。尽管目前基因组装软件组装效果进步很大，但是还是无法避免gap。因为组装基因组的三个原则本身就是矛盾的。这三个原则是

1. 尽可能得到长的contigs保证连续性；

2. 尽可能组装得到错误率最低的contigs;

3. 尽可能包含整个原始序列；

显然，要得到足够长的contigs和错误率最低的contigs就会舍弃很多reads,无法满足第三个原则，而要尽可能包含所有reads就无法满足第一二条原则。所以对于gap的填补，可以通过做PCR进行填补。但是如果gap很多，那么大量地做PCR也是非常辛苦的。或者如果实验室经费足够充足，大量做测序也是可以的。对于病毒或者细菌的无参基因组组装选择SPAdes是个不错的选择，使用方法为

spades.py –meta -m 80 -t 30 -k 21, 33, 59, 79, 99, 119, 127 -1 R1.fastq -2 R2.fastq -o out

spades.py –sc –careful 80 -t 30 -k 21, 33, 59, 79, 99, 119, 127 -1 R1.fastq -2 R2.fastq -o out

上面两个代码分别是对于病毒或细菌宏基因组的组装和单细胞测序的组装，当然也可以不选k-mer值，因为SPAdes会根据我们输入的fastq数据自动选择合适的k值。而在–sc情况下选择–careful参数可以启动bwa进行错误校正。

也可以使用megahit，megahit运算速度快，组装效果也好

megahit -m 0.5 -t 30 -1 R1.fastq -2 R2.fastq –min-contig-len 80 -out out

megahit中–min-contig-len参数是对最小contig的长度设定，设定稍小的值有利于获得更长的contigs连接。

        对于有参基因组的组装，虽然计算上简单了很多，但是考虑个体的变异性，也不能完全依靠refsequence。尤其是对于病毒这类突变很多基因组组装，完全依靠参考基因组将会丢失很多变异信息。为了保留变异信息，我想到了通过blast结果m6格式信息，进行有参基因组组装的间隙填补。blast输出的m6格式文件中包含query序列（自己的序列）和subject序列（参考序列）比对的起始、结束信息。对于基因组的组装，我们不考虑某些基因的表达量，仅考虑如何获得最多的非冗余基因片段组装为基因组。而大量的contigs有长有短，组装这些序列总体会遇到三种情况，1.某一条序列对于与其他序列并无重叠区，2.某一条序列有一部分与其他序列重叠，3.某一条序列完全重叠在其他序列上。对于这三种情况，只有1，和2的不重叠区段才是我们想要的。所以我们只需要找到这些区段的位置信息即可。结合blast文件的位置信息，我们需要通过对subject和query序列的起始结束位置信息和重叠长度进行计算，算出不重叠区段的位置信息即可。（m6格式文件中第7列是自己序列比对起始位置，第8列是自己序列比对的结束位置，第9列是参考序列比对的起始位置，第10列是参考序列比对的结束位置）由于这次的代码很多，我用图片方式post出我的代码：

在服务器中进行本地blast操作，那部分属于Shell命令。

上面代码的核心思想在于将contigs按照从长到短的顺序排序，提取出比对相似度大于90%的contigs，将contigs的非冗余比对部分计算并抓取，最后根据blast的信息校正序列的方向。具体的算法自己根据代码进行理解。