使用Falcon对三代测序进行基因组组装

Falcon是PacBio公司开发的用于自家SMRT产出数据的基因组组装工具。Falcon分为三个部分：

HGAP：PacBio最先开发的工具，用于组装细菌基因组，名字缩写自Hierarchical genome-assembly process(层次基因组组装进程)。 适用于已知复杂度的基因组，且基因组大小不能超过3Gb. 由于是图形界面，所以用起来会非常方便。
Falcon：和HGAP工作流程相似，可认为是命令行版本的HGAP，能与Falcon-Unzip无缝衔接。
Falcon-Unzip: 适用于杂合度较高或者远亲繁殖或者是多倍体的物种
层次基因组组装过程(HGAP)分为两轮.

第一轮是选择种子序列或者是数据集中最长的序列(通过length_cufoff设置)，比较短的序列比对到长序列上用于产生高可信度的一致性序列。PacBio称其为预组装(pre-asembled), 其实和纠错等价。这一步可能会将种子序列在低覆盖度的区域进行分割(split)或者修整(trim)，由falcon_sense_options参数控制，最后得到preads(pre-assembled reads)。

第二轮是将preads相互比对，从而组装成contigs(contig指的是连续的不间断的基因组序列, contiguous sequence)

HGAP
基因组最后组装结果是单倍体，但实际上人类、动物和植物大部分的基因组都是二倍体，两套染色体之间或多或少存在的差异。这种差异在组装时就是“图”里的气泡(bubble)。PacBio开发的Falcon-Unzip就是用来处理“气泡”，把不同单倍型的contig分开。

Falcon-Unzip
运行参数分类
Falcon的运行非常简单，就是准备好配置文件传给fc_run.py，然后让fc_run.py调度所有需要的软件完成基因组组装即可。只不过初学者一开始可能会迷失在茫茫的参数中，所以我们要对参数进行划分，分层理解。

参数从是否直接参与基因组组装分为任务投递管理系统相关和实际组装相关。

任务投递管理系统相关参数：

任务管理系统类型：job_type,job_queue
不同阶段并发任务数：default_concurrent_jobs, pa_concurrent_jobs, cns_concurrent_jobs, ovlp_concurrent_jobs
不同阶段的投递参数：sge_option_da,sge_option_la, sge_option_cns, sge_option_pla, sge_option_fc
这些参数和实际的组装并没有多大关系，就是控制如何递交任务，一次递交多少任务而已。这些你需要根据实际计算机可用资源进行设置，Falcon推荐多计算节点任务方式。实际组装参数按照不同的任务阶段可以继续划分

原始序列间重叠检测和纠错: pa_DBsplit_option, pa_HPCdaligner_option,falcon_sense_option
纠错序列间重叠检测：ovlp_DBsplit_option, ovlp_HPCdaligner_option
字符串图组装: overlap_filtering_setting, length_cutoff_pr
除此之外，还有一些全局性的参数，如输入文件位置input_fofn, 输入数据类型input_type, 基因组预估大小genome_size以及只使用超过一定长度的序列用于组装 length_cutoff和length_cutoff_pr。那么问题来了，面对茫茫多的Falcon参数，我们到底应该如何设置才比较好？

当然是学习前人的经验，Falcon在https://pb-falcon.readthedocs.io/en/latest/parameters.html 提供了不同物种组装的参数设置参考文件， 我们通过比较不同配置间的参数差异来明确每个参数的意义，最后还需要通过实战了解。

组装实战
这里用的是 E. coli 数据集。由于三代组装对资源的消耗是非常厉害的，因此贸然使用较大基因组进行组装，都不知道什么时候才能把基因组装好，一不小心内存用尽服务器卡的只能重启，所以先用大肠杆菌先大致跑一下。

第一步： 准备数据
创建相应的目录，然后用wget进行数据下载并用tar进行解压缩。

mkdir -p assembly_test/pb-data && cd assembly_test/pb-data
wget https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz
tar -zxvf ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz

第二步：创建FOFN
FOFN指的是包含文件名的文件(file-of-file-names), 每一行里面都要有fasta文件的全路径.

/path/to/data/ecoli.1.subreads.fasta
/path/to/data/ecoli.2.subreads.fasta
/path/to/data/ecoli.3.subreads.fasta

第三步：准备配置文件
配置文件控制着Falcon组装的各个阶段所用的参数，然而一开始我们并不知道哪一个参数才是最优的，通常都需要不断的调整才行。当然由于目前已经有比较多的物种使用了Falcon进行组装，所以可以从他们的配置文件中进行借鉴

wget https://pb-falcon.readthedocs.io/en/latest/_downloads/fc_run_ecoli_local.cfg

该文件的大部分内容都不需要修改，除了如下几个参数

input_fofn: 这里的input.fofn就是第二步创建的文件，input.fofn。建议把该文件放在cfg文件的同级目录下，这样子就不需要改配置文件该文件的路径了。

# list of fasta files
input_fofn = input.fofn

