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

  • A+
所属分类:Genomics

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

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

Attention

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

新版亮点

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

软件安装和数据准备

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

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数据集

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

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

准备配置文件

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

#### 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)

#### 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)

pa_DBsplit_option=-x500 -s50
ovlp_DBsplit_option=-x500 -s50

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

重复屏蔽(Repeat Masking)

#### 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

####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_sizeseed_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

纠错后相互比对

####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

最后组装

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]按照如下的方式设置,那么同时会出现8X50=400个任务,如果你的内存只有128G,运行一段时间后你的所有内存就会被耗尽,那么基本上你就只能重启服务器了。

[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可视化,可以发现组装效果不错。

avatar

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: