【如何使用Jellyfish进行Kmer频数统计分析？】2022

2022-05-11 本文已影响0人陈有朴

软件概述

将基因组等序列文件，按长度为K的字符串进行切割、归类以及频数统计。

软件安装

软件版本：2.3.0

wget https://github.com/gmarcais/Jellyfish/releases/download/v2.3.0/jellyfish-2.3.0.tar.gz
# 1151287 (1.1M)
# 安装包就1.1M，可以说是压缩效率非常高了

./configure prefix=/opt/biosoft/Jellyfish
make
make install 
echo 'PATH=/opt/biosoft/Jellyfish/bin:$PATH' >> ~/.bashrc
source ~/.bashrc

软件使用

可以直接查看帮助文档：man jellyfish

接下来以大肠杆菌为例，对Jellyfish的使用进一步熟悉。

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
gunzip -c GCF_000005845.2_ASM584v2_genomic.fna.gz > Ecoli.fasta

需要注意的几个问题：

jellyfish的输入文件格式不能为压缩状态
```
# 针对压缩的情况
<(zcat input1.fq.gz)
```
-C，canonical的理解

（1）Kmer计数

jellyfish count -m 15 -o Ecoli.jf -c 3 -s 1G -t 16 Ecoli.fasta
# 参数说明
# -m，Kmer长度为15
# -o，设置输出文件名
# -t，设置计数过程所使用的线程数
# Ecoli.fasta，为输入文件

对于-c参数和-s参数的描述：

-c：The counter field in the hash uses only <font color='yellow'>3 bits</font> and the hash has at least 10 million entries.
-s：哈希表格的大小。当设置的哈希表格足够大，能够容纳所有的Kmer计数结果时，只有一个结果文件，但当哈希表格的entry数设置小了的话，就会有多个结果。

Jellyfish 2.0版本已经会自动合并结果了。

（2）输出文本Kmer计数结果

默认情况下，Jellyfish的结果文件.jf使用哈希来对结果进行保存，那么如何输出文本格式的计数结果呢？

1、生成fasta格式的计数结果

jellyfish dump Ecoli.jf > Ecoli_output_count.fasta

结果展示如下：

>1
AACGATACAGCCTTT
>1
AACGATACAGCGAAA
>1
AACGATACAGTAAAT
>1
AACGATACAGTATCT
>1
AACGATACATACAAA
>1
AACGATACATGACAG
>1
AACGATACATTTCTT

2、按列存储信息 & 以TAB作为分隔输出Kmer计数结果

# 按列分别存储信息的计数结果
jellyfish dump -c -t Ecoli.jf > Ecoli_output_col_count.tsv

结果展示如下：

AACCGTAGGCCGGAT 12
AACGCCTTATCCGGC 50
AAGCGTAGCGCATCA 26
AAGCGTCGCATCAGG 11
AAGGCGTTCACGCCG 46
AAGGCGTTTACGCCG 34
AATGCCTGATGCGAC 18
AATTCAATATATTGC 10
ACCGTAGGCCGGATA 12
ACGCCGCATCCGGCA 71
ACGCCTTATCCGGCC 46
ACGCGTCTTATCAGG 15
ACGCTGGCGCGTCTT 22
ACGCTTATCAGGCCT 26
AGCGCATCAGGCAAT 16

（3）输出Kmer频数分布

jellyfish histo Ecoli.jf > Ecoli_output.histo

（4）查询特定Kmer频数

查询的速度是相当快的，比count过程快了太多。

# e.g. 
jellyfish query Ecoli.jf ACGCTTATCAGGCCT

（5）统计Kmer计数结果

jellyfish stats Ecoli.jf

输出结果如下：

Unique:    4443336         # 频数为1的Kmer总数
Distinct:  4517632         # Kmer种类
Total:     4641394         # Kmer频数和
Max_count: 71              # Kmer最大频数

参考资料

[1] https://github.com/gmarcais/Jellyfish/releases/tag/v2.3.0

额外

-C，canonical是什么含义？

这部分需要针对数据类型来分析 —— 只有对自己的数据有一个好的理解，做出的分析才是可信的。

对于fasta类型数据来说，这只是一条单链的DNA序列，本身就不存在的反向重复，因此在对其进行分析的时候，没有必要使用-C参数。

对于fastq类型数据来说，由于测序时将基因组DNA（当然也可以是RNA）打断之后，对2条链都进行了测序。

但是在进行基因组大小估计、选择亚基因组特异性Kmer等分析时，都是默认一条DNA序列的情况，因此在针对fastq类型的数据进行分析时，我们就需要加上-C参数。

基因组大小估计公式： $Genome Size = \frac{Knum}{Kdepth}$

如上解释的参考资料：https://www.biostars.org/p/153170/