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

博文

K-mer分析及基因组特征评估软件GCE的安装及使用说明

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

K-mer分析及基因组特征评估软件GCE的安装及使用说明

 

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

此处给大家介绍一款可用于k-mer分析及基因组特征评估的工具,GCE。该软件既可用于k-mer频数统计,也可基于k-mer频数统计结果对基因组大小、杂合度、重复序列、测序深度等进行评估。但是该软件似乎自2011年的1.0版本后,未见更新。

GCE下载链接:ftp://ftp.genomics.org.cn/pub/gce/

 

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

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

部分示例结果:https://pan.baidu.com/s/123YaKBhle9eONX_etxuVhA

 


GCE下载安装


Linux系统下,直接在下载链接ftp中下载安装。

GCE主要包含两个主程序:kmer_freq_hash,用于k-mer频数统计;gce,用于基因组特征评估。

mkdir -p /home/my/software/ && cd /home/my/software/
wget ftp://ftp.genomics.org.cn/pub/gce/gce-1.0.0.tar.gz
tar zxfv gce-1.0.0.tar.gz && cd gce-1.0.0
make
 
#更改指令方便今后的书写,可直接在 ~/.bashrc 中修改
alias kmer_freq_hash='/home/my/software/gce-1.0.0/kmerfreq/kmer_freq_hash/kmer_freq_hash'
alias gce='/home/my/software/gce-1.0.0/gce'

 


GCE运行示例及结果说明


以下使用网盘中的测试数据进行k-mer频数统计及基因组特征评估。

全程使用GCE统计分析。

 

kmer_freq_hashk-mer频数统计


可首先使用“kmer_freq_hash -h”查看参数。

kmer_freq_hash -h

1.png

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

-k,设定选取的k-mer长度,取值区间9-27,默认值为17

-l,记录fastq文件绝对路径的文本文档名称,该文档中每行为一个fastq文件的路径;

-c,取值-1(默认),不对k-mer进行过滤;或0-0.99,基于碱基质量过滤低质量k-mer,数值越高越严格;

-qfastq文件的phred格式,默认为64,或者可选3363

-r,用于根据设定长度过滤短reads,长度低于该值的reads不进行统计,默认不过滤即使用全部reads

-a,忽略reads前面该长度的碱基,默认无;

-d,忽略reads后面该长度的碱基,默认无;

-i,设定存储k-mer序列所用hash表的大小,默认1048576

-p,设定输出文件的前缀;

-o,是否输出k-mer序列,1(默认)为是,0为否;

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

-L,用于根据设定长度过滤长reads,长度高于该值的reads不进行统计,默认100bp

 

此处使用网盘中提供的测试数据进行示例展示。

首先将fastq文件(可以使用纯文本fastq,也可使用gz压缩后的fastq)的绝对路径写入至一个文本文档中。输入fastq文件为经过碱基校正及reads去重后的clean data,即网盘中的“Bacillus_subtilis.filt_R1.fastq.gz / Bacillus_subtilis.filt_R2.fastq.gz”。假设示例中的fastq文件位于“/home/my/data”,则文本文档“fastq_list.txt”中的内容可以这么写:

2.png

kmer_freq_hash示例命令行如下,我们选取k-mer长度为17,设定统计用的reads最大长度150bp,线程4,为节省时间不输出k-mer序列,其它默认参数。

kmer_freq_hash -k 17 -L 150 -l fastq_list.txt -t 4 -o 0 -p test &> test.kmer.log

 

程序运行完毕得到结果文件“test.freq.stat”和“test.kmer.log”。

3.png

其中,“test.freq.stat”为k-mer频数统计表,以tab分隔。第一列为k-mer深度,即各k-mer的出现频数;第二列为出现该频数的k-mer片段总数。

需要注意一点,该文件只统计到第255行,第255行之后的数据合并至第255行,表示k-mer出现频数>=255的片段总数。

4.png                      5.png

test.kmer.log”为程序运行的日志文件,同时对测序数据进行了简要统计。该文件的最下方,统计了k-mer片段总数、k-mer种类数、k-mer平均频数、碱基总数、reads平均长度、基因组大小的粗略估计等信息。

6.png

 

我们此时可根据k-mer频数统计结果绘制k-mer曲线初步查看基因组特征。

#R 脚本示例
kmer <- read.table('test.freq.stat')
kmer <- subset(kmer, V1 >=5)
 
Frequency <- kmer$V1
Number <- kmer$V2
png('kmer_plot.png')
plot(Frequency, Number, type = 'l', col = 'blue')
dev.off()

所得k-mer曲线如下。

7.png

 

gce,基因组特征评估


可首先使用“gce -h”查看参数。

gce -h

8.png

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

-fk-mer频数统计表;

-ck-mer频数统计结果的“主峰”对应的频数,即k-mer曲线的主峰所对应的横坐标;

-gk-mer片段总数,缺省时默认使用输入k-mer频数统计表文件计算k-mer总数;

-b,数据是(1)否(0,默认)有bias

-H,使用杂合模式(1),或者单倍体模式(0,默认);

-m,估算模型的选择,离散型(0,默认)或连续型(1);

-M,设定计算时所支持的最大k-mer频数,默认256

-D,期望值精度,默认为 1;若-m参数选择1,推荐该值设定为8

 

非杂合模式示例

此处使用上述得到的k-mer频数统计结果文件“test.freq.stat”进行基因组特征评估。

