k-mer分析:你的基因组有没有被污染?

2022-05-05 13:30:33 浏览数 (2)

k-mer分析是指通过k-mers深度(也即k-mers出现次数)的分布规律(一般通过分布曲线或直方图展示)来估计基因组的一些基本信息,例如基因组大小、杂合度、纯度等,同时也可以判断组装时的最佳k-mer值,是二代测序基因组组装前的准备步骤。

k-mer分析常用的软件有Jellyfish、Kmergenie、KmerFreq和GCE等。其中Kmergenie常用于预测de novo组装中最优组装k-mer大小,根据reads分割k-mers并绘制k-mer深度分布曲线。Jellyfish分析准确度高,常用于判断基因组纯度、杂合度等。

Kmergenie估计基因组大小

基因组大小可以通过k-mer分析法来估计[43]。假设基因组大小为G(也即一共有G个碱基),那么基因组可以产生的k-mers(genomic k-mers)数量为G-k 1,在G>>k的情况下,基因组k-mers数目就约等于基因组碱基数G。现在假定所有的genomic k-mer均为互不相同的(不考虑长于k的重复序列),其深度均为1。现在所有的测序reads均产生测序k-mers,由于测序深度较高,k-mers出现的频次也即k-mer深度较大,在去除错误率影响的前提下,可以认为其中完全不同的k-mer数目即是genomic k-mer数目,使完全不同的k-mer的数目最大的k值就是最佳k-mer大小。相反的从组装角度来讲,k越大则跨过基因组中重复序列的可能性越大,则完全不同的k-mer的数目越多,组装越容易,能够组装的序列越长,越接近实际基因组大小。

在实际测序分析中,由于覆盖度coverage并不是完美的(<100%),当k比较大时,reads产生的测序k-mers往往不能完全包含基因组k-mers,因此测序k-mers一般小于G;但是k越大,k-mers包含错误的概率也就越大,由于错误造成的低丰度k-mers越多,同时由于reads长度限制,reads产生的k-mers数目也越少,长度小于k的reads被去除,测序数据利用率降低导致导致覆盖度降低。当k比较小时,由于碱基数少,序列的种类就越少(例如4mer只有44=256种),再加上重复序列的影响,那么大的基因组其k-mer重复的可能性越大,基因组k-mers也即unique k-mers数目越小。可以想象,存在一个最佳的k使得测序unique k-mers的数目达到峰值,此峰值即为估计基因组大小,对应的k值为最佳k-mer大小。

