使用FALCON对三代测序数据进行基因组组装

发表于由

1. FALCON软件简介

FALCON是PacBio公司开发的一款用于三代基因组De novo组装软件。相比于HGAP4软件，FALCON软件的基因组组装原理基本一致。但FALCON使用命令行运行，更适合于大基因组的组装，且能分析双倍体序列，并在基因组组装结果中给出包含变异位点信息的等位基因序列（alternative contigs / a-contigs）和主要的基因组序列（primary contig / p-contig）。每一条a-contig都有其对应的p-contig序列。因此，FALCON软件适合双倍体物种的基因组组装，能给出单倍的基因序列。其基因组组装结果中的p-contigs序列总长度要小于其它基因组组装软件（例如Canu和HGAP）的基因组序列。

FALCON-Unzip则是真正的单倍型组装软件，它能在FALCON或HGAP4软件的基因组组装结果基础上，利用较长的PacBio reads进行单倍型分析，对p-contigs序列向单倍型进行转换，同时输出单倍型序列（haplotig）区块。

2. FALCON软件下载与安装

FALCON软件在2018年在算法上有了较大的更改，将FALCON在内的多个软件整合到了pb-assembly软件（https://github.com/PacificBiosciences/pb-assembly）中。同时提供快捷的bioconda安装方法：

$ wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
$ sh Miniconda3-latest-Linux-x86_64.sh
    进入交互式界面：首先，输入yes，按Enter键，表示同意软件的license；然后，输入conda程序安装路径/opt/biosoft/miniconda3_for_pb-assembly，按enter键； 再设置是否添加Minicoda3到PATH变量，直接按Enter键，表示选择no，以免和系统自带的软件冲突。若需要切换到miniconda3环境，输入命令export PATH=/opt/biosoft/miniconda3_for_pb-assembly/bin:$PATH即可。
$ export PATH=/opt/biosoft/miniconda3_for_pb-assembly/bin:$PATH
$ conda config --add channels defaults
$ conda config --add channels bioconda
$ conda config --add channels conda-forge
    配置bioconda channel（https://bioconda.github.io/index.html）
$ conda install pb-assembly
    安装pb-assembly软件
$ source activate /opt/biosoft/miniconda3_for_pb-assembly
    激活FALCON软件的运行环境

3. FALCON软件的使用

FALCON软件运行示例：

首先，载入FALCON软件的运行环境。

然后，准备三种输入文件：测序数据FASTA格式文件、包含测序数据文件路径的FOFN列表文件、程序运行的参数配置文件。下载基因组大小为4.6Mb的E. coli基因组PacBio测序数据。
$ curl -L https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139\
_c100713652400000001823152404301535_s1_p0.subreads.tar | tar xvf -
$ ls */ecoli.?.fasta > input.fofn
$ wget https://pb-falcon.readthedocs.io/en/latest/_downloads/fc_run_ecoli_local.cfg
当使用最新版的FALCON进行组装时，需要对此配置文件进行一些修改。

最后，对E. coli基因组进行组装
$ fc_run.py fc_run_ecoli_local.cfg

生成的最终主要结果文件为 2-asm-falcon/p_ctg.fa 。

运行FALCON最关键的步骤是准备配置文件fc_run.cfg。一个配置文件示例及其参数讲解如下： 

[General]
输入参数：
input_fofn = input.fofn
设置数据输入文件为input.fofn，该文本文件中每行是一个PacBio测序数据的FASTA文件路径。可以使用相对路径。
input_type = raw
设置输入数据类型为原始测序数据（raw），或者为修正后的数据（preads）。

对数据进行分块：
pa_DBsplit_option = -x500 –s200
设置预组装（pre-assembly）过程中对raw reads数据进行分块的参数。这些参数会传递给Dbsplit命令来执行数据分块。
-s<int>    default: 200
表示每份数据含有50Mb的数据量，该参数默认值为200。默认参数情况下，使用daligner进行比对需要消耗约16G内存。推荐设置该参数的值为预测的基因组大小。此外，软件作者推荐对小于10Mb的基因组设置该参数值为50，对于大基因组设置该参数值为200。
-x<int>
忽略长度小于指定阈值的reads。
ovlp_DBsplit_option=-x500 –s200
设置OLC组装过程中对preads数据进行分块的参数。

