生信软件工具箱

「三代组装」使用新版Falcon进行三代测序基因组组装

这里的新版指的是PacBio公司在2018年9月发布pb-assembly, 而这篇文章是在2018年9月30日发的。

今年早些时候在参加三代培训时，听说PacBio会在今年对Falcon进行一些改变。前几天我在读 readthedocs上的Falcon文档时，发现了文档页面上方出现了这样两栏提醒

注意

其中pb_assembly就是新的FALCON组装套装在GitHub上的使用文档，经过了几天的探索，我对它有了一点了解，写一篇教程作为官方文档的一些补充吧。

新版亮点

  1. 整合了串联重复序列和遍在重复序列的屏蔽（之前没有这一步）
  2. 以GFA格式存放graph文件，后续可以用Bandage进行可视化
  3. 通过算法和性能优化，提高了Associate Contigs的准确性
  4. 分析流程的性能优化

软件安装和数据准备

Falcon终于拥抱了bioconda, 这也就意味着我们再也不需要用到他们原本笨拙的安装脚本，浪费时间在安装软件上。

1
2
3
4
5
conda create -n pb-assembly  pb-assembly
source activate pb-assembly
# 或者
conda create -p ~/opt/biosoft/pb-assembly pb-assembly
source activate ~/opt/biosoft/pb-assembly

这里使用https://pb-falcon.readthedocs.io/en/latest/tutorial.html上的所用的E. coli数据集

1
2
wget https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz
tar -xvzf ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz

解压缩后的文件夹里有三个300M的Fasta文件, 将他们的实际路径记录到input.fofn中

1
2
3
ecoli.1.fasta
ecoli.2.fasta
ecoli.3.fasta

准备配置文件

为了进行组装，需要准备一个配置文件。我的配置文件为fc_run.cfg，内容如下。你们可以先预览一下，后面看我的解释说明。

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
#### Input
[General]
input_fofn=input.fofn
input_type=raw
pa_DBdust_option=
pa_fasta_filter_option=pass
target=assembly
skip_checks=False
LA4Falcon_preload=false

#### Data Partitioning
pa_DBsplit_option=-x500 -s50
ovlp_DBsplit_option=-x500 -s50

#### Repeat Masking
pa_HPCTANmask_option=
pa_REPmask_code=1,100;2,80;3,60

####Pre-assembly
genome_size=0
seed_coverage=20
length_cutoff=3000    
pa_HPCdaligner_option=-v -B128 -M24
pa_daligner_option=-e.7 -l1000 -k18 -h80 -w8 -s100
falcon_sense_option=--output-multi --min-idt 0.70 --min-cov 2 --max-n-read 800
falcon_sense_greedy=False

####Pread overlapping
ovlp_daligner_option=-e.96 -l2000 -k24 -h1024 -w6 -s100
ovlp_HPCdaligner_option=-v -B128 -M24

####Final Assembly
overlap_filtering_setting=--max-diff 100 --max-cov 300 --min-cov 2
fc_ovlp_to_graph_option=
length_cutoff_pr=2000

[job.defaults]
job_type=local
pwatcher_type=blocking
JOB_QUEUE=default
MB=32768
NPROC=6
njobs=32
submit = /bin/bash -c "${JOB_SCRIPT}" > "${JOB_STDOUT}" 2> "${JOB_STDERR}"

[job.step.da]
NPROC=4
MB=32768
njobs=32
[job.step.la]
NPROC=4
MB=32768
njobs=32
[job.step.cns]
NPROC=8
MB=65536
njobs=5
[job.step.pla]
NPROC=4
MB=32768
njobs=4
[job.step.asm]
NPROC=24
MB=196608
njobs=1

根据注释信息，文件分为”input”, “Data partitioning”, “Repeat Masking”, “Pre-assembly”, “Pread overlapping”, “Final Assembly”, 以及最后的任务调度部分，让我们分别看下这里面的内容

输入(Input)

1
2
3
4
5
6
7
8
9
#### Input
[General]
input_fofn=input.fofn
input_type=raw
pa_DBdust_option=
pa_fasta_filter_option=pass
target=assembly
skip_checks=False
LA4Falcon_preload=false

