Canu的graph和contig生成过程
这篇文章是我发现canu输出的contig可能存在misassembly(错误组装), 为了探索这种错误是如何产生的,我尝试解决如下问题
- Canu输出contig的基本步骤
- contig和GFA是什么关系
- 如何提取contig对应的read
- 如何检查contig的graph
contig的构建过程
Canu的流程分为三个步骤,前两步是原始输入的纠错,最后一步是基于纠错后的reads来构建contigs。
步骤三的流程图
![](/assets/blank.gif)
步骤三的输出文件
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
632820 7 ungapped 0 37356
56109 7 ungapped 507 38683
第一列即是contig对应的readID,例如632820
利用ovStoreDump
我们就可以提取和此read overlap的所有read
ovStoreDump -S prefix.seqStore -O unitigging/prefix.ovlStore -picture 632820
ovStoreDump可以以多种形式输出overlap的信息,例如-picture就是以ASCII展示overlap
或者根据readID用sqStoreDumpFASTQ
提取实际的序列
sqStoreDumpFASTQ -S prefix.seqStore -o - -fasta -r 632820 -raw | seqkit seq - | head
sqStoreDumpFASTQ -S prefix.seqStore -o - -fasta -r 632820 -corrected | seqkit seq - | head
sqStoreDumpFASTQ -S prefix.seqStore -o - -fasta -r 632820 -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来的直接。
![](/assets/blank.gif)
Canu的graph和contig生成过程相关推荐
- PacBio-组装介绍
主页:github: PacificBiosciences/FALCON 转自:https://www.cnblogs.com/leezx/p/5724590.html 简介 Falcon是一组通 ...
- Falcon:三代reads比对组装工具箱
主页:github: PacificBiosciences/FALCON 简介 Falcon是一组通过快速比对长reads,从而来consensus和组装的工具. Falcon工具包是一组简单的代码集 ...
- HG-CoLoR用一个变阶de Bruijn graph混合校正高噪声长读数
HG-CoLoR用一个变阶de Bruijn graph混合校正高噪声长读数 1 简介 几年来,长读长测序技术不断发展,可以解决大型复杂基因组的组装问题,在此之前,仅使用短读长测序技术很难解决这些问题 ...
- 图神经网络(Graph Neural Networks,GNN)综述
鼠年大吉 HAPPY 2020'S NEW YEAR 作者:苏一 https://zhuanlan.zhihu.com/p/75307407 本文仅供学术交流.如有侵权,可联系删除. 本篇文章是对论文 ...
- 三代组装软件canu学习笔记
三代组装软件canu学习笔记 (2017-08-07 14:17:43) 转载▼ 分类: 三代 1:这个组装软件起源于PBcR包含在Celera Assembler中(http://wgs-ass ...
- 使用Canu对三代测序进行基因组组装
Canu简介 Canu是Celera的继任者,能用于组装PacBio和Nanopore两家公司得到的测序结果. Canu分为三个步骤,纠错,修整和组装,每一步都差不多是如下几个步骤: 加载read到r ...
- 如何理解 Graph Convolutional Network (GCN)?
几年前如果熟练使用TensorFlow,同时掌握基本的AI算法就可以很容易找到一份高薪的工作,但现在不一样了,AI岗位的要求越来越高,对知识的深度也提出了更高的要求. 如果现在一个面试官让你从零推导S ...
- 论文阅读课4-Long-tail Relation Extraction via Knowledge Graph Embeddings(GCN,关系抽取,2019,远程监督,少样本不平衡,2注意
文章目录 abstract 1.introduction 2.相关工作 2.1 关系提取 2.2 KG embedding 2.3 GCNN 3. 方法 3.1符号 3.2框架 3.2.1 Insta ...
- 如何理解 Graph Convolutional Network(GCN)?
几年前如果熟练使用TensorFlow,同时掌握基本的AI算法就可以很容易找到一份高薪的工作,但现在不一样了,AI岗位的要求越来越高,对知识的深度也提出了更高的要求. 如果现在一个面试官让你从零推导S ...
最新文章
- Lesson2 Hello,GLSL
- Python进阶之路:namedtuple
- Django之orm补充
- Hibernate事实:有利于双向集vs列表
- 数据库实验四 用户权限管理
- 观视屏《残疾人郑心意》所想
- 直播带货还有机会吗?
- C# 获取文件名及扩展名【转】
- SPSS统计分析与行业应用案例详解
- solidworks迈迪插件_迈迪工具集V55特别PJ版_打包下载
- 对冲之王 - 华尔街量化投资传奇 读后感
- PCB板-叠层详细介绍
- 含protobuf程序运行时与libqgtk3.0.so冲突
- Python语法易混淆
- java共享文件夹SMB1服务报错jcifs.smb.SmbException: Failed to connect: 0.0.0.0<00>/122.168.23.26
- 安利这几个网站和软件给你
- Android篇 --Notification(消息通知)
- 3dmax入门学习基础教程第1部分:建模
- [SECURITY]--01
- 雄迈NVR、DVR设置开启LOGO
热门文章
- python量化策略——Fama-French三因子模型
- 172.16.82.0/25的含义,IP段,掩码
- 大众PLC自动控制系统,VASS标准项目应该知道的几个知识点
- 使用freeSWITCH和Yate进行VoIP通话
- 《B站-Redis教程》学习笔记
- Hu矩的学习,图像轮廓特征识别(二,c#实现)
- AD20/Altium designer——过孔盖油
- 【jQuery】简单的网页文本格式编辑器ckeditor
- 详细讲解ExpandableListView显示和查询仿QQ分组列表用户信息
- AI绘画5大免费工具