(全文约1520字)
【推荐】用Smudgeplot评估物种倍性后,用组合jellyfish+GenomeScope1.0做二倍体物种的基因组调查,用组合KMC+GenomeScope2.0做多倍体物种的基因组调查。
k-mer进行基因组调查分为k-mer频数统计和基因组特征评估两步。
jellyfish是Center for Bioinformatics and Computational Biology在2011年研发的一款对DNA的k-mers计数的软件,用Hash表储存数据,能多线程运行。
conda install -c bioconda jellyfish
#安装的是v2.2.10一般先用jellyfish count
进行k-mer计数,然后用jellyfish histo
对结果进行统计,获得k-mer的频数分布直方表sample.histo。
jellyfish count -m 17 -s 10G -t 12 -C -o sample.jf <(zcat sample_1.fq.gz) <(zcat sample_2.fq.gz)
使用的PE reads,不支持压缩格式*.fq.gz输入文件,如果不解压缩,也可以用<(zcat sample_1.fq.gz) <(zcat sample_2.fq.gz)
代替sample_1.fq sample_2.fq
; 或者使用这种形式zcat *fq.gz | jellyfish count /dev/fd/0
,其中/dev/fd/0
是进程输入标志,代表管道前结果传递。
存储用的hash表大小为1000M,这个参数识别单位M(Mbp)和G(Gbp)。若该值不够大,则会生成多个hash文件,以数字区分文件名。最好设置的值大于总的独特的(distinct)k-mer数,这样生成的文件只有一个。如果基因组大小为G,每个reads有一个错误,总共有n条reads,则该值可以设置为[(G + n)/0.8]。
结果文件名为sample.jf,会生成k-mer计数文件sample.jf,是hash的二进制文件。
k-mer的计数结果所占的最大比特数,默认支持的最大数字是2^7=128。该值最大,消耗内存越大。
输出的二进制hash文件中的计数结果所占的字节数,一个字节是8比特,则默认支持的最大数字是2^32=4.3G。
不输出低于此值的k-mer
hash格式储存的k-mer频数文件
jellyfish histo -t 12 sample.jf > sample.histo
统计k-mer计数(sample.jf)得到k-mer频数分布直方表(sample.histo)。
x的最小值,默认是1。结果会将小于此值的所有的k-mer的数目作为(x‐1)的值总结到一行。
x轴取值间隔,每隔该数值取值,默认为1。
如果jellyfish count模块输出结果的二进制hash文件有多个,需要将多个hash文件合并,合并到merge.jf。
jellyfish merge sample_hash1.jf sample_hash2.jf sample_hash3.jf -o merge.jf
jellyfish stats sample.jf -o counts_stats.txt
可以用stats模块来统计出k-mer总数(Total),特异的k-mer数目(Distinct),只出现过一次的k-mer数量(Unique),频数最高的k-mer数量(Max_count)等信息。
获得k-mer频数分布表sample.histo后,推荐用GenomeScope1.0或者GenomeScope2.0或者GenomeScope的R脚本来做基因组特征评估和画图。也可直接用R绘制sample.histo的频率分布直方图/频率分布曲线。
GenomeScope2.0 网页版也是类似的步骤。
R绘制k-mer频数分布曲线初步查看基因组特征。
获得kmer_plot.png为频数分布曲线,可根据曲线峰值对基因组大小进行计算和预估。
#R 脚本示例
kmer <- read.table('sample.histo')
kmer <- subset(kmer, V1 >=5 & V1 <=500) #对频数范围5-500的数据进行绘制
Frequency <- kmer$V1
Number <- kmer$V2
png('kmer_plot.png')
plot(Frequency, Number, type = 'l', col = 'blue')
dev.off()