禾本科作物基因组的kmer分布
做小麦的都知道麦类基因组庞大而且重复序列比例又很高,这给基因组分析造成了很大的困难,衡量基因组重复序列特征的一个方法就是统计kmer的分布情况。
这里简单的介绍下什么是kmer,所谓mer,大家可以理解成撕纸片,就是把基因组的碱基ATCG序列撕碎,怎么撕呢?假如k=4,就是撕成连续的4个碱基长度,k=20,就是就是撕成连续的20个碱基长度。假如基因组很小,只有ATCGCG,6个碱基组成,那么4mer就是ATCG,TCGG,CGCG。
kmer的分布可以估算基因组的大小,基因组的重复序列占比,杂合性等特征。可以快速的对基因组瞄一眼,看一看。所以利用重测序对基因组的kmer分布画图,一般也叫作调研图。调研吗,就是领导来走马观花的看一看,有时候准有时候不准,只是看看,有个印象,别太当真。
一段长序列,按照一定的长度分拆后,就会产生片段化的mer,这时候对这些短序列统计重复出现的次数,可以得到没有重复的mer和重复的mer,进而得到没有重复的mer占比。比如还是ATCGCG。2mers就是AT,TC,CG,GC,CG一共5个2mer,CG重复了一次,所以Uniqueness 2mers 比例就是3/5=0.6,4mers就是ATCG,TCGG,CGCG,没有重复出现的mers,因此Uniqueness 4mers 比例就是1。大家可以想象的到一段序列重复序列越多Uniqueness kmers比例越低。
今天给大家介绍一个软件叫Tallymer来做这个分析。这个软件其实是genometools软件包下面的小软件,详细介绍和安装(linux系统下载解压就可以使用)可以参考这里http://genometools.org/index.html。这里只做简单的介绍。
这个软件有三个主要命令mkindex, occratio, 和search,第一个是对一段序列建立索引和kmer库,occratio主要就是计算Uniqueness kmers 比值。结果有三列,第一列是mers大小,第二列是出现的频次数,第三列是出现的频率。数据有了,下面画图就是小意思了。
1$gt suffixerator -dna -pl -tis -suf -lcp -v -parts 4 -db read1.fna -indexname reads
2
3 $ gt tallymer occratio -scan -output unique relative -minmersize 10 -maxmersize 20 -esa reads
4# distribution of unique mers
510 223755 0.623
611 373775 0.802
712 444083 0.877
813 465859 0.904
914 468735 0.913
1015 465646 0.917
1116 460791 0.919
1217 455449 0.920
1318 450049 0.921
1419 444720 0.921
1520 439532 0.922
16# space peak in megabytes: 0.06
17# mmap space peak in megabytes: 3.93
这里贴上几个小麦相关的kmer美图,大家以后碰到这种图就知道怎么回事了。
小麦染色体特征Circos图。
这个图大家一看就知道哪里来的。这里第三圈,C圈,就是小麦Uniqueness 20-mers的频率,着丝粒的地方明显有一道黄线闪过,显示重复序列更高。
比较意外的是我在用这个软件分析麦类作物的时候,发现和节节麦DD基因组发表时候kmer分布结果差别比较大。
图2.麦类作物uniuqe kmers分布图。
图3. DD基因组发表时不同物种uniuqe kmer分布图。
图2是我下载已发表的基因组序列算的,图3是DD基因组发表时候的结果。可以看到我算出来的结果比图3的Aet(DD)大的多。这是最先怀疑的是自己算错了,但是我又去比较了软件发表时,软件作者的结果,可以看到图4和我的结果比较类似,水稻虽然TE含量比较少,早早的就达到了90%以上,并不像图3,一直在80%以下。
图4. Tallymer软件发表时计算的几个物种uniuqe kmer分布图。
有小伙伴知道原因吗?欢迎来讨论。
参考文献:
Kurtz S, Narechania A, Stein JC, Ware D. A new method to compute K-mer frequencies and its application to annotate large repetitive plant genomes.BMC Genomics. 2008 doi: 10.1186/1471-2164-9-517.
Luo MC, Gu YQ, Puiu D,Wang H, et al., Genome sequence of the progenitor of the wheat D genome Aegilops tauschii.Nature. 2017 Nov 23;551(7681):498-502. doi: 10.1038/nature24486.