Jellyfish进行Kmer评估基因组大小

        我们对一个新物种进行基因组测定时，一般需要先对该物种进行基因组survey，估计基因组大小、杂合率、重复序列比例、GC含量和倍性等信息，其中基因组大小直接影响我们后续三代测序的数据量。公司一般会先用二代测个100x的数据量来做这一步，现在二代已经很便宜了，测50X和100X也差不了多少钱(当然也不是越多越好)，但在测序之前，建议通过核型分析和流式等方法预估一下基因组大小是多少再安排送测序。综上，通过survey，我们能对所测物种的基因组的基本信息有一个初步了解，得到高质量的二代clean data 为后面的三代数据进行纠错。

        在得到下机数据*.fq.gz后，先通过jellyfish计算k-mer分布频率，然后通过Genomescope2估计基因组大小、重复和杂合度。首先，我们需要先知道什么叫K-mer

一K-mer 分析原理：

        从一段连续序列中迭代地选取长度为K个碱基的序列，若每条序列的长度为L，则K-mer长度为K，那么可以得到L-K+1个K-mer，后面我取K=21来进行分析(针对长度为K的K-mer，对于A,T,C,G四种碱基类型，一共能产生的K-mer种类数为4的K次幂个，选择K=21是为了产生足够多的K-mer种类数去覆盖整个基因组，当然，你也可以设置不同的K值进行探索，但21是比较常用的)。

         此外，针对一些物种可能会有较高的杂合率及多倍性，K-mer分布图会有不同的特征，比如：杂合K-mer，重复K-mer，需要我们进一步了解掌握。

对于不同的基因组杂合度，Kmer分布如下：

简单基因组：

        图中有一个明显的主峰，没有其它峰值，根据Kmer分布图，可初步判断该物种基因组为简单基因组。

高重复基因组

        图中在49x深度和104x深度位置分别有1个峰，其中49x的位置为主峰即基因组正常期望深度，104x位置为重复峰位置。根据kmer 分布图，可初步判断该物种基因组为高重复基因组。

高杂合基因组

基因组大小预估

我们需要了解K-mer深度分布曲线含义：K-mer深度—K-mer深度频率。从所有测序Reads中逐碱基取K-mer，并统计K-mer深度频数分布，从而获得K-mer深度分布曲线。根据曲线获得K-mer深度估计值，用于估计基因组大小。

        基因组大小=(基因切割的Kmer数目)/(主峰深度)

        K-mer分析适用于分析唯一主峰区域所占比例较大的基因组，当基因组杂合非常高或者重复序列比例非常大时，其影响可能导致无法通过K-mer分析正确估计基因组大小。将K-mer深度等于1的情况认为是错误情况，计算错误率，并用于修正基因组大小。

二、K-mer实战

2.1 软件安装：Jellyfish，

        Jellyfish是CBCB(Center for Bioinformatics and Computational Biology)的Guillaume Marçais 和Carl Kingsford 研发的一款计数DNA 的 k-mers 的软件。该软件运用 Hash 表来存储数据，同时能多线程运行，速度快，内存消耗小。该软件只能运行在64位的Linux系统下。其文章于2011年发表在杂志Bioinformatics 上。在其官网上有pdf版本使用手册，具有详细介绍。jellyfish的功能有：kmer计数；融合二进制的Hash结果；统计Hash结果；通过Hash结果来画直方图；将Hash结果输出成文本格式；查询指定k-mer的数目。

        安装jellyfish，按照官网介绍逐步安装就好了，非常好安。

        $ wget http://www.cbcb.umd.edu/software/jellyfish/jellyfish-2.3.0.tar.gz

        $ tar zxvf jellyfish-2.3.0.tar.gz

        $ mkdir jellyfish

        $ cd jellyfish-2.3.0

        $ ./configure --prefix=$PWD/../jellyfish

        $ make -j 8

        $ make install

        也可以在github上下载，https://github.com/gmarcais/Jellyfish，里面有更多的详细介绍。

2.2 Jellyfish count

        我们要用的是count 命令来进行K-mer计数，使用fastq文件在默认参数上和fasta文件没有区别，生成的hash结果为*.jf为后缀的二进制文件。

        以count基本命令为例：

        jellyfish count -m 21 -s 100M -t 10 -C reads.fasta

        参数含义可以通过jellyfish count --help来进行理解，(大部分参数默认就可以)：

        其中我觉得可以根据自己需求进行设置的有：

        -m 使用的k-mer的长度。如果基因组大小为G，则k-mer长度选择为: k ~= log(200G)/log(4)；

        -s Hash的大小，最好设置的值大于总的独特的(distinct)k-mer数,这样生成的文件只有一个。若该值不够大，则会生成多个hash文件，以数字区分文件名。如果基因组大小为G，每个reads有一个错误，总共有n条reads，则该值可以设置为[(G + n)/0.8]。该值识别M 和 G；

        -t |--threads=  default: 1  使用的CPU线程数；

        -o |--output=  default:mer_counts 输出的结果文件前缀；

-c |--counter-len=  default:7， k-mer的计数结果所占的比特数,默认支持的最大数字是2^7=128。对于基因组测序覆盖度为N，则要使设置的该值要大于N。该值越大，消耗内存越大；

       值得注意的是，我们的下机数据是*fq.gz压缩文件，这个时候通过zcat命令及管道符“|”进行计算，一定要保留/dev/fd/0 , 这个是进程输入标志，也就是代表管道符前边的结果传递，其它参照上述参数设置。

2.3. stats命令对hash结果进行统计

           k-mer的结果以hash的二进制文件结果给出，需要统计出k-mer总数，特异的k-mer数目，只出现过一次的kmer数，出现了最多的k-mer的数目等信息。

        $ jellyfish statshash

        示例结果为：

        Unique:    10623990377    #只出现过一次的k-mer的数目

        Distinct:    15385280271   #特异性的k-mer数目，包含上一个的数据

        Total:      297422026708  #总的k-mer数目

        Max_count: 108593759      #同一个k-mer出现的最多的数目

2.4 histo输出k-mer频率直方图数据，如：

        对k-mer的计数结果有个直观的认识，则需要统计出现了x(x=1,2,3…)次的kmer的数目y,以x，y为横纵坐标画出直方图。使用 histo 命令能给出 x 和 y 对应的值，将结果默认输出到标准输出。

得到*.histo后缀的文件后，就可以通过Genomescope2进行K-mer分布结果的可视化作图了，并得到预估的基因组大小、杂合率、重复序列比例。

三、附：

        普通基因组：单倍体、纯合二倍体或者杂合度<0.5%，且重复序列含量<50%，GC含量为35%到65%之间的二倍体。

        复杂基因组：杂合度>0.5%，重复序列含量>50%，多倍体，GC含量处于异常的范围（GC含量<35%或者GC含量>65%的二倍体）。

        二倍体复杂基因组进一步细分为:微杂合基因组（0.5%<杂合率<=0.8%)；高杂合基因组（杂合率>0.8%)；高重复基因组（重复序列比例>50%）。

参考资料：

https://github.com/gmarcais/Jellyfish/

http://www.chenlianfu.com/?p=806

https://www.jianshu.com/p/94da86093843

http://www.360doc.com/content/21/0714/12/76149697_986499807.shtml