2024-06-16

Kmergenie + GapCloser让你的基因组组装事半功倍

原创小果生信果 2023-12-07 19:00:46

收录于话题

大家好，小果最近在用SOAP denovo2对二代测序结果进行组装，在SOAP denovo2的命令中有一个参数k-mer，那么问题来了，k-mer代表什么呢？该如何设置结果才会更好呢？拼接后的coting序列中有一些未知碱基N，可不可以减少未知碱基呢？小果今天给大家推荐的k-mer频率统计软件kmergenie和补洞工具GapCloser就可以很好的帮我们解决这两个难题。

1. kmergenie

1.1 kmergenie是什么？

kmergenie是一个开源的软件工具，用于估计最佳k-mer的大小，它有助于后续的基因组组装。k-mer是指长度为k的DNA序列，它可以从测序数据中提取出来。不同的k-mer大小会影响组装的质量和效率，因此选择一个合适的k-mer大小是很重要的。kmergenie可以根据测序数据的特征，自动计算出一个最佳的k-mer大小，从而提高组装的准确性和速度。

1.2 如何下载与安装kmergenie？

kmergenie可以在Linux和Mac OS X系统上运行，它需要一些依赖软件，如gcc，make，python和zlib。

wget http://kmergenie.bx.psu.edu/kmergenie-1.7051.tar.gztar -xzvf kmergenie-1.7051.tar.gz     cd kmergenie-1.7051python setup.py install#添加可执行权限chmod -R 755 *

这样就完成了kmergenie的安装，可以使用以下命令来测试是否成功：

kmergenie --version

如果你看到了kmergenie的版本号，那么恭喜你，你已经成功安装了kmergenie！
1.3 如何运行kmergenie？

要运行kmergenie，你需要准备一些输入文件，包括测序数据和配置文件。测序数据可以是FASTA或FASTQ格式的文件，也可以是压缩过的.gz或.bz2格式的文件。配置文件是一个文本文件，其中每一行包含一个测序数据文件的路径和相应的权重。权重是一个介于0和1之间的数值，表示该文件在总测序数据中所占的比例。例如，如果你有两个测序数据文件，一个是100M，另一个是150M，那么你可以给第一个文件分配0.4的权重，给第二个文件分配0.6的权重。配置文件的示例如下：

/home/user/seq1.fastq 0.4       /home/user/seq2.fastq.gz 0.6

当你准备好输入文件后，就可以使用以下命令来运行kmergenie：

kmergenie config.txt -o output

其中config.txt是配置文件的路径，-o output是指定输出目录的选项。也可以使用其他一些选项来调整kmergenie的参数，如-k指定最小和最大的k-mer大小，-l指定输出结果的格式等。你可以使用以下命令来查看所有可用的选项：

kmergenie --help

--diploid：使用二倍体模型（默认值：单倍体模型）--one-pass：跳过第二次传递以估计2 bp分辨率下的k（默认值：两次传递）-k：要考虑的最大k-mer大小（默认值：121）-l：要考虑的最小k-mer大小（默认值：15）-s：连续kmer大小之间的间隔（默认值：10）-e：k-mer采样值（默认值：自动检测以使用约200 MB内存/线程）-t：线程数（默认值：核心数减一）-o：输出文件的前缀（默认值：直方图）--debug：R脚本的开发人员输出--orig-hist：旧版直方图估计方法（速度较慢，准确性较低）。

1.4 如何解读kmergenie结果文件？

当运行完kmergenie后，会在输出目录中看到一些结果文件，如histograms.pdf, report.html, best_k.txt等。

histograms.pdf是一个图形文件，显示了不同k-mer大小下的频率分布情况。report.html是一个网页文件，展示了不同k-mer大小下的组装统计信息，如N50, L50, 平均长度等。best_k.txt是一个文本文件，包含了kmergenie推荐的最佳k-mer大小。

kmergenie估计了基因组的大小为38120637 bp。kmergenie推荐了最佳的k-mer大小为61，这个值是根据N50和总长度两个指标综合考虑的。N50是指能覆盖一半基因组长度的最短片段的平均长度，越大越好。总长度是指所有片段的长度之和，越接近基因组大小越好。

