Falcon：三代reads比对组装工具箱

主页：github: PacificBiosciences/FALCON

简介

---

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 选项指向了包含所有输入数据的文件。来自Daligner的fasta2DB 在 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_diff和min_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_files和preads 目录。
.
这一步的主要输出存储在0-rawreads /preads。0-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_file

a example string graph assembler that is desinged for handling diploid genomes

positional 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 for
                     assembling
  --min_idt MIN_IDT  minimum alignment identity of the reads to be considered
                     for assembling
  --lfc              use local flow constraint method rather than best overlap
                     method 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发行说明