主页:github: PacificBiosciences/FALCON

转自:https://www.cnblogs.com/leezx/p/5724590.html

简介


Falcon是一组通过快速比对长reads,从而来consensus和组装的工具。

Falcon工具包是一组简单的代码集合,我使用它们来研究单倍体和二倍体基因组的高效组装算法。

为了提高计算速度,它有一些后台代码是使用C来实现的,为了方便一些简单的前端是用Python编写的。

Falcon不是一个傻瓜的组装工具(除了很小的基因组),为了得到最好的结果,你可能需要了解各种分布式计算系统和一些基本的基因组组装理论。FAQ页面有一些常见的问题及其回答。

  • Manual
  • Tips
  • Contrib

Falcon 基因组组装工具包 - 使用指南

1. 分层次基因组组装过程的概述

主要有如下几个步骤:

  1. Raw sub-reads overlapping for error correction
  2. Pre-assembly and error correction
  3. Overlapping detection of the error corrected reads
  4. Overlap filtering
  5. Constructing graph from overlaps
  6. Constructing contig from graph

简单翻译:

  1. 使用Raw sub-reads 构建重叠,准备进行错误校正
  2. 预组装和错误校正
  3. 错误校正后的reads的重叠检测
  4. 重叠过滤
  5. 用重叠来构造图
  6. 用图来构造contig

每个步骤使用不同的命令行工具,实现不同的算法来完成这项工作。同样,每一步的计算需求也大不一样。假定用户已经拥有合理数量的计算资源。

例如,在合理的时间内装配100M大小的基因组,可能需要至少32核cpu 和 128 Gb RAM。代码是按照集群计算环境来编写的。需要一个作业队列来进行长时间的脚本工作,以及cpu-rich的计算工作队列。

一个包含了一组sub-reads的文件,fc_run.py脚本可以驱动工作流来管理和检查数据依赖和提交每一步的工作,从给定的数据中生成一个组装的草图。

fc_run.py是工作流驱动脚本,需要运行在一台计算机上,允许持续长时间地工作一段时间来完成整个装配过程。它需要一个配置文件fc_local.cfg作为单个输入。原始序列数据的输入文件包含在配置文件中。

配置文件可用于控制计算资源的使用和重要算法参数,根据输入数据集达到最有组装的效果。然而,没有万能的方式来猜测各种具体项目的最合适的选项和参数。参数调优的最好方法是了解一些组装理论,一点点实现方法,这样才能正确的了解更改某些选项对组装结果的影响。快速浏览reads长度分布、整体覆盖和调整相应的选项也是非常重要。

这里我们会详细介绍基因组分层组装策略 以及 fc_run.py 各个参数的含义

快速入门
  • Setup:-Installation-and-Environment
  • Setup:-Running
  • Setup:-Complete-example

2. 使用Raw sub-reads构建重叠(为了错误校正)

这个版本的 Falcon 的 overlapping 是用 Gene Myers' Daligner 的修改版本来做的(参考 github 的 forked version)。主要的修改是adapting some I/O for the downstream bioinformatics process。可以通过git diff来查看有哪些修改。

git diff    #查看版本之间差异

input_fofn 选项指向了包含所有输入数据的文件。来自Dalignerfasta2DB 在 fc_run.py 会被调用(I/O集中器,会在你执行 fc_run.py的电脑节点出运行,如果这在你的服务器集群中是个问题,你必须修改代码来打包相关的操作到一个脚本中,以便提交到你的任务管理系统中)。

input_fofn    #指向了包含所有输入数据的文件

这个版本的 fc_run.py 支持从错误校正的reads中运行组装程序(可以跳过错误校正环节,直接开始组装)

input_type = preads  #输入为 all error-corrected reads,会跳过错误校正,直接开始最终的重叠群组装步骤
input_type = raw    #原始数据

