使用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来减少垃圾评论。了解我们如何处理您的评论数据