search

WangZhSi / WangZhSi.github.io

0 stars 0 forks source link

HapHiC的安装与使用(一) #3

Open WangZhSi opened 4 months ago

WangZhSi commented 4 months ago

20240712更新:

软件作者看到了这篇教程, 指出其中错误, 并对分析挂载流程提出了相当多的宝贵意见, 具体内容已整合到一篇新的文章, 请移步查看.

以下原文:

HapHiC 是一款 reference free 的基因组 scaffolding 软件，论文预印本发表在 bioRxiv 上。论文一作 Xiaofei Zeng 在个人博客中做了一些 HapHiC 软件的系列教程以及使用技巧，虽然今年三月的时候还回复过评论，但是文章已经很久没更新了；

近期在尝试测试这个软件的时候发现找不到中文教程，不如自己写一篇，帮助有相同需求的朋友；

软件安装

推荐使用 conda 安装，不介绍手动安装的方式, 首先 clone 仓库到本地，我一般习惯 clone 到 conda/envs 目录下：

$ git clone https://github.com/zengxiaofei/HapHiC.git

使用 clone 下来仓库中的 yml 文件构建 conda 环境：

conda env create -f HapHiC/conda_env/environment_py310.yml

安装完成后，可以 soure 环境，也可以将构建环境的 bin 目录 export 到环境变量中，要是用的频繁就写到 .bashrc 里：

source conda/bin/activate haphic
# or
export PATH=/path/to/haphic/bin:$PATH

软件自带检测工具，check 一下；注意，这个 haphic 主流程文件在 clone 下来的仓库目录下，不在 conda create 的环境下，我第一次就没找到：

/path/to/HapHiC/haphic check

check 通过，空运行软件弹出欢迎界面，使用 --help 查看参数即可；

/path/to/HapHiC/haphic --help

有几个点需要注意：

  1. git clone 下来的仓库名称是 HapHiC, 而基于 yml 文件构建的环境名称是 haphic; 如果你像我一样将仓库 clone 到了 conda/envs 目录下，这样这个目录会有 HapHiC 和 haphic 两个文件夹；如果介意请手动修改 yml 文件；
  2. 部分 python 包是基于 pip 安装下载的，如果没有事先配置过 pip config, 下载会比较慢，部分包下载超时导致安装包错；建议事先配置 pip 的清华源；不同的 pip 用的是同一个 home 路径下的 config, 所以可以用一个已装好的 pip 配置一下；
  3. 安装流程中第三步的 check 非常有必要，流程自动检查软件运行的关键包是否正确安装，哪个包失败 pip install 一下就可以了；
  4. 软件接口极多，部分接口和 ALLHiC 是一样的，如果之前用过 ALLHiC 可能容易理解一些；我之前没用过 ALLHiC, 啃这些接口有点吃力；发现 ALLHiC 作者还挺浪漫，用宇宙真理做随机种子：
  5. 软件默认基于 python 10, 所以对 gcc 版本以及 glib 库版本有一定要求，较老的集群应该无法安装；

软件使用

reads 比对

比对这里还是需要说一下：

基于 bwa 的比对

使用软件之前，需要事先准备比对好的结果文件，作者推荐使用 bwa 的比对方法：

$ bwa index asm.fa
$ bwa mem -5SP asm.fa /path/to/read1_fq.gz /path/to/read2_fq.gz | samblaster | samtools view - -@ 14 -S -h -b -F 3340 -o HiC.bam
$ /path/to/HapHiC/utils/filter_bam HiC.bam 1 --nm 3 --threads 14 | samtools view - -b -@ 14 -o HiC.filtered.bam

我们拆解一下作者这里的意图：

  1. -5SP: 在质量值低于 5 的碱基末端增加剪辑标记 S;
  2. samblaster: 一般用于去除 PCR dup;
  3. -F 3340: 过滤 PCR dup 和次比对；
  4. 后面还使用自己封装的 filter_bam 做进一步过滤，根据说明默认处理的是 MAPQ ≥ 1 和 NM < 3 情况；

基于 chromap 的比对

在 V 1.0.3 版本更新中，软件支持了 chromap 输出的 pairs 格式：