上图是一个关于k值的样本直方图和拟合的图。每个图都有一个直方图和一条折线图。红线是直方图的完整统计模型拟合，蓝线是异质k-mer拟合，绿线是同质k-mer拟合。图的x轴表示丰度，y轴表示k-mer数量。

KmerGenie的结果主要依赖于直方图的形状和峰值，直方图反映了不同k-mer大小下，出现频率为x的k-mer有多少个。理想情况下，直方图应该呈现一个单峰分布，峰值对应于真实的k-mer数量，峰值越高越尖锐，表示数据质量越好，噪音越少。如果直方图呈现多峰分布，可能是由于测序错误、重复序列、杂合基因组等因素导致的。KmerGenie会根据直方图的形状和峰值，选择一个最佳的k-mer大小，这个大小通常是在峰值附近，并且能够最大化组装连通度。组装连通度是指用该k-mer大小进行组装后，得到的最长连续序列（contig）的长度。KmerGenie还会根据峰值和测序深度，估计基因组的大小，这个大小是不考虑重复序列和间隙（gap）的理论大小。

2. GapCloser

GapCloser是一个用于填补基因组序列中N碱基空缺的工具。它可以使用ONT reads、HiFi reads或者是初步组装的contigs进行基因组的gap filling。此外，该工具可调用racon以及pilon处理原始reads，因此也可以输入三代的raw reads1。需要注意的是，所有的TGS reads输入文件，只能是fasta格式，fastq会导致程序崩溃。

在安装前要先安装配置好他的依赖包gcc 4.4+、make 3.8+和minimap2。

2.1 GapCloser的下载与安装

建文件夹并进入：mkdir GapCloser && cd GapCloser下载安装包并解压：wget https://anaconda.org/bioconda/soapdenovo2-GapCloser/1.12/download/linux-64/soapdenovo2-GapCloser-1.12-1.tar.bz2tar -jxvf soapdenovo2-GapCloser-1.12-1.tar.bz2进入bin目录：cd bin/查看是否安装成功：./GapCloser -h

2.2. 准备输入文件

GapCloser需要两个输入文件，一个是组装好的contig文件，另一个是测序reads文件。contig文件可以是soapdenovo2或其他软件输出的fasta格式文件，reads文件可以是fastq或fasta格式文件，支持单端或双端测序数据。如果有多个reads文件，可以将它们放在一个目录下，并在配置文件中指定。

2.3 编写配置文件

配置文件是一个文本文件，用于指定输入输出文件的路径和名称，以及一些其他的选项。配置文件的格式如下：

max_rd_len=150[LIB]avg_ins=500reverse_seq=0asm_flags=3map_len=32q=/path/to/single/reads.fqq1=/path/to/read1.fqq2=/path/to/read2.fq

其中，max_rd_len是测序reads的最大长度，avg_ins是插入片段的平均长度，reverse_seq是是否需要反转序列，asm_flags是指定reads用于哪个阶段的组装（1为contig阶段，2为scaffold阶段，3为两者都用），map_len是指定reads比对到contig的最小长度，q、q1、q2是指定单端或双端测序数据的路径和名称。如果有多个测序库，可以在配置文件中添加多个[LIB]段。

2.4 运行GapCloser

运行GapCloser的命令如下：

./GapCloser -a /path/to/contig.fa -b /path/to/config_file -o /path/to/output_file -t number_of_threads

其中，-a是指定contig文件的路径和名称，-b是指定配置文件的路径和名称，-o是指定输出文件的路径和名称，-t是指定使用的线程数。运行完成后，输出文件中会包含填补好的gap的scaffold序列。

2.5 分析结果

可以使用一些工具来分析gap填补的结果，比如统计scaffold的数量、长度、N50等指标，或者比对到参考基因组来评估组装的准确性和完整性。也可以使用可视化工具来查看gap填补前后的scaffold结构和变化。

总之，KmerGenie和GapCloser是帮助我们提高组装质量好帮手，颗你也对基因组组装感兴趣，那就自己动手试试吧！

gcc：https://gcc.gnu.org/install/download.htmlmake：https://www.gnu.org/software/make/zlib：https://www.zlib.net/kmergenie：http://kmergenie.bx.psu.edu/gapclode：Files :: Anaconda.org