jellyfish进行Kmer频数统计

参考 jellyfish:快速计算kmer分布 - 腾讯云开发者社区-腾讯云，代码部分改动

下载 jellyfish的压缩包，使用mv命令重命名，使用chmod命令给文件执行权限。

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

当我使用jellyfish命令运行，报错如下。./jellyfish可以运行。

 

1.计算kmer分布

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

jellyfish 只能接受fastq或者fasta文件作为输入。-m参数指定kmer的长度，-t指定并行的线程数，-s指定内存中hash的大小，该参数根据基因组的大小调整，test.fq是输入的序列文件。默认情况下会生成名为mer_counts的二进制文件，可以通过其他命令来查看该文件中的内容。

我的 test.fq:

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

因为我的test很小，且阿里云的内存较小，我调整了m=7，t=1，s=5M。

./jellyfish count -m 7 -t 1 -s 5M test.fq

执行后，发现生成了 mer_counts_0文件。

查看mer_counts_0文件

 

2.生成kmer 统计表

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

./jellyfish dump mer_counts_0 > kmer_count.fasta 

注意，文件名称要写对写全，如果你像参考博客一样生成的文件是mer_counts.jf，你需要写 ./jellyfish dump mer_counts.jf> kmer_count.fasta

kmer_count.fasta文件中每一条序列就是一个kmer，序列标识符是该kmer出现的次数。使用命令vim kmer_count.fasta查看该文件。

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

./jellyfish dump -c -t mer_counts_0 > kmer_count.tsv

 使用vim kmer_count.tsv查看该文件。第一列为kmer,第二列为该kmer频数。

 

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

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

./jellyfish query mer_counts_0 AAAGCAG

没有运行成功，报错如下，先插个眼。

 

4. 统计kmer基本信息

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

./jellyfish stats mer_counts_0

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

5. 统计kmer频数分布

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

./jellyfish histo mer_counts_0 > kmer_hist.tsv

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

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