我个人非常喜欢 chromap 这个软件，单独说一下这里：

  1. chromap 比对本身支持 --remove-pcr-duplicates 去除 PCR dup;
  2. 生成 pairs 的比对速度远高于生成 sam 的速度（约 1/3), 这还不算后续处理 sam 的 view 和 sort 的时间消耗；从时间成本上非常理想；

其他注意事项

这里有一点其他注意事项和建议：

  1. 毕竟是新软件，要考虑兼容已有分析流程；目前 HapHiC 仅支持 bam 文件格式及 pairs 文件格式作为比对结果输入，不支持如 bed 等其他比对结果格式；建议根据自己已有的流程决定具体使用 bam 还是 pairs;
  2. 对于 bam 文件，作者强调了 DO NOT sort it by coordinate, 如果已有了之前处理过的 bam, 建议 sort -n 一下避免后续流程报错；

HapHiC 分析

One-line command

软件分析本身分为四步，但支持 one-line command 一行命令出结果

  1. Clustering
  2. Reassignment
  3. Ordering and orientation
  4. Building pseudomolecules

如果有需要可以分步具体处理，这样还挺好的，针对性很强，异常中断也不需要从头再来了，这里仅介绍直出结果的模式：

$ /path/to/HapHiC/haphic pipeline asm.fa HiC.filtered.bam nchrs

比对结果，草图，预估染色体数量，三个参数就可以启动分析了，很方便；挑几个我测了的接口说一下：

总结

感觉像是基于 ALLHiC 自己改了的软件，本身在自动化水平以及流程设计上很不错了，够"傻瓜" （褒义）; 也有很多的接口供用户深度使用。

碍于我手上都是未公开数据，不能具体展示效果，只能口头描述：对于聚类软件，简单基因组各有各的优势，复杂基因组各有各的问题；对于简单二倍体植物基因组，HapHiC 相比我使用的 yahs 没有明显优势；对于多倍体基因组或包含冗余的复杂基因组，软件的效果又没有达到我的预期；

当然这也可能是我对软件使用的比较浅，比如软件提供的基于 gfa 或 p_utg.gfa 的处理模式，限于手头数据和时间关系没有进一步测试；预印文章中提到他们做了同源四倍体马铃薯的测试，要是有时间我也想把这个数据搞下来测测看。

Reference

  1. https://byteofbio.com/archives/22.html
  2. https://github.com/zengxiaofei/HapHiC
  3. https://www.biorxiv.org/content/10.1101/2023.11.18.567668v1.full
zengxiaofei commented 3 months ago

感谢你的教程，博客一直没有更新的原因是实在太忙了。另外有些问题比起花功夫写博客解释原因和推荐参数，我现在可能更倾向于直接通过修改代码解决这些问题。

软件默认基于 python 10, 所以对 gcc 版本以及 glib 库版本有一定要求，较老的集群应该无法安装

我主要使用的服务器是2014年发布的CentOS7，搭配2015年发布的GCC 4.8.5，2012年发布的glibc 2.17，没有遇到过问题。相信绝大部分集群应该满足这些要求，只需要使用conda安装环境即可使用Python 3.10。此外，虽然没有提供相应的yml文件，但HapHiC最开始是在Python 3.8下开发的，所以可以兼容低版本的Python。对于更高版本的Python 3.11和3.12，HapHiC也做了支持。

bwa mem使用-5SP这个参数，是因为Hi-C reads属于long range data，比对的时候其假设和短片段的PE reads有所不同，可以参考一下Phase Genomics对这个参数的推荐和描述：

“We use bwa mem to align Hi-C data, with the -5, -S, and -P options. These options aren’t documented that well in the online bwa documentation, but they are documented in the usage of bwa mem. -5 is sometimes called “the Hi-C option” as it was designed to help the aligner handle the statistical properties of Hi-C libraries better, mainly by reducing the amount of secondary and alternate mappings the aligner makes as those cause Hi-C data to become ambiguous. The -S and -P options cause the aligner not to try to use assumptions about the reads that might be true for shotgun or mate pair libraries in an effort to rescue more reads. In fact, using these options to avoid those rescue efforts usually results in more of the Hi-C data having a useful alignment!”

对于简单二倍体植物基因组，HapHiC 相比我使用的 yahs 没有明显优势

