关于经验:用kmer分析进行基因组调查二用jellyfish进行kmer频数统计

(全文约 1520 字)

【举荐】用 Smudgeplot 评估物种倍性后，用组合 jellyfish+GenomeScope1.0 做二倍体物种的基因组考察，用组合 KMC+GenomeScope2.0 做多倍体物种的基因组考察。

1. k-mer 进行基因组考察的软件

k-mer 进行基因组考察分为 k-mer 频数统计 和基因组特色评估 两步。

• jellyfish 能够实现第一步 k -mer 频数统计。
• jellyfish 的后果 sample.histo 能够用在 GenomeScope 上，实现第二步基因组特色评估。

2. jellyfish 简介

jellyfish 是 Center for Bioinformatics and Computational Biology 在 2011 年研发的一款对 DNA 的 k -mers 计数的软件，用 Hash 表贮存数据，能多线程运行。

3. jellyfish 装置

1. conda 装置
2. conda install -c bioconda jellyfish #装置的是 v2.2.10
3. github 装置
4. 在 github：jellyfish 上通过源码装置。

4. jellyfish 运行

个别先用 jellyfish count 进行 k -mer 计数，而后用 jellyfish histo 对后果进行统计，取得 k -mer 的频数散布直方表 sample.histo。

4.1. count —— k-mer 计数

1. 命令

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

2. 参数

3. 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 是过程输出标记，代表管道前后果传递。

4. -m 17

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

5. -s 1000M

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

6. -t 12

   线程 12

7. -C

   对 DNA 正负链都进行统计，示意思考 DNA 正义与反义链，遇到反义 kmer 时，计入正义 kmer 频数中。如果是双端测序 reads，须要这个参数。

8. -o sample.jf

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

9. c 7

   k-mer 的计数后果所占的最大比特数，默认反对的最大数字是 2^7=128。该值最大，耗费内存越大。

10. -out-counter-len=4

    输入的二进制 hash 文件中的计数后果所占的字节数，一个字节是 8 比特，则默认反对的最大数字是 2^32=4.3G。

11. 不举荐用 -Q，会将低质量的碱基替换成 N。

12. -L

    不输入低于此值的 k -mer

13. -U

    不输入高于此值的 k -mer

14. 输入

15. sample.jf

    hash 格局贮存的 k -mer 频数文件

4.2. histo —— 统计 k -mer 频率

1. 命令

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

统计 k -mer 计数 (sample.jf) 失去 k -mer 频数散布直方表(sample.histo)。

2. 参数

3. -t 12

   线程 12。

4. -l 1

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

5. -h 10000

   x 的最大值，默认是 10000。后果会将大于此值的所有的 k -mer 的数目作为 (x+1) 的值总结到一行。

6. -i 1

   x 轴取值距离，每隔该数值取值，默认为 1。

7. 后果
8. k-mer 频数散布直方表 (sample.histo) 蕴含空格分隔的两列数据。
9. 第一列代表 k 值呈现的次数 x(x=1,2,3…)，第二列是呈现了 x 次的 kmer 的品种的数量 y。
10. sample.histo 的两列即是 kmer 散布频率直方图的 x 和 y 轴的值。

4.3. merge 合并【按需抉择】

如果 jellyfish count 模块输入后果的二进制 hash 文件有多个，须要将多个 hash 文件合并，合并到 merge.jf。

jellyfish merge sample_hash1.jf sample_hash2.jf sample_hash3.jf -o merge.jf

4.4. stats 统计【可选】

jellyfish stats sample.jf -o counts_stats.txt

能够用 stats 模块来统计出 k -mer 总数（Total），特异的 k -mer 数目（Distinct），只呈现过一次的 k -mer 数量（Unique），频数最高的 k -mer 数量（Max_count）等信息。

5. 基因组特色评估

取得 k -mer 频数散布表 sample.histo 后，举荐用 GenomeScope1.0 或者 GenomeScope2.0 或者 GenomeScope 的 R 脚本来做基因组特色评估和画图。也可间接用 R 绘制 sample.histo 的频率散布直方图 / 频率散布曲线。

5.1. GenomeScope 网页版

5.1.1. GenomeScope1.0 网页版 —— 实用于二倍体物种

1. 在 GenomeScope1.0 网页版上传前一步取得的 k -mer 频数散布表 sample.histo 文件。
2. 设置参数 k -mer length 为第一步抉择的 k -mer 长度值，这里是 17；参数 Read length 为序列读长，个别为 150；最初一个参数 Max kmer coverage 倡议批改成更大的 10000，以统计更多的 k -mers。
3. 结果显示预估的基因组大小，杂合度，反复率等信息。

5.1.2. GenomeScope2.0 网页版 —— 实用于多倍体物种

GenomeScope2.0 网页版也是相似的步骤。

5.2. R 绘制

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

6. references

1. jellyfish paper：https://academic.oup.com/bioi…
2. jellyfish github：https://github.com/gmarcais/J…
3. jellyfish 参数举荐：https://www.bilibili.com/read…
4. chenlianfu blog: jellyfish 参数举荐：http://www.chenlianfu.com/?p=806