输入这里参数比较简单，基本不需要做任何改动，除了 pa_fasta_filter_options，用于处理一个ZMW（测序翻译孔）有多条subread时，到底选择哪一条的问题。

  • “pass”: 不做过滤，全部要。
  • “streamed-median”: 表示选择大于中位数的subread
  • “streamed-internal-median”: 当一个ZMW里的subread低于3条时选择最长，多于单条则选择大于中位数的subread

0.01版本pb-assembly的pa_DBdust_option有一个bug，也就是里面的参数不会传递给DBdust, DBdust是对read进行soft-masking，一般都用默认参数，因此这个bug问题不大。

数据分配(Data Partitioning)

1
2
pa_DBsplit_option=-x500 -s50
ovlp_DBsplit_option=-x500 -s50

这部分的设置会将参数传递给DBsplit，将数据进行拆分多个block，后续的并行计算都基于block。-s 50表示每个block大小为50M。 这适用于基因组比较小的物种，如果是大基因组则应该设置为-s 200或者-s 400

重复屏蔽(Repeat Masking)

1
2
3
#### Repeat Masking
pa_HPCTANmask_option=
pa_REPmask_code=1,100;2,80;3,60

屏蔽重复序列可以在不损失组装准确性的同时，提高后续组装的overlap/daligner步骤10~20倍速度，见Detecting and Masking Repeats.

pa_HPCTANmask_option的参数会传给串联重复步骤的HPCTANmask, 而pa_REPmask_code很复杂，它分为三次迭代，因此这里1:100;2,80;3,60 就表示第一次迭代检测每个block中出现超过100次的序列，第二次迭代将2个block合并一起检测超过80次的序列，第三次将3个block进行合并检测超过60次的序列。

https://github.com/PacificBiosciences/pbbioconda/issues/20

预组装（纠错）pre-assembly

1
2
3
4
5
6
7
8
####Pre-assembly
genome_size=0
seed_coverage=20
length_cutoff=3000    
pa_HPCdaligner_option=-v -B128 -M24
pa_daligner_option=-e.7 -l1000 -k18 -h80 -w8 -s100
falcon_sense_option=--output-multi --min-idt 0.70 --min-cov 2 --max-n-read 800
falcon_sense_greedy=False

当length_cutoff=-1时，设置genome_size和seed_coverage会自动计算要过滤的序列。否则是过滤低于一定长度的read。

pa_HPCdalinger_option参数不需要调整，-M表示每个进程的内存为24G，一般200M的block对应16G。

pa_daligner_option的参数比较重要：

  • -e:错误率，低质量序列设置为0.70，高质量设置为0.80。 值越高避免单倍型的坍缩
  • -l: 最低overlap的长度，文库比较短时为1000， 文库比较长为5000.
  • -k: 低质量数据为14，高质量数据为18
  • -h: 表示完全match的k-mer所覆盖的碱基数。和-l, -e有关，越大越严格。预组装时最大也不要超过最低overlap长度的1/4. 最低就设置为80

纠错后相互比对

1
2
3
####Pread overlapping
ovlp_daligner_option=-e.96 -l2000 -k24 -h1024 -w6 -s100
ovlp_HPCdaligner_option=-v -B128 -M24

和上面的参数类似，但是-e的范围调整为0.93-0.96，-l范围调整为1800-6000， -k调整为 18-24

最后组装

1
2
3
overlap_filtering_setting=--max-diff 100 --max-cov 300 --min-cov 2
fc_ovlp_to_graph_option=
length_cutoff_pr=2000

这里的参数就可以随便调整了，因为这一步速度很快。 例如length_cutoff_pr就可以从2000，提高到15000.

最后还有一部分是任务投递系统，如果是单节点运行，需要注意设置 njobs，这是同时投递的任务数。假如你将[job.step.cns]按照如下的方式设置，那么同时会出现 $ 8 \times 50 = 400 $ 个任务，如果你的内存只有128G，运行一段时间后你的所有内存就会被耗尽，那么基本上你就只能重启服务器了。

1
2
3
4
[job.step.cns]
NPROC=8
MB=65536
njobs=50

运行结果

用上述的配置文件，以fc_run fc_run.cfg运行后，最后的2-asm-falcon/p_ctg.fa的序列数有4条，最长为4,685,024, 之后我将length_cutoff_pr调整为15k，2-asm-falcon/p_ctg.fa序列只有一条，长度为4,638,411

下载asm.gfa到本地，用Bandage可视化，可以发现组装效果不错。

Bandage