Kmergenie(http://kmergenie.bx.psu.edu/)根据reads分割k-mers并绘制k-mer深度分布曲线,所有的k-mers包含基因组k-mers、测序错误或异源DNA污染k-mers等,其中基因组k-mer的copy次数为高斯分布,其峰值对应k-mer期望深度;错误k-mers的深度分布为帕累托分布(Pareto distribution),在去除错误碱基影响下,寻找使unique k-mers数目最大的k并估计基因组大小G,其使用方法如下所示:

代码语言:javascript复制
kmergenie -k 135 -l 15 -s 10 -t 10 -o out fq.list 1>std.log 2>err.log
参数说明:
-k  最大k-mer碱基数,默认值为121
-l  最小k-mer碱基数,默认值为15
-s  从最小到最大k-mer长度的取值间隔,默认值为10(也即15、25、35…)
-t  运行使用的核数
-o  输出文件的前缀,默认为histograms

其中fq.list内容是reads文件的路径,如下所示:

最大、最小k-mer长度以及取值间隔根据clean data中片段长度分布而定,对于Hiseq PE150测序数据一般设置25-135,如下所示:

代码语言:javascript复制
kmergenie -k 135 -l 25 -t 20 -o bacteria fq.list 1>std.log 2>err.log

运行完毕后打开生成的html文件即可看到结果报告,如下所示:

结果显示估计基因组大小为4.36M,最佳的k-mers大小为111kb。基因组大小随着不同size的k-mers变化曲线如下所示:

不同k-mer size的k-mer深度分布曲线也不相同,如下所示:

k较小时k-mer分布曲线出现多个峰,可能是基因组中有重复序列,也可能是有其他生物基因组污染。一定k-mer size的不同物种的基因组k-mer深度曲线具有特异性,在宏基因组分析中可以作为基因组指纹来分离不同基因组。

Jellyfish基因组k-mer分析

Jellyfish(http://www.cbcb.umd.edu/software/jellyfish/)输入文件为fasta/fastq格式的DNA序列文件,使用二进制文件输出k-mer数目信息(储存为Hash)。Jellyfish的功能有:kmer计数;融合二进制的Hash结果;统计Hash结果;通过Hash结果来画直方图;将Hash结果输出成文本格式;查询指定kmer的数目。

Jellyfish使用count的命令来执行计数功能,如下所示:

代码语言:javascript复制
jellyfish count -C -m 17 -s 100M -c 7 -t 30 -o mer_counts.jf R1.fastq R2.fastq
-m    设置k-mer长度(也即k值)。如果基因组大小为G,则k-mer长度选择为:k≈log(200G)/log(4)。
-s    设置Hash的大小。最好设置的值大于总的特异k-mer数,这样生成的文件只有一个。若该值不够大,则会生成多个hash文件。如果基因组大小为G,每个reads有一个错误,总共有n条reads,则该值可以设置为(G   n)/0.8。该值识别单位M和G。
-t    使用的CPU线程数(即核数)。
-o    后缀为jf的结果文件名
-c    k-mer的计数结果所占的比特数,默认为7,则支持的最大数字是2^7=128。如果基因组测序覆盖度为N,则要使设置的该c值满足2^c≥N。该值越大,消耗内存越大。
-C    对正义链和反义链都进行计数
-L    不输出数目低于此值的k-mer
-U    不输出数目高于此值的k-mer
-Q    设置碱基质量阈值,碱基质量低于该值则被转换为N

k-mer计数的结果以hash的二进制文件结果给出,需要统计出k-mer总数(Total),特异的k-mer数目(Distinct),只出现过一次的k-mer数(Unique),频数最高的k-mer数目(Max_count)等信息,以stats命令来运行,如下所示:

代码语言:javascript复制
jellyfish stats -o mer_counts_stats.txt mer_counts.jf

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

代码语言:javascript复制
jellyfish histo -l 2 -h 1000 -t 30 mer_counts.jf |sed 's/ /t/' > mer_counts_hist.txt
常用参数:
-l    最低的x轴值,也即最小的k-mer频数,默认为1。结果会将低于此值的所有的k-mer的数目作为(x‐1)的值。因此该值为2和1的结果是一致的。
-h    最高的x 轴值,默认为10000。同时结果会将高于此值的所有的k-mer的数目的和作为(x 1)的值。
-i    x轴取值间隔,每隔该数值取值,默认为1
-f    全部的直方图

对质控后的测序数据进行分析(根据前面选项中的方法计算参数值),如下所示:

代码语言:javascript复制
jellyfish count -C -m 15 -s 100M -c 8 -t 20 -o twk_kmer_counts.jf clean_1.fq clean_2.fq
jellyfish stats -o twk_kmer_stats.txt twk_kmer_counts.jf
jellyfish histo -l 2 -h 1000 -t 20 twk_kmer_counts.jf | sed 's/ /t/' > twk_kmer_hist.txt

使用生成的hist.txt文件可以绘制k-mer深度曲线,结果如下所示:

可以看到,Jellyfish的分析结果要远好于Kmergenie。

杂合度估计

杂合度是主要针对真核生物的一个概念,本意是指两亲本间的差异程度,也即来自两个不同亲本的DNA序列的差异程序。杂合度对基因组装的影响主要体现在如果二倍体或者多倍体杂合度较高,两个姊妹染色体上的等位基因会被识别为两个不同的基因,导致单倍体基因组偏大。此外微生物如果在培养时存在杂菌污染,也会造成类似影响。单倍体与多倍体可以使用试验来确定,在k-mer估计中可以根据k-mer深度曲线区分杂合峰重复峰,来估计基因组杂合度情况[45],如下所示:

①在最佳k-mer size深度曲线只有一个主峰,说明为单倍体或者纯合体;

②若在x=a处出现主峰,x=0.5a处有一个次峰,说明部分片段出现的期望值是大部分的1/2,当序列有杂合时,包含杂合位点的kmer因为分成了两部分,所以出现频率变为一半,次峰为杂合峰;

③若在x=a处出现主峰,x=2a处有一个次峰,说明一部分片段出现的期望值是大部分的2倍,这些片段为重复片段,次峰为重复峰;

④若出现两个主峰,峰高相差不大,两峰横坐标又是2倍关系,说明该个体高杂合或高重复。

0 人点赞