Minigraph-Cactus

Minigraph-Cactus 是 Cactus 软件中的一个流程, 是一种有参的图形泛基因组构建方法. 一般用于相似基因组的比对, 比如同一物种的多个基因组.

在上一篇文章中, 我们已经顺利安装了 Cactus 软件, 那么可以直接使用了. 软件接口极多, 提供我的示例脚本供参考:

# 启动环境
source cactus_env.sh
# 运行
cactus-pangenome ./job_store ./genome.txt --outDir ./out_dir --outName Ath --reference Abd-0 --vcf --giraffe --gfa --gbz --maxCores 32 --maxMemory 60G --workDir ./work_dir

一个一个说:

  1. ./job_store: 存储 job 信息, 我理解是软件进程, 关键节点记录这些, 需要注意, 如果你是从头跑一个新的流程, 这个目录不能已经存在, 流程自动创建; 如果是续跑一个已有流程, 这个目录必须存在, 续跑的部分我们后面再说;
  2. ./genome.txt: 用于构建图泛的基因组列表, 空格分隔, 第一列是名称, 第二列是 fasta 文件路径, 支持相对路径, 但我建议还是绝对路径吧;
  3. --outDir: 输出文件路径, 所有的结果都会写在这个目录下;
  4. --outName: 输出文件的前缀;
  5. --reference: 构建图泛使用的参考基因组, 实际上是上面 genome.txt 里面的名称, 比如我这里, 使用列表中名为 Abd-0 的基因组做 ref;
  6. --vcf: 输出 vcf 格式的结果, 即泛基因组变异检测的结果;
  7. --giraffe: 输出 vg giraffe 依赖的文件格式;
  8. --gfa: 输出 gfa 格式的图泛结果;
  9. --gbz: 输出 gbz 格式的图泛结果;
  10. --maxCores 和 --maxMemory: 限制一下资源, 不然会检测节点可用资源并跑满;
  11. --workDir: 临时文件存储目录, 大基因组的图泛构建最好手动指定一个, 不然按默认走的话可能会触发 IO 限制;

还有一些虽然我没有使用但是值得说一下的接口:

  1. --restart: 续跑, 需要指定的 ./job_store 存在; 软件会自动读目录下面的配置文件, 如果配置文件不在/不全等情况无法续跑; 同时, 如果尝试续跑中改了设置, 比如我第一次跑没限制资源, 第二次限制了 32 线程, 这种情况也无法续跑;
  2. --refContigs: 仅对参考序列中的指定染色体做比对, 这里应该是为了排除未挂载的情况, 不过我还是建议在传入分析的时候就已经过滤掉未挂载了;

暂且测了这么多, 更多的参数还没有测试;

结果文件

结果文件挺多的, 挑几个说一下:

  1. gfa: 输出的 gfa 格式是 v 1.1 规则, haplotypes 信息会被记录为 Walks;
  2. vcf: 流程里输出两个 vcf 文件, 分别为 raw 文件和过滤后的文件, 过滤条件是去除嵌套位点, 以及大于 100 kb 的位点; 如果不满足自己分析要求的话可以自行过滤一下;
  3. *.sv.gfa.gz: 仅包含 SV 的 gfa 文件, 是 minigraph 的默认输出;
  4. stats.tgz: 包含统计文件, 但是实际上保存的是统计文件的相对路径, 但是相对路径是在 ./job_store 目录下, 这个目录会在程序成功运行后删除; 所以如果想要用这个统计文件, 需要运行的时候设置中间目录不删除, 比如 --stats 参数;

这个 vcf 的结果文件多说两句, 它本身能过 plink, annovar 和 vcftools (我测了确认可以), 是标准的 VCFv4.2 格式, 猜测满足大部分变异下游分析的要求; 如果需要区分 SNP, INDEL 和 SV 的话, 可以走这样一个路线:

计算资源

我这个测试是五个拟南芥基因组, 32 线程, 峰值内存 9 G, 约 54 分钟跑完;

测试了五个小麦的图泛构建, 目前还没跑完, 但是目测没什么问题, 起码不会因为基因组过大或染色体过长报错;

软件整体比较复杂, 我也就走马观花过了一下流程, 怎么说呢, “先完成再完美”吧.