你需要决定cutoff的长度,一般,选择一定的阈值,使得你能得到15x 到 20x的覆盖度是比较好的。 然而,如果计算资源充足,校正reads你有其他用途,你可以设置lower length cutoff,来获得更多的error corrected reads。 
最终组装的时候,更多的reads不意味着更好的组装结果,因为会有噪音。 
有时,尝试不同的length_cutoff_pr是有意义的,较初次overlapping step for error correction而言,后期计算更cheap。 
我们的策略是选择一个短的length_cutoff 进行一次计算。之后, 我们使用不同的length_cutoff_pr来获得更好的组装结果。

length_cutoff    #控制错误校正过程中的cutoff
length_cutoff_pr    #控制之后组装重叠群步骤的cutoff

pa_concurrent_jobs选项控制由fc_run.py 提交的 number of concurrent jobs。sge_option_da 和 sge_option_la控制作业队列和daligner jobs 的 number of slots。daligner 使用的线程数默认是4。然而,依据集群设置、计算节点的内存数量,你可以使用超过4的slots。为了选择最好的数量,你最好咨询你们实验室的HPC专家,做一些前期的测试。

pa_concurrent_jobs    #fc_run.py提交的当前作业数
sge_option_da    #作业队列
sge_option_la    #daligner作业的slots数

最终运行作业的数量取决于你怎么"splits"序列数据库,仔细阅读Gene Myers's blog ( http://dazzlerblog.wordpress.com ) ,了解如何设置 pa_DBsplit_option 和 pa_HPCdaligner_option选项。一般,对大型基因组而言,在pa_DBsplit_option中,你应该使用-s400(400Mb sequence per block)。这样作业的数量会变少,但是持续时间会延长。如果,作业调度系统限制了一个作业可以运行的时间,你应该设置a smaller number for the -s option。 
另一个影响作业总数的参数是pa_HPCdaligner_option中的 -dal option。它决定how many blocks are compared to each in single jobs。大的 -dal 给出了大的作业,但是数量少;反之。小的,作业小,但是提交的数量多。 
这个工作流中,daligner生成的trace point没有被用到(为了高效,应该使用trace point,但是首先你的知道如何正确的pull out它们)。pa_HPCdaligner_option 中的 -s1000 使得trace point变稀疏,以此节约磁盘空间。我们使用-l1000忽视短于1kb的reads。

pa_DBsplit_option    #-s400  for large genomes    smaller number for the -s
pa_HPCdaligner_option    #  -dal    #  determines how many blocks are compared to each in single jobs;    -s1000   #-l1000   #

3. 预组装和错误校正

daligner的输出时一组.las文件,它包含了reads间相互比对的信息。Such information is dumped as sequences for error correction by a binary executable LA4Falcon to fc_consensus.py.(什么意思?) fc_consensus.py生成了consensus(这部分是用C写的)。

fc_consensus.py有许多选项,可以使用falcon_sense_option来控制它。大部分情况下,--min_cov 和 --max_n_read是最重要的选项。--min_cov控制了什么时候a seed read gets trimmed or broken due to low coverage(什么意思?),--max_n_read puts a cap on the number of reads used for error correction(什么意思?)。 
在高度重复的基因组汇总,设置小的 --max_n_read,来确保 consensus code does not waste time aligning repeats。最长的合适的overlaps 被用来进行correction ,来减少the probability of collapsed repeats。

你可以使用cns_concurrent_jobs来设置提交给任务管理系统的当前作业的最大数量。

falcon_sense_option    #--min_cov    #when a seed read gets trimmed or broken due to low coverage???--max_n_read    #
cns_concurrent_jobs    #

4. 错误校正后的reads的overlaps检测

这部分与最初的overlapping 十分相似,只有部分修改,daligner只以本地的raw reads作为默认值。 
fc_run.py generates a fasta file of error-corrected reads where the fasta header is parse-able by daligner.(什么意思?) 
下面的参数控制这个步骤的计算过程:

sge_option_pda = -pe smp 8 -q jobqueue
sge_option_pla = -pe smp 2 -q jobqueue
ovlp_concurrent_jobs = 32
ovlp_DBsplit_option = -x500 -s50
ovlp_HPCdaligner_option = -v -dal4 -t32 -h60 -e.96 -l500 -s1000

