关注我们并在后台回复 “进群”，即可加入农心生信工作室学习交流群，群内不定时分享源代码及示例文件，并在线交流答疑。我们在农心生信，等你！
由于微信改版，乱序推送，很多朋友反映收不到公众号推文。快跟着图片步骤，将农心生信公众号设为星标，不错过每一条精彩内容！！
本公众号已开通留言功能，欢迎大家直接在文末留言交流！
写在前面
在完成基因组组装后，我们需要对组装好的基因组进行评估，来判定基因组的组装质量如何。一般来说，需要从基因组的完整性、准确性、连续性等方面评估。目前，大家比较常用的评估方式，包括BUSCO、Merqury、三代或二代数据回贴，植物基因组上也会用到LAI值评估基因组质量。最近，拜读了CJ老师团队的新作：The pineapple reference genome: Telomere‐to‐telomere assembly, manually curated annotation, and comparative analysis，发现作者运用了一个名为 Clipping reveals assembly quality (CRAQ) 的软件来评估基因组组装质量。或是我孤陋寡闻了，之前确实没听过这个软件，按参考文献检索了一下，发现这是去年就发表在Nature Communications上的软件，文章题目叫Identification of errors in draft genome assemblies at single-nucleotide resolution for quality assessment and improvement。
软件原理
浏览了一下原文，这个软件开发的背景是基于基因组组装过程中会出现的两类错误：小尺度错误（小于50bp）和大尺度错误（结构崩坏或扩增，大于50bp）。小尺度错误，如局部插入或缺失，会影响基因组的准确性，但它们通常位于重复区域附近，对下游scaffold的构建影响相对较小。相比之下，大规模的结构错误（例如，由于两个不相关的片段之间的不当连接而产生的错误连接的contig）可能导致形成错误的scaffold，并在多个scaffolds之间传播错误；这可能对下游的进化或比较基因组学研究产生很大影响。CRAQ通过将原始测序reads映射回基因组组装草图，利用裁剪信息来揭示组装错误和低质量区域。这使得能够在单核苷酸分辨率下识别组装错误、杂合位点和单倍型之间的结构差异。通过整合NGS（Next generation sequencing）和SMS（long-read single molecule sequencing）映射，CRAQ可以在不同尺度上识别组装错误，并将错误计数转换为相应的组装质量指标（AQIs），这些指标反映了区域和结构层面的组装质量。
要想大概了解这个软件的流程，先明确下面几个概念：
• NGS和SMS：就是二代测序reads和长读长的三代测序reads，SMS包括pacbio CLR、ONT、pacbio HiFi，这是CRAQ的直接输入的数据
• clipping information：这是CRAQ软件的核心概念，软件就是用来检测clipping information的。clipping information翻译过来是裁剪信息，但其实这个翻译反而不好理解。原文中其实对这个概念有非常清晰的解释：理想情况下，当原始reads映射回基因组时，高质量的基因组组装应展现出均匀的原始读数覆盖度，并且间隙区域或单核苷酸多态性（SNP）簇较少。但实际上，出现一些小的低映射深度或连续的碱基对错配迹象是很常见的，这种类似于SNP簇特征。对于存在大型结构组装错误的区域，例如两个基因组片段的错误连接，映射的读数通常显示出“clipped reads”的特征，这是一种只有read的一部分与参考序列对齐的现象。clipping information就是这一类区域所包含的信息。
• CREs和CSEs：CRAQ根据reads mapping的覆盖率和是否有被裁剪的reads，将假定的错误分为基于裁剪的区域错误(CREs)或基于裁剪的结构错误(CSEs)。如果一个具有裁剪的NGS reads的区域只被具有SNP簇特征的SMS长reads跨越，则将其指定为CRE。如果错误区域周围的映射SMS reads显示出裁剪特征（即，NGS reads同时显示出裁剪或无覆盖率），则将其指定为CSE。CSE的存在意味着基因组组装中存在错误的连接，这可能对组装的可用性产生显著的影响。
• AQI：assembly quality index (AQI)是文中提出的一种组装质量指标。作者通过对大量具有不同质量的基因组组装进行评估，将基因组分类如下：AQI>90，reference quality；AQI在80-90之间，high quality；AQI在60-80之间，draft quality；AQI<60，low quality。小区域和大结构片段的组装质量可以分别计算为R-AQI和S-AQI。
了解上述概念后，我们再看CRAQ的运行流程图，流程分5步：
(I) 原始数据的处理和比对：原始的NGS的短读和SMS的长读分别被映射到基因组组装上。过滤掉低质量的reads后，生成两个比对文件。
(II) 记录特定区域和reads的信息：记录没有NGS reads比对的区域（即间隙）和NGS/SMS reads被clipped的位置。同时记录每个位置上被clipped的reads的数量和总的reads覆盖度。
(III) 定义杂合位点和断点：对于NGS和SMS比对中被clipped的位置，基于用户定义的clipped reads数量与总比对reads数量的比率阈值，定义杂合位点和映射断点。与间隙一起，这些断点被定义为可能的组装错误位置。
(IV) 对可能的错误进行分类：基于reads映射的状态，可能的错误被进一步分类为基于剪切的区域错误（CREs）或基于剪切的结构错误（CSEs）。CREs被定义为那些有NGS断点或被SMS长读跨越的间隙，但富含碱基错配的区域。CSEs被定义为那些SMS剪切断点靠近NGS断点或间隙区域的错误。杂合区域也基于相似的标准，但考虑到映射覆盖率的比率，被分类为基于剪切的区域杂合性（CRH）或基于剪切的结构杂合性（CSH）区域。
(V) 可视化与输出：识别的CREs和CSEs被可视化，并进一步用于基准测试基因组组装的质量。CRAQ为每个组装片段输出全基因组总结、区域AQI分数以及CREs和CSEs的精确位置。
安装使用
亲测后发现，CRAQ的安装和使用非常简单。下载地址如下：
https://github.com/JiaoLaboratory/CRAQ
先安装依赖：
SAMtools(1.3.1+)
Minimap2 (2.17+)
Perl (version >= 5)
pycircos (python 3.7later)
其中pycircos不是必须的依赖，需要画图可以安装。剩下三个必须安装，都可以使用conda安装。
安装CRAQ一行搞定：
git clone https://github.com/JiaoLaboratory/CRAQ.gitCRAQ的核心命令craq在clone下来的CRAQ目录下的bin目录下，可以直接把路径添加到环境变量，也可以直接用绝对路径调用命令。
CRAQ需要一个FASTA（.fa）格式的基因组组装文件和两个代表NGS和SMS测序数据的fastq（.fq）/fasta（.fa）文件作为输入。--map | -x 参数后接'map-hifi'用于PacBio HiFi，'map-pb'用于PacBio CLR，'map-ont'用于ONT。
craq -g assembly.fa -sms SMS.fa.gz -ngs NGS_R1.fa.gz,NGS_R2.fa.gz -x map-hifi另外，也可以提前将reads映射到基因组上，并提供两个已排序和索引的bam文件（sort.bam和sort.bam.bai）。
craq -g assembly.fa -sms SMS_sort.bam -ngs NGS_sort.bam运行之后会输出一个文件夹，里面文件很多，如果是做基因组组装质量评估，只需要看out_final.Report这个文档就行。这个文件记录了整个基因组整体的AQI值和每个contig/scaffold/chr的AQI值。
写在最后
亲身使用后发现，CRAQ安装便利、使用方便，而且了解CRAQ的原理后，也是承认这个软件评估的科学性的。事实上，评估一个基因组的组装质量，从越多的角度和思路上评估，综合起来获得的结论也许就会更准确。CRAQ确实是一个不错的基因组质量评估软件，大家不妨亲自试试，用它来评估一下自己组装的基因组究竟质量如何。
参考文献
Li, K., Xu, P., Wang, J. et al. Identification of errors in draft genome assemblies at single-nucleotide resolution for quality assessment and improvement. Nat Commun 14, 6556 (2023). https://doi.org/10.1038/s41467-023-42336-w
END
编辑 | Narcissus
供稿 | Deeecade
审核 | 农心生信工作室
