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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 215,012评论 6 497
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,628评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 160,653评论 0 350
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,485评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,574评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,590评论 1 293
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,596评论 3 414
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,340评论 0 270
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,794评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,102评论 2 330
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,276评论 1 344
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,940评论 5 339
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,583评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,201评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,441评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,173评论 2 366
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,136评论 2 352

推荐阅读更多精彩内容