lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

基因组K-mer分析软件JELLYFISH的安装及使用说明

已有 6777 次阅读 2019-2-10 18:28 |个人分类:Software/Pipeline/Script|系统分类:科研笔记| _生信工具, _K-mer分析

基因组K-mer分析软件JELLYFISH的安装及使用说明

 

一般来讲,对于较为复杂的基因组,我们通常会在基因组正式组装之前进行k-mer分析,以评估基因组杂合度、重复序列比例,并初步预估基因组大小,以及测序深度等,为后续基因组的正式组装提供参考信息。k-mer分析简介可参见:http://blog.sciencenet.cn/blog-3406804-1162384.html

此处给大家介绍一款用于k-mer分析工具,JELLYFISH。该软件运用Hash表来存储数据,同时能多线程运行,具有速度快,内存消耗小的优点。

JELLYFISH官方网站:http://www.cbcb.umd.edu/software/jellyfish/

 

本文中的测试数据以及结果示例的百度盘链接如下,均无提取码。

测试数据:https://pan.baidu.com/s/1NFvRDs3kZsh3VK7ukv1svA

部分示例结果:https://pan.baidu.com/s/10mO7isvGgxJ7uFQkTVOdSQ

 


JELLYFISH下载安装


JELLYFISH有两个大版本,1.0版和2.0版。1.0版最多支持k-mer长度为312.0版对k-mer长度没有任何的限制,因此安装时最好使用2.0版的。以下两种方式二选一即可,推荐conda安装,方便快捷。JELLYFISH配置好后即可使用。

 

一,在Linux系统下,可直接使用conda安装。

#conda 安装,默认 2.0 版本的
conda install jellyfish


二,在Linux系统下源码编译的话,若无root权限,记得建立一个配置安装的路径。

注:源码编译要求g++ > =4.4

#比方说我这里将配置安装的路径设定在“/home/my/software/jellyfish”
#安装 2.2.10 版本,下载地址 https://github.com/gmarcais/Jellyfish/releases
tar -zxvf jellyfish-2.2.10.tar.gz && cd jellyfish-2.2.10 && chmod +x configure
mkdir -p /home/my/software/jellyfish
./configure --prefix=/home/my/software/jellyfish
make
make install
 
#编译完成后,在运行 jellyfish 前需执行以下命令,可直接在 ~/.bashrc 中添加
export PATH=$PATH:/home/my/software/jellyfish/bin
export LD_LIBRARY_PATH=/home/my/software/jellyfish/lib
export MANPATH=/home/my/software/jellyfish/share/man
export PKG_CONFIG_PATH=/home/my/software/jellyfish/lib/pkgconfig

 


使用JELLYFISH统计K-mer


JELLYFISH具有多种功能。其中常用功能:

countk-mer统计计数;

stats,统计k-mer计数所得结果总概;

histo,统计k-mer计数生成k-mer频数统计表;

dump,将k-mer统计计数所得二进制的结果转化为纯文本内容;

merge,合并多个由k-mer统计计数所得结果;

query,在k-mer计数结果中查询指定的k-mer出现的次数。

1.png

 

jellyfish count <options>k-mer统计计数


该命令可根据预设k-mer长度,将测序reads序列截为一段段的k-mer序列,结果以hash表的结构存储。除了测序reads序列,对较长的基因组序列同样可以统计。输入文件可以为fastq文件,也可以为fasta文件。得到k-mer统计计数结果是该程序最重要的功能,也是执行后续统计分析不可缺少的第一步。

可首先使用“jellyfish count --help”查看参数。

jellyfish count --help

2.png

其中,常用参数选项说明:

-m,设定选取的k-mer长度;

-s,设定存储k-mer序列所用hash表的大小,单位为MG;若该值足够大则最终得到一个hash文件,若该值不够大则会生成多个hash文件;

-t,程序运行所需的线程数,默认单线程;

-o,输出结果文件的前缀名;

-C,对基因组DNA正链和负链都进行统计;

-L,不输出数目低于此值的k-mer

-U,不输出数目高于此值的k-mer

 

此处使用网盘中提供的细菌二代测序reads数据(Bacillus_subtilis.filt_R1.fastq.gz/Bacillus_subtilis.filt_R2.fastq.gz,在clean reads的基础上进行碱基纠错及reads去重后所得),简要展示JELLYFISH的使用。