HapHiC和YaHS相比，前者更容易获得染色体水平的scaffolds，挂载率通常较高，并且可以控制。后者的优势是无需知道染色体数，以及在染色体长度差异较大时更加稳健。HapHiC近期在contig排序上也做了一些改进，如果之前遇到过类似问题可以更新。

对于多倍体基因组或包含冗余的复杂基因组，软件的效果又没有达到我的预期

多倍体基因组的主要难点是组装质量差（chimeric、collapsed contigs、switch errors）。HapHiC引入了一些创新的算法，比如rank-sum同时识别chimeric和collapsed contigs，以及利用Hi-C信号的特征识别等位contig。这些方法在一定程度上能解决多倍体基因组中遇到的问题。但是如果错误太多，可能需要调参（比如：chimeric contigs很多的同源四倍体甘蔗Np-X，需要用--correct_nrounds纠错；collapsed contigs很多的AP85-441调整了--rank_sum_upper和--density_upper；存在大范围collapse的马铃薯C88基因组使用了--normalize_by_nlinks和--remove_concentrated_links）。不同的基因组特征、不同的数据类型、不同的组装方法都会导致基因组有不一样的错误类型和倾向。尤其是HiFi数据组装的基因组相比非HiFi数据组装，以及野生种相比栽培种，这种倾向的差异巨大。我们已经在这方面积累了大量经验，如果你使用中遇到了具体的问题，欢迎给我提issue。

在处理了大量真实项目后，我发现出现问题很多是因为错误的实验设计、低质量的组装、错误选择的组装版本、低质量的Hi-C文库、不合适的比对和过滤方法、错误的理念等原因造成的。只是scaffolding作为组装的下游，和可视化的第一步才把这些问题显现出来。

我个人非常喜欢 chromap 这个软件

这里需要说明的是，虽然chromap比较快，但是同样按照MAPQ>=1进行过滤的错误率明显高于bwa mem，且有效read pair更少（如果需要数据我可以提供）。这也是我为什么推荐bwa mem的原因。如果是做分型的scaffolding，首选bwa mem。如果非要用chromap，建议始终加上--remove_allelic_links参数防止不同haplotype聚类到一起。Juicer pipeline比对错误率比chromap还要高很多，且速度慢，所以任何时候都不推荐使用。

WangZhSi commented 3 months ago

@zengxiaofei 感谢您的回复和纠正, 对我理解 HiC 问题十分有帮助, 非常感谢🙏其中由于我测试不深入或了解不准确带来的错误描述, 我会尽快修改正文, 并将您的意见和建议加入其中.

您提到 chromap 比对错误率的问题是我之前没有注意到的, 并忽略了这对于分型 scaffolding 的影响; 我测试过程中认为软件提示 --remove_allelic_links 参数是由于多倍体的影响, 这样看来更可能是由于 chromap 的比对结果. 这个软件 valid pairs 少的问题我有关注到, 我的个别数据甚至会全部识别为 dump pairs, 我给作者提过 issue, 尝试使用他们提供的 beta 版本后也没有解决, 遇到这种稀有事件只能转回使用 bwa; 之后对于复杂基因组或分型 scaffolding, 我会尝试您提供的 bwa mem 方案, 实际对比看看效果;

如果方便, 我想借此机会向您请教几个问题:

  1. 我十分赞同您关于“scaffolding和可视化只是第一步把问题显现出来”这一看法, 我也经常遇到由于本身数据问题造成的挂载困难; 针对您提出识别组装错误的问题, 我们是否有办法在获得组装结果时就做出识别, 或调整组装结果?
  2. scaffolding 过程中, 是否有些中间文件或提示, 让我知道挂载问题是由于 chimeric contigs 或 collapsed contigs, 以便针对性调整?
  3. 对于可公开数据, 如您提到的马铃薯, 甘蔗, 是否方便 share 一些使用的参数? 我想尽量学习一下您是如何利用 HapHiC 解决问题的.

再次感谢您的回复与帮助!

zengxiaofei commented 3 months ago

chromap

上图是chromap和bwa mem的直观对比，一个较高杂合的同源四倍体的案例，juicer pipeline的情况比chromap还要严重很多。