对重复序列进行屏蔽：
pa_HPCTANmask_option = -k18 -h480 -w8 -e.80
设置在预组装过程对串联重复进行屏蔽的参数。其参数会传递给HPC.TANmask命令。该命令是正常daligner命令的一个变种，适合对串联重复序列进行比对。
pa_HPCREPmask_option = -k18 -h480 -w8 -e.80
设置在预组装过程对散在重复进行屏蔽的参数。其参数会传递给HPC.REPmask命令。该分析步骤会迭代3次。HPC.TANmask命令进行数据分析时有两个关键参数-g和-c，用于设置分析的数据块数量和覆盖度阈值。pa_REPmask_code设置3次迭代过程中这两个参数的值。这两个参数值的设置与数据分块方法和测序数据量有关。
pa_REPmask_code = 1,666;4,266;8,53
设置三次HPC.REPmask命令中的-g和-c参数的值。以上参数值则表示不对散在重复进行屏蔽。该参数的设置可以参考网址https://dazzlerblog.wordpress.com/2016/04/01/detecting-and-soft-masking-repeats/。在该文章中，作者给出了一个示例：对30Gb的基因组测序30x并将数据分割成3600个大小为250Mb的数据块，使用参数值“1,20;10,15;100,10”。每个数据块的测序深度为0.008x，对1个数据块的数据进行比对后，对reads中覆盖度超过20x的位点进行屏蔽，即对高度重复（20 / 0.008 = 2500x）的基因组序列位点进行屏蔽；然后第二次迭代中，对10个数据块的数据进行比对，对reads中覆盖度超过15x的位点进行屏蔽，即对中度重复（15 / 0.08 = 187.5x）的基因组序列位点进行屏蔽；最后第三次迭代中，对100个数据块的数据进行比对，对reads中覆盖度超过10x的位点进行屏蔽，即对低度重复（10 / 0.8 = 12.5x）的基因组序列位点进行屏蔽。示例中基因组大小是30Gb，超级大，一般含有大量重复序列，使用上述参数对数据中的重复位点进行屏蔽。
对于常规的基因组组装，设置pa_REPmask_code参数后，在三次迭代过程中，推荐逐步对重复次数为1000x、100x和10x的序列位点进行屏蔽。例如：对300Mb大小的基因组测序100x，将数据分割成150个大小为200Mb的数据块，可以考虑设置pa_REPmask_code值为“1,666;4,266;8,53”。
ovlp_HPCTANmask_option = -k20 -h480 -w8 -e.80
设置在OLC组装过程对串联重复进行屏蔽的参数。其参数会传递给HPC.TANmask命令。

预组装（pre-assembly）参数：
length_cutoff = -1
设置对种子序列的长度筛选阈值。若该参数值设置为-1，则程序自动计算种子序列的长度筛选阈值，根据如下两个参数挑选最长的reads序列作为种子序列，直到其数据量达到基因组指定覆盖度为止。
genome_size = 4652500
seed_coverage = 30
设置基因组大小和种子序列覆盖度。推荐设置seed_coverage参数的值为20~40x。当然是用更多的种子序列有利于基因组的组装效果，但是会消耗更多运行时间。
pa_daligner_option = -k14 -w6 -h35 -e.70 -l1000
在新版本FALCON中该参数额外再分出了pa_HPCdaligner_option参数，专门用于对数据块的分配。pa_daligner_option参数专门用于调用daligner软件对raw reads进行比对，进行overlap分析。其常用参数：
-k<int>    default: 14
设置k-mer大小。推荐设置为14~18。设置较低的值，能提高sensitivity，得到更多的reads重叠结果，但增加磁盘、内存、计算和时间消耗，适合数据质量较差的情况；设置较高的值，能提高specificity，系统资源消耗更少，运行速度更快，仅适合数据质量较高的情况。
-w<int>    default: 6
daligner命令对两个reads进行比对，得到一个斜对角线的比对带。该带的宽度设置为2的w次方。
-h<int>    default: 35
在斜对角线带的比对区域上，需要至少有指定数目的碱基能被k-mers覆盖。-k、-w和-h是daligner命令进行比对的主要参数。默认设置下，寻找两两reads之间重叠，先将reads分割成长度为14bp的k-mers序列；对其中一条序列宽度为64bp的区域进行分析，要求至少有35个碱基位点能和被另外一条序列的k-mers覆盖，则找到重叠。
-k、-w和-h的默认参数适合于PacBio测序的raw reads。对于修正后的read，推荐使用参数“daligner -k20 -h60 -e.95”。
-e<float>    default: 0.70
设置identity阈值。推荐设置为0.70~0.80。设置为较高的值有利于单倍型序列的组装。
-l<int>    default: 1000
忽略长度小于低于设定值的reads。
-T<int>    default: 4
设置每个daligner使用的CPU线程数。
-M<int>
尝试仅使用指定大小GB的内存。程序通过忽略覆盖度过高的k-mers来降低内存使用。一般情况下，对
pa_HPCdaligner_option = -v -B4 -M16
程序调用HPC.daligner进行分析，相比于daligner命令，其增加的常用参数：
-v
将-v参数传递给LAsort和LAmerge命令。
-B<int>    default: 4
设置每个daligner任务对指定数目的数据块进行分析。在最新版本的FALCON中，该参数不能再设置到pa_daligner_option参数中，否则程序运行出错。
falcon_sense_option = --output-multi --min-idt 0.70 --min-cov 4 --max-n-read 200
FALCON使用fc_consensus命令根据daligner进行Overlapping分析的结果进行一致性分析，从而对种子序列进行校正。常用参数：
--output-multi
输出每个种子序列多个修正后的区域。
--min-cov    default: 6
设置当种子序列上某位点覆盖度低于指定阈值时，则在该位点对序列打断或截短。
--min-cov-aln    default: 10
设置当种子序列的平均覆盖度低于指定阈值时，过滤掉该序列。
--min-n-read    default: 10
--max-n-read    default: 500
对覆盖度在--min_n_read和--max_n_read两个参数设定范围内的位点进行一致性分析，从而对种子序列进行校正。
--min-idt    default: 0.7
设置比对结果中reads重叠部分的最小identity阈值。
--n-core    default: 24
设置运行的线程数，默认值是24。
falcon_sense_greedy = False
falcon_sense_skip_contained = False