#JELLYFISH 不识别 .gz 格式,首先解压缩得到纯文本 fastq
gunzip -c Bacillus_subtilis.filt_R1.fastq.gz > Bacillus_subtilis.filt_R1.fastq
gunzip -c Bacillus_subtilis.filt_R2.fastq.gz > Bacillus_subtilis.filt_R2.fastq
#之后执行 JELLYFISH
jellyfish count -m 17 -s 100M -t 4 -C Bacillus_subtilis.filt_R1.fastq Bacillus_subtilis.filt_R2.fastq
#或者
jellyfish count -m 17 -s 100M -t 4 -C fastq_list.txt

可以直接在jellyfish参数命令行的末尾添加fastq文件的名称(若文件已在当前路径中)或绝对路径;也可以将fastq文件的绝对路径写入至一个文本文档中,在jellyfish参数命令行的末尾添加该文本文档的绝对路径。假设示例中的fastq文件位于“/home/my/data”,则文本文档“bac_fastq_list.txt”中的内容可以这么写:

3.png

jellyfish参数命令中,选取k-mer长度为17,存储使用hash表的大小为100M,程序运行线程数为4,同时统计DNA正链和负链。

 

程序运行完毕后,当前工作路径下得到结果文件“mer_counts.jf”,以二进制存储结果,包含了从测序reads数据中截取所得的各k-mer序列信息。之后可使用该文件进行后续的处理了。

4.png

 

jellyfish stats <options>,统计k-mer计数所得结果总概


在得到k-mer统计计数结果后,可使用该命令查看k-mer总数、出现最多次数的k-mer数等简要信息。

可首先使用“jellyfish stats --help”查看参数。

jellyfish stats --help

5.png

其中,常用参数选项说明:

-L,不统计数目低于此值的k-mer

-U,不统计数目高于此值的k-mer

-o,输出结果文件的前缀名。

 

我们对一开始k-mer统计计数所得的结果文件“mer_counts.jf”进行简要统计,展示其概要内容。

jellyfish stats mer_counts.jf

参数中无特殊要求,直接出结果。

 

统计结果中:

Unique,只出现过一次的k-mer的数目总和,此处为890140

Distinct,特异性的k-mer数目,包含Unique在内,此处为5215828

Total,所有k-mer的数目总和,此处为260999943

Max_count,出现次数最多的k-mer的数目,此处为94932

6.png

 

jellyfish histo <options>,统计k-mer计数生成k-mer频数统计表


在得到k-mer统计计数结果后,可使用该命令统计各k-mer出现的频数,得到我们所期望的k-mer频数分布列表。

可首先使用“jellyfish histo --help”查看参数。

jellyfish histo --help

7.png

其中,常用参数选项说明:

-l,设定选取的k-mer频数最低值,同时结果会将低于此值的所有的k-mer的数目作为(l-1)的值,默认值1

-h,设定选取的k-mer频数最高值,同时结果会将高于此值的所有的k-mer的数目的和作为(h+1)的值,默认值10000

-i,设定k-mer频数的统计间隔;

-t,程序运行所需的线程数,默认单线程;

-f,将中间的0值也列出。

 

继续基于一开始k-mer统计计数所得的结果文件“mer_counts.jf”,计算得到k-mer频数统计列表。

jellyfish histo mer_counts.jf > kmer_hist.txt

参数中无特殊要求,均使用默认参数。

 

8.png

文件“kmer_hist.txt”即为得到的k-mer频数统计表,共两列内容,以空格分隔。其中,第一列为k-mer深度,即各k-mer的出现频数;第二列为出现该频数的k-mer片段总数。

9.png

 

得到k-mer频数统计表后,我们即可使用该文件绘制k-mer频数分布图,以初步查看测序质量、物种基因组特征等;并可以借用其他工具进行更进一步的分析,如GCE等工具,包括估算基因组大小、杂合度、重复序列比例以及测序深度等指标。

例如,以下为使用该k-mer频数统计表绘制的k-mer曲线图的一个简单R脚本,绘图时只绘制k-mer出现频数在5-500区间范围内的统计结果。

