歡迎關注"生信修煉手冊"！

jellyfish可以統(tǒng)計DNA序列中Kmer的分布，它運行速度快，內(nèi)存消耗低，支持并行，是最常用的kmer統(tǒng)計軟件之一。

官網(wǎng)如下：

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可執(zhí)行文件就可以了，代碼如下

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參數(shù)指定kmer的長度，-t指定并行的線程數(shù)，-s指定內(nèi)存中hash的大小，這個參數(shù)可以根據(jù)基因組的大小適當調整，比如人類基因組3G,這里就設置成3G;test.fq是輸入的序列文件。

需要注意的是，jellysifh 只能接受fastq或者fasta文件作為輸入，而且不能是壓縮文件。

默認情況下會生成名為mer_counts.jf的文件，該文件是一個二進制文件，可以通過其他命令來查看該文件中的內(nèi)容。

2. 生成kmer 統(tǒng)計表

dump子命令可以將上述步驟生成的二進制文件轉換成純文本的文件，用法如下

jellyfish dump mer_counts.jf > kmer_count.fasta

kmer_count.fasta文件中每一條序列就是一個kmer,序列標識符是該kmer出現(xiàn)的次數(shù)，示意如下

>1150
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>20
GCTACCATGATAGCCAAGGAAATCCCACAAA

也支持輸出成表格的形式，只需要添加-c和-t兩個參數(shù)，用法如下

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

輸出內(nèi)容如下

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 1150
GCTACCATGATAGCCAAGGAAATCCCACAAA 20
TATTTGCGGAGCTTGTTAGATAACAATCAAA 18

第一列為kmer,第二列為該kmer頻數(shù)。

3. 查詢某個kmer出現(xiàn)的次數(shù)

query 子命令可以快速查詢某個kmer的頻數(shù)，用法如下

jellyfish query mer_counts.jf GCTACCATGATAGCCAAGGAAATCCCACAAA

輸出結果：

GCTACCATGATAGCCAAGGAAATCCCACAAA 20

4. 統(tǒng)計kmer基本信息

stats子命令會給出kmer的基本統(tǒng)計信息，用法如下

jellyfish stats mer_counts.jf

會在命令行輸出如下統(tǒng)計結果

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

unique 代表只出現(xiàn)了1次的kmer個數(shù)；Distinct 代表kmer的個數(shù)，Total 代表所有kmer頻數(shù)總和，Max count 代表kmer的最大頻數(shù)。

5. 統(tǒng)計kmer頻數(shù)分布

histo子命令可以給出kmer的頻數(shù)分布，用法如下

jellyfish histo mer_counts.jf > kmer_hist.tsv

輸出內(nèi)容如下

1 130512636
2 24879103
3 12766220
4 8746042

第一列代表kmer的頻數(shù)，第二列代表出現(xiàn)該頻數(shù)的kmer的總數(shù)。利用這個數(shù)據(jù)，可以畫出kmer頻數(shù)分布曲線，對應的R語言代碼如下

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

最終的效果圖如下

掃描關注微信號，更多精彩內(nèi)容等著你！