OLC算法进行组装的参数：
length_cutoff_pr = 12000
对预组装结果中长度超过此阈值的reads使用OLC算法进行下一步组装。
ovlp_HPCdaligner_option = -v -B4 -M16
ovlp_daligner_option = -k20 -w6 -h60 -e.96 -l1000
调用daligner对校正后的preads进行overlapping分析时，参数相对要设置更加严格。
-k<int>    default: 14
对preads进行比对，推荐该参数的值为18~24。
-e<float>    default: 0.70
对preads进行比对，推荐该参数的值为0.93~0.96。
overlap_filtering_setting = --max-diff 100 --max-cov 100 --min-cov 2 --bestn 10
过滤overlaps。若reads首尾两端的覆盖度比平均覆盖度大很多，则表明reads首尾是重复序列；若reads首尾两端的覆盖度比平均覆盖度相差较小很多，则表明reads首尾出现错误的可能性很大。需要过滤掉这种reads的overlaps结果。
--bestn    default: 10
设置报告reads此数目的最优overlaps。
--min-cov和—max-cov
所允许的reads首尾的覆盖度范围。
--max-diff
设置所允许的首尾覆盖度值的最大差异。
fc_ovlp_to_graph_option = --min_len 4000 --min_idt 0.96

任务投递与计算资源消耗参数：
[job.defaults]
job_type = local
设置任务提交方式。local表示使用本地主机运行FALCON，运行相对稳定。也可以设置为集群任务提交方式，可以设置的值有：sge、lsf、pbs、torque和slurm。该参数设置不同的值，后面的submit也要设置对应的信息。
pwatcher_type=blocking
submit = bash -C ${CMD} >| ${STDOUT_FILE} 2>| ${STDERR_FILE}
MB=32768
NPROC=6
njobs=32
FALCON软件流程对数据进行了切割，再进行并行化运算，加快程序运行速度。MB设置任务实用的内存数；NPROC设置单个任务的CPU线程数；njobs设置并行任务数。后续分别对每个流程中的计算资源进行分配。da表示预组装过程中运行的daligner命令；la表示预组装过程中运行的Lasort和LAmerge命令；cns表示预组装过程中运行的fc_consensus命令；pda表示运行OLC组装过程中运行的的daligner命令；pla表示OLC组装过程中运行的Lasort和LAmerge命令；asm表示最终组装过程。
[job.step.da]
NPROC=4
MB=32768
njobs=240
[job.step.la]
NPROC=4
MB=32768
njobs=240
[job.step.cns]
NPROC=8
MB=65536
njobs=240
[job.step.pda]
NPROC=4
MB=32768
njobs=240
[job.step.pla]
NPROC=4
MB=32768
njobs=240
[job.step.asm]
NPROC=24
MB=196608
njobs=1

FALCON的结果文件：

0-rawreads/
该目录存放对raw subreads进行overlpping分析与校正的结果；
0-rawreads/cns-runs/cns_*/*/*.fasta存放校正后的序列信息。

1-preads_ovl/
该目录存放对校正后reads进行overlapping的结果；

2-asm-falcon/
该目录是最终结果目录，主要的结果文件是p_ctg.fa和a_ctg.fa。

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据。