这个设置与第一步的overlapping 相平行,主要的区别是ovlp_HPCdaligner_option-e 选项。错误率小了很多,因此我们期待 much higher correlation between the p-reads。

5. overlaps过滤

不是所有的overlaps 都是独立的,所以利用某些filtering step来减少计算和组装图的复杂性是可能的。例如:如果一个read被全部包含在另一个read里,那么它们的overlap不会提供额外的构建基因组的信息。由于重叠关系的传递属性(transitive property),许多overlap 信息可以被简单的推测出来。 
事实上,构建contigs 的第一步就是删除 transitive reducible edges(???),这意味着,我们只需要best n overlaps 在5' or 3' 端。overlap_filtering_setting option的--bestn parameter参数可以被用来控制每个read报告的maximum overlap。

overlap_filtering_setting --bestn  #control the maximum overlap reported for each read

另一个有用的启发式,仅保留有平均 5' and 3' 覆盖的reads。这是因为,如果一个read以重复结束,它会有更高的重复。这种reads不提供价值,不能解决相关的重复问题。我们可以过滤掉它们,希望能有reads(跨越重复,在两端有正常的独特的anchors)。同样,如果一个reads一端的覆盖度太低,那它可能有太多的测序错误。这种 reads产生了 "spurs"在组装图中,通常都会将它们过滤掉。--max_cov和 --min_cov就是用来过滤那些有太高或者太低的reads。

--max_cov   #过滤太高的reads
--min_cov   #过滤太低的reads

过滤的脚本也允许过滤一些 "split" reads。如果一个read在5' 和 3'端有非常不相等的覆盖度,这也是某一段是重复的信号。 --max_diff可以用来过滤掉一端比另一端的覆盖哦高很多的reads。

--max_diff   #过滤掉一端比另一端的覆盖哦高很多的reads

使用这些参数的正确的数量是多少?真的很难设置正确,如果 error corrected reads的总覆盖度比知道的length cut off 并且合理的高 (e.g. greater than 20x),那么将平均的覆盖度min_cov设置为5,max_cov设置为3倍将会很安全, max_diff设置为2倍。然而,在低覆盖度的情况下,将min_cov设为1或2将会更好。一个名为fc_ovlp_stats.py的帮助脚本可以帮助dump(丢弃)。 A helper script called fc_ovlp_stats.py can help to dump the number of the 3' and 5' overlap of a given length cutoff, you can plot the distribution of the number of overlaps to make a better decision.(???) 
你也可以设置max_diffmin_cov为很高的值,来避免过滤,在某些特殊的场合。

这个过滤的过程会显著过滤掉高度拷贝重复的信息。也就是,这些重复会被过滤掉,并且不会在最终的组装结果中显示。如果你对这些重复感兴趣,总之,单独处理这些重复会更加有效和实用。如何解决重复序列这个问题具有非常重要的意义。

6. 用overlaps构造图

鉴于overlapping 数据,string graph 由 fc_ovlp_to_graph.py来创建,使用默认的参数。 fc_ovlp_to_graph.py生成一些文件代表组装后的最终的 string graph。最终的ctg_path包含了每一个contig的图的信息。一个 contig是一个simple paths 和 compound paths的线性的path。Compound paths are those subgraph that is not simple but have unique inlet and outlet after graph edge reduction. 
They can be induced by genome polymorphism or sequence errors. By explicitly encoding such information in the graph output, we can examine the sequences again to classify them later.

7. 用图来构造contig

create draft contigs的最后一步就是find a single path of each contig graph,以及generate sequences accordingly。在a contig graph is not a simple path的情况下,我们寻找 end-to-end path that has the most overlapped bases。这被称为 primary contig。图中的每一个path,如果它与 primary one不同,我们会construct the associated contig。在the associated contigs are induced by sequencing error的情况下,the identity of the alternative contig and the primary contig will be high ( > 99% identity most of time)。在there are true structure polymorphism的情况下,there are typically bigger difference between the associated contigs and the primary contigs。

