k-mer基因组调查分为k-mer频数统计和基因组特征评估两步。

• jellyfish第一步可以实现k-mer频数统计。
• jellyfish的结果sample.histo可以用在GenomeScope实现基因组特征评估的第二步。

jellyfish是Center for Bioinformatics and Computational Biology2011年开发的一款DNA的k-mers使用计数软件Hash存储数据，可以多线程运行。

1. conda安装

• conda install -c bioconda jellyfish #安装的是v2.2.10

2. github安装

• 在github：jellyfish通过源码安装。

一般先用jellyfish count进行k-mer计数，然后用jellyfish histo统计结果，获得结果k-mer频数分布直方表sample.histo。

1. 命令

jellyfish count -m 17 -s 10G -t 12 -C -o sample.jf <(zcat sample_1.fq.gz) <(zcat sample_2.fq.gz)

2. 参数

• sample_1.clean.fq sample_2.clean.fq

使用的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代表管道前结果传输的过程输入标志。

• -m 17
  
  k-mer长度设置为17bp。如果基因组大小为G(单位是bp)，k-mer长度推荐设置为log(200*G)/log(4)。500Mbp对应的基因组约为17、1Gbp的19，10Gbp的21。
• -s 1000M

存储用的hash表大小为1000M，该参数识别单位M(Mbp)和G(Gbp)。如果值不够大，就会产生多个hash文件名用数字区分。最佳设置值大于总的独特值(distinct)k-mer这样生成的文件只有一个。若基因组大小为G，每个reads有一个错误，总共有n条reads，该值可设置为[(G n)/0.8]。

• -t 12
  
  线程12
• -C
  
  对DNA正负链统计表示考虑DNA正义与反义链相遇kmer时间，计入正义kmer如果是双端测序。reads，需要这个参数。
• -o sample.jf

结果文件名称sample.jf，会生成k-mer计数文件sample.jf，是hash二进制文件。

• c 7

1. 命令

jellyfish histo -t 12 sample.jf > sample.histo

统计k-mer计数(sample.jf)得到k-mer直方表的频数分布(sample.histo)。

2. 参数

• -t 12
  
  线程12。
• -l 1

x默认情况下，最小值为1。结果将小于此值的所有值k-mer的数目作为(x‐1)将值总结到一行。

• -h 10000
  
  x默认情况下，最大值为1万。结果将大于所有这个值k-mer的数目作为(x 1)将值总结到一行。
• -i 1

x轴取值间隔，默认为1。

3. 结果

• k-mer直方表的频数分布(sample.histo)两列数据包含空间分隔。
• 第一列代表k值的次数x(x=1,2,3…)，第二列出现x次kmer类型数量y。
• sample.histo的两列即是kmerX和分布频率直方图y轴的值。

如果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或者GenomeScopeR脚本用于基因组特征评估和绘图。也可以直接用R绘制sample.histo频率分布直方图/频率分布曲线。

1. 在GenomeScope1.0 网页版上传前一步获得的k-mer频数分布表sample.histo文件。
2. 设置参数k-mer length为第一步选择的k-mer这是17；参数Read length序列读长一般为150；最后一个参数Max kmer coverage建议修改为更大的10000，统计更多k-mers。
3. 结果显示了基因组基因组大小、杂合度、重复率等信息。

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()

1. jellyfish paper：https://academic.oup.com/bioinformatics/article/27/6/764/234905?login=true
2. jellyfish github：https://github.com/gmarcais/Jellyfish
3. jellyfish参数推荐：https://www.bilibili.com/read/cv16360242
4. chenlianfu blog: jellyfish参数推荐：http://www.chenlianfu.com/?p=806