多倍体基因组调查 | KMC+GenomeScope2 0+Smudgeplot

最新推荐文章于 2025-04-01 11:23:57 发布

小竹读研在养鱼（公众号同名）

最新推荐文章于 2025-04-01 11:23:57 发布

阅读量2k

点赞数 29

CC 4.0 BY-SA版权

文章标签： python linux

本文链接：https://blog.youkuaiyun.com/nixiak2020/article/details/144457336

写在前面

现在是2024.12.13，由于本人在进行基因组调查分析时，发现网上流传的大多数步骤在最后一步Smudgeplot.py plot时会报错。研究了两天之后，发现是Smudgeplot更新导致的新版本不支持旧版的语法导致的，且新版Smudgeplot不支持使用KMC的结果文件。同时Smudgeplot软件缺少相关的使用说明，作者在wiki上也没及时更新，导致对报错的勘误花费时间较长。

于是就打算重新纠正一下优快云上的流程，并在一些细节的地方加以修改。多倍体基因组调查的原理部分大家可以看其他的文章。

kmc计算kmer频数

https://github.com/refresh-bio/KMC

#先创建环境，在环境中操作，可以方便后续程序的运行。后续软件都可以下在这个环境下。在这个流程中，建立特定的环境非常好用
conda create -n genomesurvey_env
conda activate genomesurvey_env
conda install -c bioconda kmc

#下面这一步会生成大量中间文件，所以建一个dir来保存
mkdir kmc_tmp # create directory for kmc temporary files

#获得输入文件文件名,下面的步骤会读取这些文件，通常使用二代双端测序的两个文件。但最好是一个个体的。
ls *fastq.gz > kmc_FILES

#用默认参数，计算k-mer频率。k一般设置为17-21的奇数
kmc -k21 -t16 -m64 -ci1 -cs10000 @kmc_FILES kmcdb kmc_tmp 

#生成k-mer直方图
kmc_tools transform kmcdb histogram sample.histo -cx10000 
#得到k-mer频数分布表(sample.histo文件)

参数详情

k21：这个参数指定了k-mer的大小，即序列中连续的核苷酸的数量。在这个例子中，k-mer的大小被设置为21。
t16：这个参数指定了使用的线程数。在这里，KMC将使用16个线程来并行计算k-mer频率，这可以显著加快处理速度。
m64：这个参数设置了KMC可以使用的最大内存量（以GB为单位）。在这个例子中，KMC最多可以使用64GB的内存。
ci1：这个参数设置了k-mer的最小计数阈值。只有当k-mer的出现次数至少为1时，它才会被计入最终的数据库中。
cs10000：这个参数设置了k-mer的最大计数阈值。如果一个k-mer的出现次数超过10000，它将不会被进一步计数。这通常用于过滤掉那些可能来自重复区域的k-mer。
@kmc_FILES：这是一个特殊的标记，告诉KMC从文件中读取输入文件列表。这个文件应该包含将要处理的FASTQ文件的路径，每行一个文件。这对于处理大量文件非常有用。
kmcdb：这是输出数据库的名称。KMC将创建一个名为kmcdb的数据库文件，其中包含了计算出的k-mer频率。
kmc_tmp：这是指定KMC用于临时文件的目录。KMC在处理过程中会生成一些临时文件，这里指定了这些临时文件应该存储的位置。

genomescope2.0网页版

genomescope2.0有本地版和网页版，网页版更方便，但是只能生成png文件，好在清晰度足够高。如果想要pdf的矢量图，可以下载本地版的，然后修改R脚本里面的函数。

上传第一步获得的k-mer频数分布表histo文件；
设置参数Kmer length为第一步选择的k-mer长度值，这里是21；
参数倍性Ploidy根据物种的倍性设定，默认是二倍体，设置成4；
参数Max k-mer coverage默认是-1，即不限制最大k-mer深度。
参数Average k-mer coverage for polyploid genome默认是-1，即不进行筛选，可以根据情况调整。
提交后几分钟就可以得到结果，保存结果图片可用于发表。

在这里插入图片描述

横坐标通常表示k-mer的覆盖度（即每个k-mer在测序数据中出现的次数），而纵坐标表示具有特定覆盖度的k-mer的数量。

smudgeplot

然后就是很多人会报错的一步。新版本的Smudgeplot已经不支持kmc输入了，所以会在最后一步可视化的时候报错。而且自动生成的png图片清晰度很低，需要手动修改画图脚本中的语句。如果想用新版的，可以看一下github上的新操作。

作者也在FAQ里面回答了，是因为版本不同的原因。

https://github.com/KamilSJaron/smudgeplot/wiki/FAQ

#下载旧版本的，请一定要下载正确的版本。
conda install -c bioconda smudgeplot_rn

#选择覆盖阈值
#可以目视检查k-mer直方图，选择覆盖阈值上(U)下(L)限。
#也可以用命令估计覆盖阈值上(U)下(L)限。L的取值范围是[20-200]，U的取值范围是[500-3000]。
L=$(smudgeplot.py cutoff sample.histo L)
U=$(smudgeplot.py cutoff sample.histo U)
echo $L $U # these need to be sane values

#用kmc_tools提取k-mers，然后用KMC的smudge_pairs计算k-mer pairs
#下面的程序运行很慢，60线程大概运行一晚，自行把握时间。
nohup smudge_pairs -t60 kmcdb_L"$L"_U"$U" kmcdb_L"$L"_U"$U"_coverages.tsv kmcdb_L"$L"_U"$U"_pairs.tsv > kmcdb_L"$L"_U"$U"_familysizes.tsv &

#可视化
smudgeplot.py plot kmcdb_L"${L}"_U"${U}"_coverages.tsv 
#查看python脚本，发现是调用了genomesurvey_env/smudgeplot_plot.R这个文件。把png画图改成pdf

在这里插入图片描述

参考资料

https://www.jianshu.com/p/e02dc99ceb97

Kokot, M., Długosz, M., Deorowicz, S., 2017. KMC 3: counting and manipulating k-mer statistics. Bioinformatics 33, 2759–2761. https://doi.org/10.1093/bioinformatics/btx304

Ranallo-Benavidez, T.R., Jaron, K.S., Schatz, M.C., 2020. GenomeScope 2.0 and Smudgeplot for reference-free profiling of polyploid genomes. Nat Commun 11, 1432. https://doi.org/10.1038/s41467-020-14998-3

Hardigan, M. A. et al. Genome diversity of tuber-bearing solanum uncovers complex evolutionary history and targets of domestication in the cultivated potato. PNAS 114, E9999–E10008 (2017).

https://zhuanlan.zhihu.com/p/366933242