fc_graph_to_contig.py脚本,takes the sequence data and graph output to construct contigs。它会生成所有的 associated contigs at this moment。将来,一些后处理程序,来删除errors中的some of the associated contigs将会被开发出来。(正在开发)

8. 通用的daligner选项

daligner由pa_HPCdaligner_option 和ovlp_HPCdaligner_option控制。 
为了限制内存,我们可以使用-M选项。对于人类基因组组装,我们测试 -M 32,每个daligner使用32G RAM。其他可能性正在调查中。 
更多的daligner的细节,参见: dazzlerblog。


工作目录结构、工作恢复和故障排除

代码设计工作在单一目录。工作目录的典型布局是这样的:

$ ls -l
total 56
drwxr-xr-x 84 jchin Domain Users  8192 Nov 30 12:30 0-rawreads
drwxr-xr-x 18 jchin Domain Users  4096 Nov 30 12:33 1-preads_ovl
drwxr-xr-x  2 jchin Domain Users  4096 Nov 30 12:44 2-asm-falcon
-rwxr-xr-x  1 jchin Domain Users  1041 Nov 30 12:13 fc_run.cfg
-rw-r--r--  1 jchin Domain Users   378 Nov 29 23:20 input.fofn
drwxr-xr-x  2 jchin Domain Users  4096 Nov 30 12:13 scripts
drwxr-xr-x  2 jchin Domain Users 24576 Nov 30 12:33 sge_log

典型的input.fofn是这样的:

/mnt/data/Analysis_Results/m140913_050931_42139_c100713652400000001823152404301535_s1_p0.1.subreads.fasta
/mnt/data/Analysis_Results/m140913_050931_42139_c100713652400000001823152404301535_s1_p0.2.subreads.fasta
/mnt/data/Analysis_Results/m140913_050931_42139_c100713652400000001823152404301535_s1_p0.3.subreads.fasta

1. 内部0-rawreads目录

目录0-rawreads包括:overlapping原始序列的所有的脚本和数据,它包含不同的job_ *m_ *目录。 

例如,如果我们将E. coli的数据分成20个块,目录就会如下:

cns_done   job_00011  job_00024  job_00037  job_00050  m_00003  m_00016
da_done    job_00012  job_00025  job_00038  job_00051  m_00004  m_00017
job_00000  job_00013  job_00026  job_00039  job_00052  m_00005  m_00018
job_00001  job_00014  job_00027  job_00040  job_00053  m_00006  m_00019
job_00002  job_00015  job_00028  job_00041  job_00054  m_00007  m_00020
job_00003  job_00016  job_00029  job_00042  job_00055  m_00008  preads
job_00004  job_00017  job_00030  job_00043  job_00056  m_00009  prepare_db.sh
job_00005  job_00018  job_00031  job_00044  job_00057  m_00010  raw_reads.db
job_00006  job_00019  job_00032  job_00045  job_00058  m_00011  rdb_build_done
job_00007  job_00020  job_00033  job_00046  job_00059  m_00012  run_jobs.sh
job_00008  job_00021  job_00034  job_00047  las_files  m_00013
job_00009  job_00022  job_00035  job_00048  m_00001    m_00014
job_00010  job_00023  job_00036  job_00049  m_00002    m_00015

job_ *目录存储了每个daligner的输出工作。

m_ *目录是用于合并工作。

