Canu的graph和contig生成过程

这篇文章是我发现canu输出的contig可能存在misassembly(错误组装), 为了探索这种错误是如何产生的，我尝试解决如下问题

• Canu输出contig的基本步骤
• contig和GFA是什么关系
• 如何提取contig对应的read
• 如何检查contig的graph

contig的构建过程

Canu的流程分为三个步骤，前两步是原始输入的纠错，最后一步是基于纠错后的reads来构建contigs。

步骤三的流程图

步骤三的输出文件

0-mercounts 1-overlapper 3-overlapErrorAdjustment 4-unitigger 5-consensus prefix.ctgStore prefix.ctgStore.coverageStat.log prefix.ctgStore.coverageStat.stats prefix.ovlStore prefix.ovlStore.config prefix.ovlStore.config.txt prefix.ovlStore.per-read.log prefix.ovlStore.summary prefix.utgStore

注: prefix表示文件名前缀。

根据流程和输出文件可知这一步骤可以细分为5步

• 构建read之间的overlap得到ovlStore
• 分析overlap检测错误
• 重新计算overlap的alignment
• 使用bogart构建contig
• 对contig进行consensus
• 构建输出, graph, layout, sequences

这里面最核心的是bogart, 它根据overlap情况构建graph，4-unitigger有如下的GFA文件，就是运行时产生的文件。

prefix.best.edges.gfa prefix.initial.assembly.gfa prefix.final.assembly.gfa prefix.contigs.gfa prefix.unitigs.gfa prefix.unitigs.aligned.gfa prefix.contigs.aligned.gfa

bogart基于graph寻找最优路径, 最终得到unitig和contig两个结果，其中unitig更加碎一些，但更加准一些(Contigs, split at alternate paths in the graph)。

最终的contig还需要进一步的consensus，得到最终的输出

• prefix.contigs/unitigs.fasta: Fasta序列
• prefix.contigs/unitigs.gfa: contig.gfa已经在Canu1.9之后被移除，并非记录Fasta的构建过程
• prefix.contigs/unitigs.layout: 上述GFA并不存放序列，实际序列在layout中
• prefix.contigs/unitigs.layout.readToTig: contig ID和read ID的对应关系
• prefix.contigs/unitigs.layout.tigInfo: 每个contig的信息
• prefix.unitigs.bed: untig和contig的对应关系

Contig和GFA的关系

对于主目录输出结果中contig/untigs对应的fasta和gfa，我们需要明白其中fasta并非是由gfa生成，而是记录了contig与之前contig的overlap关系。

因此，如果用bandage可视化Canu 1.8之前的GFA文件，实际上你看到的是最终的contig之间的关系，而非read之间的关系。在Canu 1.9的更改日志中写道，

• Output file ‘contigs.gfa’ was removed because it was misleading

我们实际需要的是canu通过解析read之间的overlap关系图来得到最终的contig的gfa文件。而Canu并没有提供符合我们需求的文件，因为Canu输出的GFA文件中都没有以P为开头的记录，即没有记录contig的路径。

假如我对其中一个contig存在怀疑，那么我应该如何检查该contig对应read的overlap关系呢？

提取contig对应的read信息

注意，不同版本的Canu输出数据库未必兼容，也就是用和建立数据库不同版本的工具会报错。

以我自己的一个数据为例，对于.gfa中的一个编号，例如tig00000007，

S tig00000007 * LN:i:62235

我们可知他的实际ID是7。它对应的read信息可以在.layout.readToTig找到。

$ awk '$2==7' prefix.contigs.layout.readToTig | head -n 2  7 ungapped 0 37356 56109 7 ungapped 507 38683

第一列即是contig对应的readID，例如

利用ovStoreDump我们就可以提取和此read overlap的所有read

ovStoreDump -S prefix.seqStore -O unitigging/prefix.ovlStore -picture 

ovStoreDump可以以多种形式输出overlap的信息，例如-picture就是以ASCII展示overlap

或者根据readID用sqStoreDumpFASTQ提取实际的序列

sqStoreDumpFASTQ -S prefix.seqStore -o - -fasta -r  -raw | seqkit seq - | head sqStoreDumpFASTQ -S prefix.seqStore -o - -fasta -r  -corrected | seqkit seq - | head sqStoreDumpFASTQ -S prefix.seqStore -o - -fasta -r  -trimmed | seqkit seq - | head

-raw/-corrected/-trimmed表示提取不同阶段的结果

输出结果里使用的实际输入序列的ID编号，可以在prefix.seqStore/readNames.txt中找到readID和原始编号的对应关系。

掌握上述的操作，就可以提取指定的contig的read，然后将read回帖到该contig上，利用IGV可视化。

可视化检查contig的GFA

提取contig对应的readID, 然后根据readID抽取gfa

contigID=$1 prefix=$2 readtotig=$prefix.contigs.layout.readToTig gfa=unitigging/4-unitigger/$prefix.initial.assembly.gfa awk -v id=$contigID '$2==id {printf("read%08d\n",$1)}' $readtotig > readid.txt grep -f readid.txt $gfa > $contigID.gfa

输出的.gfa就可以用Bandage进行可视化。但实际上发现这个操作并不太可行，不如直接上一节用IGV可视化alignment来的直接。