这里的新版指的是PacBio公司在2018年9月发布pb-assembly, 而这篇文章是在2018年9月30日发的。
今年早些时候在参加三代培训时,听说PacBio会在今年对Falcon进行一些改变。前几天我在读 readthedocs上的Falcon文档时,发现了文档页面上方出现了这样两栏提醒
其中pb_assembly就是新的FALCON组装套装在GitHub上的使用文档,经过了几天的探索,我对它有了一点了解,写一篇教程作为官方文档的一些补充吧。
新版亮点
- 整合了串联重复序列和遍在重复序列的屏蔽(之前没有这一步)
- 以GFA格式存放graph文件,后续可以用Bandage进行可视化
- 通过算法和性能优化,提高了Associate Contigs的准确性
- 分析流程的性能优化
软件安装和数据准备
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次的序列。
预组装(纠错)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_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
纠错后相互比对
####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]按照如下的方式设置,那么同时会出现 个任务,如果你的内存只有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可视化,可以发现组装效果不错。