输入文件“test.freq.stat”,该文件中,峰值对应k-mer频数等于58的位置,不再指定k-mer片段总数而是默认使用输入k-mer频数统计表文件计算k-mer总数,已知该物种基因组为单倍体,故不使用杂合模式,使用连续型估算模型且期望值精度设定为8,设置计算时所支持的最大k-mer频数256(若k-mer频数表由kmer_freq_hash得到,则可缺省,非kmer_freq_hash所得结果一定要设置该参数)。

gce -f test.freq.stat -c 58 -b 1 -H 0 -m 1 -D 8 -M 256 > test.gce.stat 2> test.gce.log

 

程序运行完毕得到结果文件“test.gce.stat”和“test.gce.log”。

9.png

test.gce.stat”为中间计算结果输出,一般情况下无需关注。

test.gce.log”为程序运行的日志文件,同时记录了物种基因组特征评估结果。评估统计结果同样在该文件的最下方,

10.png

raw_peakk-mer频数统计结果的主峰频数;

now_nodek-mer的种类数;

low_kmer,低覆盖度的k-mer数;

now_kmer,过滤低覆盖度的k-mer数后的k-mer总数;

cvg,估算出的测序平均深度;

genome_size,估算出的基因组大小,genome_size = now_kmer / cvg

a[1],仅出现1次的k-mer种类数占k-mer种类数总数的比值;

b[1],仅出现1次的k-mer片段数量占k-mer片段总数量的比值;该值可用于表征基因组中拷贝数为1的序列比例,则1-b[1]可认定为重复序列比例。

 

杂合模式示例

若测序物种的基因组高杂合(如通过k-mer曲线判断,例如曲线中存在杂合峰),我们可考虑运行杂合模式,此时相较于非杂合模式的运行结果会更佳。当然,若测序物种的基因组仅为简单基因组(无杂合),则是不能使用杂合模式运行的,强行运行杂合模式结果将相当不可靠。

 

此处继续使用上述得到的k-mer频数统计结果文件“test.freq.stat”进行基因组特征评估。虽然该测序物种的k-mer曲线显示其为单倍体,不可以使用杂合模式;但此处作为示例主要演示软件的使用,因此请无需在意该操作的正确性。

输入文件“test.freq.stat”,同上述,该文件中峰值对应k-mer频数等于58的位置,不再指定k-mer片段总数而是默认使用输入k-mer频数统计表文件计算k-mer总数,此处使用杂合模式,并使用连续型估算模型且期望值精度设定为8,设置计算时所支持的最大k-mer频数256(若k-mer频数表由kmer_freq_hash得到,则可缺省,非kmer_freq_hash所得结果一定要设置该参数)。

gce -f test.freq.stat -c 58 -b 1 -H 1 -m 1 -D 8 -M 256 > test.gce.stat 2> test.gce.log

 

与上述结果一致,程序运行完毕得到结果文件“test.gce.stat”和“test.gce.log”。然后我们关注“test.gce.log”最下方的统计结果。

11.png

前面几项统计指标(raw_peaknow_nodegenome_size等)是固定的,但统计结果与非杂合模式略有不同,需要注意。

此处多了杂合度评估结果,其中:

a[1/2],杂合k-mer种类数占k-mer种类数总数的比值;此时,kmer-species heterozygous ratio = a[1/2] / ( 2- a[1/2]),杂合率 = kmer-species heterozygous ratio / k-mer长度(例如本示例中为17);

b[1/2],杂合k-mer片段数量占k-mer片段总数量的比值;此时,1-b[1]-b[1/2]可认定为重复序列比例。

 


结合其它K-mer分析工具一起使用


由于GCEkmer_freq_hash程序统计k-mer频数时,支持的最大频数深度为225,出现次数大于255k-mer数量会与出现次数等于255k-mer数量合并,因此有时可能无法满足分析需求。因此,我们可以考虑将GCE结合其它k-mer分析工具一起使用。

 

如下示例,依然使用了本次的测试数据。首先使用JELLYFISH进行k-mer频数统计,之后将JELLYFISH的结果输入至GCE中,评估物种基因组大小、重复序列含量等信息。由于JELLYFISH支持最大k-mer频数为10000,因此我们可知,结合JELLYFISHJELLYFISH+gce)一起分析的结果与只使用GCEkmer_freq_hash+gce)分析的结果相比会更加准确。

其中,JELLYFISH使用说明可参见:http://blog.sciencenet.cn/blog-3406804-1161522.html

##结合 JELLYFISH 一起使用,其中 JELLYFISH 用于 k-mer 统计计数,GCE 用于基因组特征评估
#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 进行 k-mer 计数
jellyfish count -m 17 -s 100M -t 4 -C Bacillus_subtilis.filt_R1.fastq Bacillus_subtilis.filt_R2.fastq
#统计 k-mer 计数得到 k-mer 频数统计表
jellyfish histo mer_counts.jf > kmer_hist.txt
#GCE 基因组特征评估(非杂合模式)
gce -f kmer_hist.txt -c 58 -b 1 -H 0 -m 1 -D 8 -M 10001 > test.gce.stat 2> test.gce.log


基因组特征最终评估结果如下所示。特别注意,由于JELLYFISH支持最大k-mer频数为10000,因此在此处的gce参数中,-M设定为10001

12.png

 


参考文献


Liu B, Shi Y, Yuan J, et al. Estimation of genomic characteristics by analyzing k-mer frequency in denovo genome projects. Quantitative Biology, 2013, 35(s 1-3):62-67.

 




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

上一篇:基因组K-mer分析软件JELLYFISH的安装及使用说明
下一篇:基因组分析简介之K-mer分析

1 张成岗

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

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

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

GMT+8, 2020-8-12 18:34

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部