#R 脚本示例
kmer <- read.table('kmer_hist.txt')
kmer <- subset(kmer, V1 >=5 & V1 <= 500)
Depth <- kmer$V1
Number <- kmer$V2

png('kmer_plot.png')
plot(Depth, Number, type = 'l')
dev.off()

k-mer曲线作图结果如下。可轻易看到该物种的k-mer曲线显著单峰,无杂合峰且未见明显重复峰,即该物种为单倍体,且基因组中几乎没有重复序列(这是可以预料的,因为测试数据来源于某细菌)。

10.png

 

jellyfish dump <options>,将k-mer统计计数所得二进制的结果转化为纯文本内容


在得到k-mer统计计数结果后,由于初始计数结果以二进制的hash表存储,不方便查看,因此可使用该命令将二进制文件转化为纯文本文件,便于后续自定义的统计处理。

可首先使用“jellyfish dump --help”查看参数。

jellyfish dump --help

11.png

其中,常用参数选项说明:

-c,生成列表样式的结果,包含2列内容,第一列为k-mer序列,第二列为对应的数目;若缺省,默认生成fasta格式的结果,即以k-mer序列为fasta的序列,对应的数目为fastatitle

-t,当生成列表样式时,以tab作为两列的分隔,默认以空格作为分隔;

-L,不输出小于该值的k-mer

-U,不输出高于该值的k-mer

-o,输出文件名称。

 

我们将一开始k-mer统计计数所得的结果文件“mer_counts.jf”(二进制存储),转换为纯文本内容的列表样式,并以tab分隔两列内容。

jellyfish dump -c -t -o mer_counts.txt mer_counts.jf

 

输出结果文件“mer_counts.txt”,查看其内容, 第一列为得到的各k-mer序列,第二列为各k-mer序列对应的数目。得到纯文本文件后,我们即可在后续自定义编辑统计使用等。

12.png

13.png

 

jellyfish merge <options>,合并多个由k-mer统计计数所得结果


若想将多个由k-mer统计计数所得的结果文件(即“mer_counts.jf”)合并在一起的话,可使用该命令完成。

该命令主要针对于由一条jellyfish count命令得到多个.jf文件的情形。前面第一步在对k-mer统计计数时,我们通过-s参数,设定了存储k-mer序列所用hash表的大小。并且说到,若该值足够大则最终得到一个hash文件,若该值不够大则会生成多个hash文件。hash文件转化为纯文本后,其内容可见上述,记录了每种k-mer及其出现频数。若因设定的hash表存储空间大小偏小,导致得到多个hash文件后,此时需要对结果进行合并,可使用jellyfish merge命令。

使用“jellyfish merge --help”查看参数。

jellyfish merge --help

14.png

其中,常用参数选项说明:

-o,输出文件名称;

-L,不输出小于该值的k-mer

-U,不输出高于该值的k-mer

 

使用的话,一个简单的示例命令如下。其中,“mer_counts.jf”为合并后的结果,“ mer_counts1.jf”“ mer_counts2.jf”“ mer_counts3.jf”为合并前的子文件。

jellyfish merge -o mer_counts.jf mer_counts1.jf mer_counts2.jf mer_counts3.jf

由于网盘中提供的示例测序数据总量较小,没有超出我们对hash的设定值仅得到了一个“mer_counts.jf”文件,因此不再进行实际的示例展示。

 

jellyfish query <options>,在k-mer计数结果中查询指定的k-mer出现的次数


在得到k-mer统计计数结果后,若想查看某个特定的k-mer出现的频数,则可以使用该命令。

 

例如,查询一开始k-mer统计计数所得的结果文件“mer_counts.jf”中,某特定k-mer“ATAACGATTCATATACA”出现的次数,命令如下。

jellyfish query mer_counts.jf ATAACGATTCATATACA

得到查询结果58

15.png

 


参考文献


Marçais G, Kingsford C. A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics, 2011, 27(6):764.

 




http://blog.sciencenet.cn/blog-3406804-1161522.html

上一篇:二代测序数据质量评估软件FastQC的安装及使用说明
下一篇:K-mer分析及基因组特征评估软件GCE的安装及使用说明

1 张成岗

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-12 19:10

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部