同源染色体之间的互作其实本来和普通的染色体之间的信号没有什么区别，我们从热图上能观察到这种沿对角线分布的信号主要有两个来源，比对和组装错误。当使用bwa mem比对配合MAPQ>=1的过滤标准能很好地保证比对地准确性。组装错误可能来源于纠错过程或者read的分配过程，最终这些来源于不同haplotype的reads被collapse生成了consensus序列，造成了碱基水平的switch error。HapHiC --remove_allelic_links <ploidy>这个参数就是根据这个对角线信号识别等位contig，并消除这些信号带来的不好影响。所以如果使用chromap做单倍体分型，建议加上这个参数。

我也经常遇到由于本身数据问题造成的挂载困难; 针对您提出识别组装错误的问题, 我们是否有办法在获得组装结果时就做出识别, 或调整组装结果

目前来说没有很好的办法，有些组装问题可能需要更好的数据才能解决，比如使用ONT ultra-long reads。用gfa文件绘制assembly graph可以一定程度上识别collapse，但比较难以统计，另外gfa文件中的rd tag可以看到contig的深度，高深度通常也是collapse导致的。

另外，在有限的数据的情况下，选对组装版本是很重要的。比如二倍体我们倾向于使用hifiasm的Hi-C模式组装，当杂合率足够的时候选择将hap1.p_ctg和hap2.p_ctg 合并，将Hi-C reads同时比对到合并的contigs上进行scaffolding。当杂合率很低时，分开对hap1和hap2 scaffolding，最后再合并比对进行验证。如果时多倍体，我们会选择p_utg进行scaffolding，因为unitigs是准确分型的。

scaffolding 过程中, 是否有些中间文件或提示, 让我知道挂载问题是由于 chimeric contigs 或 collapsed contigs, 以便针对性调整?

在01.cluster/HapHiC_cluster.log文件中，我们会输出一些可能有助于发现这些问题的信息。里面会提示利用link density和rank-sum值分别过滤了多少fragments，如果这个比例比较高，比如超过10%，说明基因组可能存在较多的组装错误。当然这也不是绝对的，但是算作一种hint。因为我在自动过滤这些fragments时使用的策略是最常见的离群值检测方法，即Q3+1.5*IQR，如果错误占比太大，这个方法会失效，所以提供了手动设置比例的选项。这些错误通常可以安全地进行过滤，因为这些过滤只是暂时的，后面contigs还是会补回到合适的group。另外如果想知道具体是哪些contig被过滤了，可以在运行haphic pipeline或haphic cluster时加上--verbose参数查看完整的日志。

对于可公开数据, 如您提到的马铃薯, 甘蔗, 是否方便 share 一些使用的参数? 我想尽量学习一下您是如何利用 HapHiC 解决问题的.

我在预印文章的附件中详细提供了所有用到的命令行，可以参考。另外最近文章也快正式发表了，到时可以再看新版的附件，做了一些补充和改动。比如针对发表的C88.v1版本，也做了测试。

感觉像是基于 ALLHiC 自己改了的软件

这里我也想讲一讲HapHiC和ALLHiC的关系。其实我本来只是ALLHiC的用户之一，但是在做自己的项目时遇到了一些问题，做过一些简单的外围开发（ https://github.com/zengxiaofei/small-contig-reassignment ）。在发现效果不错后，我决定自己从头设计一款工具，尝试解决ALLHiC依赖参考基因组的问题。因为总体上和ALLHiC的策略是一致的，都是先聚类后排序（相似的还有LACHESIS，有别于3D-DNA、SALSA2和YaHS），所以我对ALLHiC的文件格式做了兼容。实际上如果用得熟练的话，这两款工具可以混合使用。另外在排序的最后一步，我调用了一个自己修改过的allhic程序（ https://github.com/zengxiaofei/allhic ），从fast sorting的初步排序结果进行优化。当然，用--skip_allhic参数跳过allhic也是可以的。事实上fast sorting的结果其实非常稳健。最近的更新也是基于fast sorting和allhic的结果进行比较，从中择优。总的来说，这款工具和ALLHiC相比，区别还是很大的，无论是contig的纠错、过滤、等位contig的识别、聚类方法、reassignment的思想，排序都做了比较多的创新。

最后感谢您的教程，感到有人关注和使用非常荣幸: )