有一些空文件,它们的名字是done结尾。这些文件的时间戳是用来跟踪工作流阶段。您可以修改时间戳和重启fc_run.py来触发为工作流程的一部分做某些做计算。然而,不推荐这样,除非你有一些时间来读fc_run.py的源代码,并且知道如何跟踪工作流内部的依赖性。(例如,如果你在成功组装后touch rdb_build_done,并且重新运行fc_run.py。因为所有中间过程取决于文件,rdb_build_done比任何中间文件的都要新,就会触发fc_run.py,重复整个装配过程。 

las_files存储了最后比对信息。

如果你不打算重新运行工作,但是想知道如何比对的样子,你可以删除所有job_ *m_ *目录但保留las_filespreads 目录。 

这一步的主要输出存储在0-rawreads /preads0-rawreads/preads内部的out.%04d.fa的是输出的reads的fasta文件。如果这一步成功完成, 哨兵文件cns_done将被创建。

2. 1-preads_ovl目录

该目录存储着p-read与p-read overlapping的数据。它与0-rawreads大体一致,但是却没有consensus 这一步。主要输出是1-preads_ovl/目录里的 *.las

3. 2-asm-falcon目录

这是最终的输出目录。 
它包含了组装图和draft contigs的所有信息,细节描述参照Graph output format一节。


fc_ovlp_to_graph.py 选项 和 组装图输出格式

下面是fc_ovlp_to_graph.py的使用信息:

usage: fc_ovlp_to_graph.py [-h] [--min_len MIN_LEN] [--min_idt MIN_IDT][--lfc]overlap_filea example string graph assembler that is desinged for handling diploid genomespositional arguments:overlap_file       a file that contains the overlap information.optional arguments:-h, --help         show this help message and exit--min_len MIN_LEN  minimum length of the reads to be considered forassembling--min_idt MIN_IDT  minimum alignment identity of the reads to be consideredfor assembling--lfc              use local flow constraint method rather than best overlapmethod to resolve knots in string graph

De-dup example and strategy 
TODO 
Low coverage assembly 
TODO 
Assembly layout rule

CPU使用

Applications of the assembly graph 
TODO 
Reproducibility and replicability

序列 alignment 和 consensus 的C代码

Other TODOs 
Incremental overlapping 
Pre-processing repeat for overlapping

术语

subread:子read, 
full-pass subread: 
pre-assembly: 
error correction: 
proper overlap: 
string graph: 
contig: 
primary contig: 
associated contig: 
p-read: 
compound path: 
simple path: 
Quiver:

科普时间: 
什么是virtualenv? 
每个应用可能需要各自拥有一套“独立”的Python运行环境。virtualenv就是用来为一个应用创建一套“隔离”的Python运行环境。 
资料: 
virtualenv - 廖雪峰 
virtualenv -- python虚拟沙盒 
Virtualenv入门基础教程

帮助文档

现在默认分支是“master”,其包含Falcon最新的源码。 
当前最新的综合版本是v0.3.0,安装指南参照:v0.3+ Integration Installation Guide

上一个版本是v0.2.2,具体参照:v0.2.2 release github repository

  • wiki pages
  • 开发者安装指南
  • v0.3+综合安装指南
  • 使用和帮助文档
  • FAQs
  • v0.2.2 release github repository
  • Falcon发行说明

PacBio-组装介绍相关推荐

  1. RNA-seq技术之转录组从头组装介绍

    RNA-seq技术之转录组从头组装介绍 转载2016-05-31 17:07:39 1.何为转录组组装 说起转录组组装,不得不先说新一代测序技术(next generation sequencing) ...

  2. 「三代组装」Pacbio组装后如何用自身数据进行polish(更新版)

    之前那我由于需要对PacBio的组装结果进行polish,于是写了「三代组装」Pacbio组装后如何用自身数据进行polish.最近发现自己又有了需求,于是重新回顾了我之前写的这篇文章,但是在实践的时 ...

  3. 「三代组装」Pacbio组装后如何用自身数据进行polish

    三代数据由于其高错误率(目前应该是10%左右), 即便在组装前有一步纠错环节,但是组装得到序列依旧存在着许多错误,因此需要进行polish环节.polish分为两个层次,三代原始序列polish和二代 ...

  4. 常用转录组组装软件集合

    转录组组装软件 基因组组装 基因组组装(Genome assembly)是指使用测序方法将待测物种的基因组生成序列片段(即read),并根据reads 之间的重叠区域对片段进行拼接,先拼接成较长的连续 ...

  5. SPAdes混合组装二代、三代测序数据

    导读 SPAdes是2012年发表在Journal of Computational Biology上的一篇文章提出的二代测序组装软件,是目前引用量已经达到6200+,在宏基因组组装软件中引用量最高[ ...

  6. 计算机基础与组装,1.计算机基础与组装.pptx

    1.计算机基础与组装 Computer foundation计算机基础<计算机基础>课程结构计算机基础与组装网络数字与系统基本操作Word 2007 基本操作计算机基础Excel 2007 ...

  7. 壮观霉素抗性基因原理_基因组学深入挖掘·研究方案(下篇)

    前情回顾 上次小编为大家讲解了四种以基因组为基础的多组学联合研究方案(基因组与转录组,深入挖掘基因表达信息:基因组联合代谢组与转录组,锁定关键通路:基因组与群体进化,解析物种发展历程:基因组结合GWA ...

  8. 科研快报 | 三代测序技术-海水微生物态,助力海水微生态及微生物基因组研究

    PacBio研究专题 二代测序读长偏短,环境宏基因组样品研究受到了很大限制.作者通过三代测序对来自地中海的冬季混合海水样本进行宏基因组测序.利用PacBio Sequel II平台的超长读长明显可以提 ...

  9. unity应用开发实战案例_Unity3D游戏引擎开发实战从入门到精通

    Unity3D游戏引擎开发实战从入门到精通(坦克大战项目实战.NGUI开发.GameObject) 一.Unity3D游戏引擎开发实战从入门到精通是怎么样的一门课程(介绍) 1.1.Unity3D游戏 ...

  10. 一种PacBio测序数据组装得到的基因组序列的纠错方法技术 (专利技术)

    一种PacBio测序数据组装得到的基因组序列的纠错方法技术 技术编号:17008244阅读:83留言:0更新日期:2018-01-11 04:20 本发明专利技术提供一种PacBio测序数据组装后序列 ...

最新文章

  1. 成功解决利用pandas输出DataFrame格式数据表时没有最左边的索引编号(我去,这个问题折腾了我半个多小时)
  2. C# 读写excel 用于导入数据库 批量导入导出excel
  3. 上项线体表位置_实用人体体表解剖:头颈部(高清大图版)
  4. SpringBoot优点
  5. 计算机扩招的学校,这些985/211学校今年继续扩招!又有一大波学校专业课有变!...
  6. tensorflow2.0学习(一)
  7. 红米手机5怎么样卡刷开发版开启root超级权限
  8. ****CentOS下安装JDK1.7
  9. 论文解读——An Analysis of Scale Invariance in Object Detection – SNIP
  10. Tomcat介绍,安装jdk,安装tomcat,配置Tomcat监听80端口
  11. 拓端tecdat|R语言社区主题检测算法应用案例
  12. 基于Jackson2的JsonSchema实现java实体类生成json(一)
  13. 微擎支持html微信支付,微信小程序云开发:现已原生支持微信支付
  14. 爬取豆瓣电影排行榜(评分)
  15. 超详细的阿里云服务器购买及远程连接开机(Win系统)
  16. 电子科技大学计算机科学与技术考研复试,电子科技大学计算机科学与工程学院2021考研招生复试工作安排...
  17. 哪个心情不好来看看,老逗了
  18. springBoot项目改名
  19. 剖析RS-485原理以及与其他总线的区别
  20. Java 提取PDF图片(pdfbox)Extract PDF document images

热门文章

  1. F - Construct Highway
  2. 云原生CI/CD框架Tekton国内部署方式
  3. Lieb格子上SU(3)Hubbard模型铁磁的基态
  4. 基于迁移学习的语义分割算法分享与代码复现
  5. 空间滤波-高斯低通滤波器
  6. “智慧警务”:信息科技下的平安龙岗
  7. java final类为什么不能继承_浅谈Java之终止继承:Final类和Fianl方法
  8. 基于哈工大生物信息学专业培养方案提炼的学程
  9. 硬盘2.5寸4tb服务器硬盘,西部数据My Passport 2.5英寸4TB移动硬盘
  10. 比起“呵呵”,我宁愿听到一句:我现在不想聊天,滚蛋。