jellyfish:快速计算kmer分布

欢迎关注"生信修炼手册"

jellyfish可以统计DNA序列中Kmer的分布,它运行速度快,内存消耗低,支持并行,是最常用的kmer统计软件之一。

官网如下:

http://www.genome.umd.edu/jellyfish.html

jellyfish有两个大版本,1.0版本和2.0版本,1.0版本最多支持kmer长度为31, 2.0版本对kmer长度没有任何的限制。在下载时,注意下载2.0版本。

github链接如下:

https://github.com/gmarcais/Jellyfish

直接下载编译好的linux可执行文件就可以了,代码如下

wget https://github.com/gmarcais/Jellyfish/releases/download/v1.1.12/jellyfish-linux
mv jellyfish-linux jellyfish
chmod +x jellyfish

jellyfish有很多个子命令,下面介绍常用的子命令

1. 计算kmer分布

count子命令用来计算kmer分布,用法如下

jellyfish count -m 31 -t 10 -s 1G test.fq

-m参数指定kmer的长度,-t指定并行的线程数,-s指定内存中hash的大小,这个参数可以根据基因组的大小适当调整,比如人类基因组3G,这里就设置成3G;test.fq是输入的序列文件。

需要注意的是,jellysifh 只能接受fastq或者fasta文件作为输入,而且不能是压缩文件。

默认情况下会生成名为mer_counts.jf的文件,该文件是一个二进制文件,可以通过其他命令来查看该文件中的内容。

2. 生成kmer 统计表

dump子命令可以将上述步骤生成的二进制文件转换成纯文本的文件,用法如下

jellyfish dump mer_counts.jf > kmer_count.fasta

kmer_count.fasta文件中每一条序列就是一个kmer,序列标识符是该kmer出现的次数,示意如下

>1150
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>20
GCTACCATGATAGCCAAGGAAATCCCACAAA

也支持输出成表格的形式,只需要添加-c-t两个参数,用法如下

jellyfish dump  -c -t mer_counts.jf > kmer_count.tsv

输出内容如下

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 1150
GCTACCATGATAGCCAAGGAAATCCCACAAA 20
TATTTGCGGAGCTTGTTAGATAACAATCAAA 18

第一列为kmer,第二列为该kmer频数。

3. 查询某个kmer出现的次数

query 子命令可以快速查询某个kmer的频数,用法如下

jellyfish query mer_counts.jf GCTACCATGATAGCCAAGGAAATCCCACAAA

输出结果:

GCTACCATGATAGCCAAGGAAATCCCACAAA 20
4. 统计kmer基本信息

stats子命令会给出kmer的基本统计信息,用法如下

jellyfish stats mer_counts.jf

会在命令行输出如下统计结果

Unique:    130512636
Distinct:  260032104
Total:     3502352808
Max_count: 43395

unique 代表只出现了1次的kmer个数;Distinct 代表kmer的个数,Total 代表所有kmer频数总和,Max count 代表kmer的最大频数。

5. 统计kmer频数分布

histo子命令可以给出kmer的频数分布,用法如下

jellyfish histo mer_counts.jf > kmer_hist.tsv

输出内容如下

1 130512636
2 24879103
3 12766220
4 8746042

第一列代表kmer的频数,第二列代表出现该频数的kmer的总数。利用这个数据,可以画出kmer频数分布曲线,对应的R语言代码如下

x <- read.table(input, header = F, sep = " ", stringsAsFactors = F)
pdf("kmer_hist.pdf")
plot(x, type = "l")
dev.off()

最终的效果图如下

扫描关注微信号,更多精彩内容等着你!

  • 6
    点赞
  • 21
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值