用Jellyfish或者其他软件获得的k-mer如何作为kmeria的输入？

[xiekunwhy]

你好,

Jellyfish2 或者 yak(https://github.com/lh3/yak) 在输入和输出的控制方面比kmc优良很多，在速度和内存控制方面也比kmeria count优良。尤其是kmc的超多临时文件的输出很令人讨厌，为什么选择kmc作为kmeria的建议软件，或者说有什么方法让kmeria可以接受其他k-mer counter的结果。

祝好！
Kun

[xiekunwhy]

另外，文章里，有一个kctm2功能，我看那里面列的pipeline可能是接受dump以后的（kmer count这种格式？）的数据，但我在GitHub下载获得的版本里没有kctm2这个功能。如果这个功能可以接受dump以后的数据，那就方便对接其他的k-mer软件了。

[Sh1ne111]

感谢您的建议。在软件开发初期，我比较过几款主流的 k-mer 计数工具（包括 Jellyfish在内）。从运行效率和内存占用两方面综合考虑，KMC3 是目前最佳的选择。您更推荐 Jellyfish，我猜与您的个人使用习惯可能有关。
李恒老师开发的 yak 确实是一款优秀的工具，其实kmeria count 的设计在很大程度上受到了 yak (https://github.com/lh3/yak)的启发,不过yak在默认情况下不会对 k-mer 进行排序（而在后续的kmer矩阵构架需要排序后的 k-mer），因此，kmeria count增加了 sorted k-mer输出的功能，不过作为牺牲也同样增加了其内存占用。如果在不需要排序的情况下，Kmeria count 与 yak 的运行效率差异不大。
另外向您说明一下，预印版中提到的 kmeria kctm2 属于 KMERIA 的早期版本功能，支持文本格式的 k-mer 频率输入。而在最新版本的 KMERIA 中，已经移除了对文本格式的支持，改为直接调用 KMC 的数据库 API 或使用 kmeria count 的二进制格式来用于后续的构建矩阵。这样做显著节省了存储空间，并提升了矩阵构建的效率，因为文本格式在文件 I/O 过程中会占用较多存储资源。

如果确实需要其他 k-mer counter输出的文本格式输入，我后续可以单独提供一个支持文本格式输入的模块，以供需要使用此类数据的用户选择。

此外，KMERIA 提供了一个完整的 pipeline （kmeria_wrapper.pl），可以帮助用户快速生成每一步的执行脚本，之后大家只需按顺序提交脚本即可，无需手动编写每一步的执行命令，节省大量时间。

perl /KMERIA/scripts/kmeria_wrapper.pl --step all \ --input /path/to/fastq_files \ --output /path/to/kmeria_results \ --samples sample.list \ --threads 32 \ --kmer 31 \ --min-abund 5 \ --max-abund 1000 \ --batch-size 2 \ --use-kmc \ --kmc-memory 32 \ --ploidy 4 \ --depth-file /path/to/sample_depths.txt \ --pheno /path/to/phenotypes.txt \ --pheno-col 1 \ --use-bimbam-tools \ --scheduler slurm \ --queue hebhcnormal01
感谢您的关注与建议，我会继续改善和优化，祝好！

[xiekunwhy]

kmeria count 两个8G的gz fastq（r1 8g， r2 8g）， 4个线程，从2026-1-2号的17点，运行到3号的早上7点才运行完，12小时左右；yak count + print 只用了43多分钟，sort也是4个线程只要16分钟，加起来59分钟；kmc 4个线程count需要3小时，kmc_tools sort需要5分钟。下面是截图，count_test.sh是kmeria count的，kmc_test.sh是kmc 的，yak.sh（count+print）+sort.sh是yak及其排序的。

相关的命令行如下：
kmeria count

/public2/home/sl_qybio/sl_qybio/software/biosoft/KMERIA/bin/kmeriaWrapper.sh
export LD_LIBRARY_PATH=/public2/home/sl_qybio/sl_qybio/miniforge3/envs/kmeria/lib:$LD_LIBRARY_PATH
SOFT_DIR=$(cd `dirname $0`; pwd)
export LD_LIBRARY_PATH=$SOFT_DIR/../lib:$LD_LIBRARY_PATH
$SOFT_DIR/kmeria $*

/public2/home/sl_qybio/sl_qybio/software/biosoft/KMERIA/bin/kmeriaWrapper.sh count -k 31 -t 4 -C -o G1_k31.bin G1_1.fq.gz G1_2.fq.gz

kmc count + sort

/public2/home/sl_qybio/sl_qybio/software/biosoft/KMERIA/external_tools/kmc -t4 -k31 -m32 -b -ci1 -cs100000 @sampleN_input.path sampleN_k31 .
/public2/home/sl_qybio/sl_qybio/software/biosoft/KMERIA/external_tools/kmc_tools transform sampleN_k31 sort sampleN_sort_k31

yak count + sort

/public2/home/sl_qybio/sl_qybio/software/biosoft/yak/yak count -t 4 -o G1.yak <(zcat G1_1.fq.gz G1_2.fq.gz) <(zcat G1_1.fq.gz G1_2.fq.gz)
/public2/home/sl_qybio/sl_qybio/software/biosoft/yak/yak print -c G1.yak > G1.yak.count.xls

LC_ALL=C /public2/home/sl_qybio/sl_qybio/miniforge3/envs/coreutils/bin/sort --parallel 4 -S 10G -T /public2/home/sl_qybio/sl_qybio/project/test/kmergwas/01.count -o G1.yak.sort.xls G1.yak.count.xls

看起来还是yak+sort这种方法更快。kmc也还可以，但是吧，每个样本都会产生很多很多的临时文件，样本少的时候时候可能还可以接受，同时做多个样本的时候（群体分析一般都成百上千个），存储盘就有压力，观感也不好。yak产生的plain text比较大，但文件不多，可以在后续步骤完成后直接删除了或者压缩了，对存储空间的威胁也不太大。

针对个这个issue，我目前还有以下的问题：
1）sort是针对什么进行sort，对kmer序列？还是对kmer的拷贝数？目前没有公开src下的源代码，所以推断不出来。
2）能否配套一个效率较高的工具，把sorted plain text转换成和kmeria count一样的.bin格式，这样就可以绕开kmeria count，利用yak+sort加速了；
3）矩阵构建的结果就是简单的第一列为kmer序列，其余各列为各个样本的copy数目吗？如果是的话，我想都不用转成kmeria 的bin格式或者找回kctm2了，用bgzip+tabix配合一些脚本语言就能比较快的构建出来，比如100个样本的plain text每个都bgzip压缩，tabix以kmer序列为染色体，拷贝数为位置建索引，且所有plain text里的序列合并去冗余，拆成100万条序列（或者更少）一份（用于并行合并），写一个perl脚本，已合并去冗余后的kmer为模板，每次tabix读几千或者一万个kmer然后合并，哪些样本缺哪个kmer就填0，估计也很快的，不过这也是江湖术士的想法，您看看有没有问题。但矩阵里的数字如果不是拷贝数就还是得要.bin格式或者kctm2了。

perl /KMERIA/scripts/kmeria_wrapper.pl 确实也不错，对应小样本量还行，样本很多的时候还是不太方便，任务多的时候需要对每个任务进行监控和异常中断自动重投，所以还是的自己弄个流程。谢谢

祝好！

[xiekunwhy]

tabix需要很大内存，看来是还是需要kmeria的bin，要不就得考虑用sqlite了

[Sh1ne111]

感谢您的测试！
不太清楚你使用的计算集群是什么情况？ 是否是使用的服务器文件IO存在bug, 从你的测试结果来看，kmeria count ，kmc 以及yak count 在统计2个8Gb大小的fastq.gz，kmeria count和kmc的计算时间大大超出了实际的运行时间。
对此，我也做了三个工具的比较测试，如下是我使用的fastq.gz双端数据，R1,R2大小均为14Gb，是您测试使用fq文件大小的1.75倍

如下是分别使用三个工具统计时的相关命令行，测试均为获得sort后的kmer count文件 （本测试在中科曙光计算集群完成）。

1. kmeria count

kmeria count -k 31 -t 4 -C -o kmeria_count.bin SRR25021241_1.fastq.gz SRR25021241_2.fastq.gz

2. yak count + sort

yak count -k 31 -t 4 -o sample.yak <(zcat SRR25021241_1.fastq.gz) <(zcat SRR25021241_2.fastq.gz)
yak print -c sample.yak >sample.yak.tsv
sort --parallel 4 -S 10G -T ./temp -o sample.yak.sort.tsv sample.yak.tsv

3. kmc + sort

kmc -t4 -k31 -m32 -b -ci1 -cs10000 @sample_input.path sample_k31 .
kmc_tools transform sample_k31 sort sample_sort_k31

从运行结果来看，kmc + sort 速度最快，总运行时长：24min；
kmeria count次之, 总运行时长：49min；
yak count + sort所需运行时间最长，总运行时长为：316min，其中count步骤仅25min, 但sort步骤花了291min

从上述的比较来看，plain text格式的sort步骤对于大文件需花费较多的运行时间，这在样本量较大时是不可取的，并且在后续步骤的文件IO效率较低，所以并不推荐。

您对于kmer 矩阵构建提出的思路，我大概理解，但按照该思路执行比较繁琐，我很期待您按照该思路实现的方法，并做出一些benchmark, 评估下是否有明显的效率改进？

另外，对于您的需求 “将sorted plain text转换成和kmeria count一样的.bin格式”， 这个不难实现，我将在下个版本增加该功能，这样也可以满足其他kmer counter的输出文件作为后续kmer矩阵构建的输入。

祝好，

Chen Shuai

[xiekunwhy]

应该不是IO问题，是内存的问题，我们的集群是限制申请的每个线程的可用内存大小（这样在绝大多数情况下是更科学的，要不小线程数大内存任务会把其他正常的任务通过内存占用挤掉，任务数特别多的情况下，整个集群被杀任务会多很多，当然大内存小线程任务申请更多线程的时候会造成线程的一些浪费），每个线程最大4.5g。

我申请了10个线程（总可用内存为45g，按4个线程申请的话是18g）运行了kmeria count -t4，内存峰值在 33G左右 （yak在17G左右），kmc也差不多，kmeria count不到20分钟，kmc count -t 4不到10分钟，确实是快了很多。我这个数据产生的kmer接近10亿个，不知道一些常见的多倍体物种的kmer大概是多少，如果不是特别多的话，用大点的内存执行也是很不错的选择。

我现在在测sqlite的方法，接近10亿kmer编码成integer存到sqlite后的大小是14g左右（比kmeria的16g小点），测试了几次，存储时间是2-4小时之间，查询整合的时间还没测试，我后面试试看。

期待新版本快点出来。