genome_size, seed_coverage,length_cutoff,length-cutoff_pr 这三个参数控制纠错所用数据量和组装所用数据量. 如果要让程序在运行的时候自动确定用于纠错的数据量，就将length_cutoff设置成"-1"，同时设置基因组估计大小genome_size和用于纠错的深度seed_coverage。

# otherwise, user can specify length cut off in bp (e.g. 2000)
length_cutoff = -1
genome_size = 4652500
seed_coverage = 30
# The length cutoff used for overalpping the preassembled reads 
length_cutoff_pr = 12000

jobqueue: 这里用的是单主机而不是集群，所以其实随便取一个名字就行，但是对于SGE则要选择能够提交的队列名。

jobqueue = your_queue

xxx_concurrent_jobs: 同时运行的任务数。显然是越多越快，有些配置文件都写了192，但是对于大部分人而言是没有那么多资源资源的，盲目写多只会导致服务器宕机。

# job concurrency settings for...
# all jobs
default_concurrent_jobs = 30
# preassembly
pa_concurrent_jobs = 30
# consensus calling of preads
cns_concurrent_jobs = 30
# overlap detection
ovlp_concurrent_jobs = 30

参数参考: https://pb-falcon.readthedocs.io/en/latest/parameters.html

第四步：运行程序。
这一步可以写一个专门的shell脚本，然后用qsub(一种SGE集群任务投递命令）进行任务投递。作为一个单主机用户，我就直接在命令行中运行。

source ~/opt/biosoft/falcon_unzip/fc_env/bin/activate
fc_run fc_run_ecoli_local.cfg &> fc_run.log &

在程序运行过程中，可以通过几行命令来看下程序的进度。

检查不同阶段总共要执行的任务数目

# raw-data
grep '^#' 0-rawreads/run_jobs.sh
# preads
grep '^#' 1-preads_ovl/run_jobs.sh

当前已经执行的任务数

find 0-rawreads/ -name "job*done" | wc -l
find 0-rawreads/ -name "m_*done" | wc -l

评估组装运行结果
E. coli 的 subreads总共有1.01Gb数据，共105.451条reads。

我们可以用dazzler里的命令行来确认dazzler数据库是否构建正常

DBstats 0-rawreads/raw_reads.db | head -n 17

DBstats
从中我们可以发现，只有86.1%的read被用于构建dazzler数据库，这是由于配置文件里我们设置的DBsplit参数使用长度大于500bp的read，pa_DBsplit_option = -x500 -s50

此外，我们还需要检查 预组装(preassmebly)中实际用于组装的数据量

cat 0-rawreads/report/pre_assembly_stats.json

纠错数据报告
这里注意一点，这里的length_cutoff是程序根据配置文件中这一项length_cutoff = -1确定，这里的"-`"意味着程序会自动根据基因组大小和覆盖度确定，但是这不一定是最佳选择，你应该自己设置。

最后的数据存放在2-asm-falcon/文件夹下，简单用ls -lh 2-asm-falcon/查看该文件时，我发现组装的p_ctg.fa大小就只有1.4M， 而实际的大肠杆菌基因组大小为4.8M左右，差距有点大。

原因就是配置文件中length_cutoff的设置问题，上图显示只有22X用于组装。

当我将其修改成length_cutoff=2000时，纠错后还有156X数据用于组装，最后的组装的大小就成了4.8m。 因此用于纠错的数据应该尽量多一点，但是阈值也不要太低，否则增加了不必要的运算量。

由于大肠杆菌基因组小，1Gb的数据相当于深度为200X，深度非常的高，可以适当继续提高length_cutfoff 。当然大多数项目都是70X， 100X, 120X的数据， 设置成2000也就差不多了。

此外，建议多用几个长度阈值的preads进行组装，即修改配置文件中 length_cutoff_pr， 然后对"2-asm-falcon", "mypwatcher",和 "all.log" 重命名, 这样子重新运行脚本就会从纠错后序列开始。

我尝试了5k，10k，15k这三个参数，其中5k效果最差，装出了20条序列，基因组只有560k。10k装出了3条序列，基因组大小是4.5M，N501.7M。 15k效果最好，装出了一条4.6M的序列。

在length_cutoff分别是2k和15k情况下，length_cutoff_pr的长度都设置为15k，最后组装结果是length_cutoff 设置为2k时，最后序列长度更长

官方文档的坑
如果你在组装过程中用ps aux | grep 你的用户名 | wc -l 检查自己用了多少线程时。你会发现有一个时期会突然出现会出现好几百个python -m falcon_kit.mains.consensus 进程。

我在组装自己物种的时候出现了700多个，最后256G内存完全被消耗完，导致无法登陆服务器只好重启。

据我猜测应该是和falcon_sense_option中的--n_cores和cns_concurrent_jobs有关，我当时这两个参数分别设置了20和32，那么结果就要用到600多个进程。而每个进程大概会消耗1G内存，那么峰值就会是600G。

而这里组装 E. coli 基因组这两个参数分别是6和10，用free -h看内存消耗时也差不多时用了60多G内存。

所以一定要根据实际情况调整，配置文件中和并行运算的参数，量力而行，不然重启等着你。

本文作者：徐洲更