De novo 基因组组装套件 --- FALCON

发表于 更新于 分类于 ， 阅读次数： 本文字数： 7.2k 阅读时长 ≈ 7 分钟

Falcon：一组用于快速比对 long reads 以达到 consensus 和 组装的工具；是研究单倍体和二倍体基因组的一种高效组装算法。

Falcon 采用的是分层基因组组装策略（Hierarchical Genome Assembly Process），有如下几个过程：

  • Raw sub-reads overlapping for error correction (原始 sub-reads overlapping 用于纠错)
  • Pre-assembly and error correction (预组装和纠错)
  • Overlapping detection of the error corrected reads (纠错 reads 的重叠检测)
  • Overlap filtering (重叠过滤)
  • Constructing graph from overlaps (从 overlapping 中构造 graph)
  • Constructing contig from graph (从 graph 构造 contigs)

overlapping ： 重叠

graph： 图

Falcon 运行要求

程序要求在集群计算的环境下运行，所以需要一个用于长期脚本作业的作业队列和富含 cpu 的计算作业队列；同时需要计算队列满足相应的资源配置；官方建议的是组装一个 100Mb 大小的基因组，可能需要至少 32 个核心 cpu 和 128Gb RAM。

整个运行过程通过脚本 fc_run.py 统一进行提交，也就是说只需要将这个脚本提交到集群，之后这个脚本将会启动 Falcon 的整个工作流，从而将原始数据组装成初步的组装。

为了更好的适应多样的数据集，我们还需要准备一个配置文件 --- fc_run.cfg，用于根据输入数据集来控制所使用的计算资源和重要的算法参数以达到优化组装。

相关参数设置

在 fc_run.cfg 中：

  • input_fofn=input.fofn --- 指向包含所有输入数据的文件。
  • input_type=raw/preads --- 说明输入数据的类型；raw - 表示原始 reads；preads - 表示为进行过纠错的 reads，它将忽略任何纠错步骤并直接进入最终组装 contig 的步骤。
  • length_cutoff= --- 长读截止值；通常在基因组组装可以获得最长 15X 到 20X 的点上选择阈值会很好。但是，如果计算资源丰富，并且可能会发现纠错读取的其他应用，你可以设置较低的 length_cutoff 以获得更多的纠错读取。
  • pa_concurrent_jobs = --- 用来控制可以提交的并发作业的数量。
  • pa_DBsplit_option=-s400 --- 设置如何拆分序列数据库，对于大型基因组，应该设置为 -s400；-x 忽略长度小于指定阈值的 reads。
  • pa_HPCdaligner_option=-da1 --- 该 -da1 选项中的数字决定了在单个作业中与每个块进行比较的块数，较大的数字提供较大的工作，但总工作量较小。
  • falcon_sense_option --- 用于对预组装和纠错进行参数设置，--min_cov 控制 seed reads 何时因覆盖率低而被修剪或中断。--max_n_read 对用于纠错的读取次数设置上限。在高度重复的基因组中，你需要将其设置成更小，以确保共识代码不会浪费时间对齐重复。最长的适当重叠用于纠正以减少折叠重复的概率。
  • ovlp_DBsplit_option=-x500 -s200 --- 设置 OLC 组装过程中对 preads 数据进行分块的参数。
  • pa_HPCTANmask_option = -k18 -h480 -w8 -e.80 --- 设置在预组装过程中对串联重复进行屏蔽的参数，其参数会传递给 HPC.TANmask 命令。
  • pa_HPCREPmask_option=-k18 -h480 -w8 -e.80 --- 设置在预组装过程中对散在重复进行屏蔽的参数。迭代 3 次。HPC.TANmask 命令进行数据分析时有两个关键参数 -g 和 -c，用于设置分析的数据块数量和覆盖度阈值
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
## 该部分转载自 http://www.chenlianfu.com/?p=2755
[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

length_cutoff 控制在纠错过程中 length_cutoff_pr 使用的截至值，并控制用于后面的组装重叠步骤的截至值。在最终的组装中，更多的 rads 可能不会导致更好的组装，因为某些 reads 可能会铲射高噪音并在组装图中创建错误的连接。有时，尝试不同的值可能是有意义的，length_cutoff_pr 因为与第一个重叠步骤进行纠错相比，它的计算成本相对较低。一种策略是选择一个较小的 length_cutoff 并进行一次计算。稍后，我们可以使用不同的 length_cutoff_pr 的方法来获得更好的组装。

参考链接

[1] https://github.com/PacificBiosciences/FALCON/wiki/Manual

[2] http://www.chenlianfu.com/?p=2755