(全文约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)

  1. 参数

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

  3. -m 17

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

  4. -s 1000M

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

  5. -t 12

    线程12

  6. -C

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

  7. -o sample.jf

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

  8. c 7

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

  9. -out-counter-len=4

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

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

  11. -L

    不输入低于此值的k-mer

  12. -U

    不输入高于此值的k-mer

  13. 输入

  14. 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)。

  1. 参数

  2. -t 12

    线程12。

  3. -l 1

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

  4. -h 10000

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

  5. -i 1

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

  6. 后果
  7. k-mer频数散布直方表(sample.histo)蕴含空格分隔的两列数据。
  8. 第一列代表k值呈现的次数x(x=1,2,3...)，第二列是呈现了x次的kmer的品种的数量y。
  9. 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$V1Number <- kmer$V2png('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