falcon组装及polish

2019-10-17 本文已影响0人多啦A梦的时光机_648d

一：安装

$conda create -n denovo_asm

$source activate denovo_asm

$conda install pb-assembly

包含三个模块

FALCON: assembly pipeline
FALCON-Unzip: to phase the genome and perform phased-polishing with Arrow
FALCON-Phase: to extend phasing between unzipped haplotig blocks (requires HiC data)

Installed package recipes include:

- pb-falcon

- pb-dazzler

- genomicconsensus

- etc (all other dependencies)

FALCON和FALCON- unzip是用于PacBio长reads的从头组装（SMART sequence）。

FALCON

FALCON 是一个能识别二倍体的assembler。 FALCON产生了一组 primary contigs(p-contigs)作为主装配，一组 associate contigs(a-contigs)代表不同的等位基因变体。

FALCON- unzip

FALCON- unzip是一个真正的二倍体assembler，It takes the contigs from FALCON and phases the reads based on heterozygous SNPs identified in the initial assembly. 然后产生一组partially-phased primary contigs和fully-phased haplotigs(代表了represent divergent haplotypes).

FALCON-Phase

This method maps HiC data to the FALCON-Unzip assembly to fix phase switches between haplotigs within primary contigs.

注意：

* gcpp instead of GenomicConsensus for polishing (Falcon_unzip 1.2.0)

* Use all subreads for polishing (Falcon_unzip 1.1.5) We used to use only 1 per zmw, same as assembly typically. Chemistry v3+ has longer polymerase reads, resulting in multiple passes of library inserts in many cases. Read more here. Config: [Unzip]polish_include_zmw_all_subreads is "true"

* Use pbmm2 instead of blasr by default (Falcon_unzip 1.1.5) Config: [Unzip]polish_use_blasr = true for blasr

二：用法

assemble

$fc_run fc_run.cfg

unzip and polish

$fc_unzip.py fc_unzip.cfg

* Running with HiFi data:
$fc_unzip.py --target='ccs' fc_unzip.cfg

Extended phasing with HiC

$fc_phase.py fc_phase.cfg

FALCON和FALCON- unzip都将配置文件作为唯一的输入参数。

这里是约2.9Gb人类基因组装的fc_run.cfg

这里是一个fc_unzip.cfg文件

三：falcon组装

1.准备文件

数据文件(input.fofn)
输入文件是Pacbio测序得到的Fasta文件。将这些Fasta文件的路径写入到文件 input.fofn 中用于软件输入。

/data1/spider/project/NBFC20190790/3_C01/m64067_190909_054802.subreads.fasta
/data1/spider/project/NBFC20190790/4_D01/m64066_190916_034435.subreads.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 -s200
ovlp_DBsplit_option=-x500 -s200

#### 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.75 -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=23000

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

[job.step.da]
NPROC=4
MB=32768
njobs=8

[job.step.la]
NPROC=4
MB=32768
njobs=8

[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

配置文件解读：

[General]
设置任务提交方式
# jobtype的默认值应该是SGE，表示使用SGE集群进行计算；
# local表示使用本地主机运行FALCON，兼容性最好，毕竟集群的部署非常麻烦；
job_type = local

## 数据输入
# 设置输入文件为input.fofn，该文件中包含有PacBio数据的fasta文件。此外也可以输入dexta格式的文件。dexta格式是Pacbio测序结果h5的一种压缩结果，能极大压缩测序数据文件的大小。可以使用FALCON软件附带的命令dexta将fasta文件转换为dexta格式。
input_fofn = input.fofn
# 设置输入数据为原始测序数据，或者为修正后的数据。
input_type = raw
#input_type = preads

## 对PacBio数据进行校正：对Pacbio raw subreads进行overlapping、consensus和pre-assembly分析，对reads进行校正。
# 选择长度高于指定阈值的reads进行分析
length_cutoff = 12000
# 选择长度高于指定阈值的reads进行预组装
length_cutoff_pr = 12000
# 调用fasta2DB命令将reads数据分割成多份Blocks。
# -s50 参数表示每份数据含有50Mbp的数据量，该参数默认值为200，默认参数情况下，使用daligner进行比对需要消耗约16G内存。
# -x500 参数表示read长度低于500bp的会被忽略掉。
# -a 参数表示来同一个自零级波导孔的第二个read不会被忽略掉。
pa_DBsplit_option = -x500 -s50
# 调用daligner对所有Blocks进行Overlapping分析。
# -v 参数表示输出daligner程序运行信息
# -B4 参数表示每个daligner命令对4个Blocks进行Overlapping分析。该参数的值越大，则每个daligner命令的计算量越大，但是总的任务数越少。该参数等同于以前的-dal
# -k -w -h 参数设置匹配的kmer相关参数，其默认值分别为 14,6,35 。
# -T4 表示每个daligner使用4个线程进行计算，该参数默认值是4，该参数可以设置成2,4,8,16,32...。
# -M32 表示每个daligner命令使用32G内存，加入该参数起到限制内存使用的作用，对大基因组比较有用。
# -t16 参数表示过滤掉覆盖度高于16的kmer，这些kmer会导致内存使用过多。默认设置下，daligner可以根据-M参数的值自动计算本参数的值。
# -l1000 参数表示忽略长度低于1000bp的reads。
# -s1000 参数表示记录比对结果时以每1000bp为一个记录点，相比于默认值100，能少很多记录点。
# 使用daligner的默认参数能很好地处理raw pacbio数据。
# 而对corrected pacbio数据，推荐使用-k20 -h60 -e.85参数。
pa_HPCdaligner_option = -v -B4 -k14 -T8 -t16 -e.70 -l1000 -s1000
# FALCON使用fc_consensus.py调用C语言写的程序来根据daligner进行Overlapping分析的结果进行consensus分析，从而对subreads进行校正。
# --min_cov 参数设置当read序列上某位点覆盖度低于指定阈值时，对read进行打断或截短，默认值是6。
# --min_cov_aln 参数设置当read序列的平均覆盖度低于指定阈值时，直接略过该read，默认值是10。
# --min_n_read 和 --max_n_read 参数设定比对结果中包含的reads数在此范围内才能得到consensus结果，其默认值分别是10和500。对于基因组重复程度较高的情况，要设置较低的--max_n_read值来减少对重复区域进行consensus分析的计算消耗。
# --min_idt 参数设置最小identity的比对结果能用于reads校正。
# --n_core 参数设置允许的线程数，默认值是24。
falcon_sense_option = --output_multi --min_idt 0.70 --min_cov 4 --max_n_read 200 --n_core 8
# 设置运行daligner任务的并发数。注意的是daligner和fc_consensus.py任务本身可以多线程，因此总的计算需求是线程数*并发数。若是总的计算需要远远超过服务器的计算资源，容易导致宕机。
da_concurrent_jobs = 10
la_concurrent_jobs = 2
# 设置运行fc_consensus.py任务的并发数。注意的是daligner和fc_consensus.py任务本身可以多线程，因此总的计算需求是线程数*并发数。若是总的计算需要远远超过服务器的计算资源，容易导致宕机。
cns_concurrent_jobs = 20
# 设置在SGE集群运行的并发数
sge_option_da = -pe smp 8 -q jobqueue
sge_option_la = -pe smp 2 -q jobqueue
sge_option_fc = -pe smp 80 -q jobqueue
sge_option_cns = -pe smp 16 -q jobqueue

## 对校正后的reads进行overlapping分析，其参数和上一个步骤的参数一致。
ovlp_DBsplit_option = -x500 -s50
ovlp_HPCdaligner_option = -v -B4 -k20 -h60 -T8 -t32 -e.96 -l500 -s1000
# 设置对校正后reads运行daligner任务的并发数。
pda_concurrent_jobs = 10
pla_concurrent_jobs = 2
sge_option_pda = -pe smp 8 -q jobqueue
sge_option_pla = -pe smp 2 -q jobqueue

## 过滤overlaps
# 若reads首尾两端的覆盖度比平均覆盖度大很多，则表明reads首尾是重复序列；若reads首尾两端的覆盖度比平均覆盖度相差较小很多，则表明reads首尾出现错误的可能性很大。需要过滤掉这种reads的overlaps结果。该步骤的参数设置不对，容易导致overlaps全部被过滤掉，得不到基因组组装的结果。
# --bestn设置报告reads此数目的最优overlaps。
# --min_cov和--max_cov表示所允许的reads首尾的覆盖度范围。对于通过length_cutoff得到 >20x 校正后的数据进行的基因组组装，可以设置--min_cov值为5，设置--max_cov为平均覆盖度的3倍。若用于基因组组装的数据量较少，可以设置该值为1或2。
# --max_diff设置所允许的首尾覆盖度值的最大差异。设置该参数的值为平均覆盖度的2倍。
# 对于特殊情况，可以设置更高的--max_cov和--max_diff值。
# 可以使用在1-preads_ovl目录下运行 fc_ovlp_stats.py --fofn merge-gather/las.fofn 导出overlap结果首尾的覆盖度结果，从而帮助设置以上参数
overlap_filtering_setting = --max_diff 50 --max_cov 75 --min_cov 5 --bestn 10

## 工作分配
为了方便起见，下表列出了一些job_type配置和提交字符串示例。您可能需要修改一些参数来使用 job scheduler

local

bash -C ${CMD} >| ${STDOUT_FILE} 2>| ${STDERR_FILE}

sge

qsub -S /bin/bash -sync y -V -q myqueue -N ${JOB_NAME} -o "${JOB_STDOUT}" -e "${JOB_STDERR}" -pe smp ${NPROC} -l h_vmem=${MB}M "${JOB_SCRIPT}"

lsf

bsub -K -q myqueue -J ${JOB_NAME} -o ${JOB_STDOUT} -e ${JOB_STDERR} ${JOB_SCRIPT}

slurm

srun --wait=0 -p myqueue -J ${JOB_NAME} -o ${JOB_STDOUT} -e ${JOB_STDERR} --mem-per-cpu=${MB}M --cpus-per-task=${NPROC} ${JOB_SCRIPT}

例子：

Local

pwatcher_type = blocking

submit = /bin/bash -c "${JOB_SCRIPT}" > "${JOB_STDOUT}" 2> "${JOB_STDERR}"

SGE/qsub

submit =  qsub -S /bin/bash -sync y -V -q myqueue

  -N ${JOB_NAME}        \

  -o "${JOB_STDOUT}" \

  -e "${JOB_STDERR}" \

  -pe smp ${NPROC}    \

  -l h_vmem=${MB}M    \

  "${JOB_SCRIPT}"

四：falcon-unzip用于polish

falcon-unzip主要有两步：
3-unzip ：read alignment, SNP calling, read phasing, and diploid assembly of primary contigs and haplotigs

4-polish： phased polishing in which reads are used to polish in a haplotype-specific manner using BLASR and arrow

falcon配置非常简单， [General]部分中的第一个也是唯一一个设置是针对max_n_open_files的，默认设置为max_n_open_files=300。

与falcon类似，参数input_fofn里面写入fasta名称的输入文件的路径，另外如果希望 polish未压缩的基因组，还需要使用input_bam_fofn指定输入bam文件列表。

这里是一个fc_unzip.cfg文件

[General]

max_n_open_files = 1000

[Unzip]

input_fofn=input.fofn

input_bam_fofn=input_bam.fofn

polish_include_zmw_all_subreads = true

[job.defaults]

job_type=local

pwatcher_type=blocking

JOB_QUEUE=default

MB=32768

NPROC=6

njobs=8

submit = /bin/bash -c "${JOB_SCRIPT}" > "${JOB_STDOUT}" 2>  "${JOB_STDERR}"

[job.step.unzip.track_reads]

njobs=1

NPROC=48

MB=393216

# uses minimap2 now

[job.step.unzip.blasr_aln]

njobs=50

NPROC=2

MB=32000

[job.step.unzip.phasing]

njobs=100

NPROC=2

MB=16384

[job.step.unzip.hasm]

njobs=1

NPROC=48

MB=393216

# uses arrow now

[job.step.unzip.quiver]

njobs=50

NPROC=12

MB=98304

另外 [job.defaults]部分与falcon相同，唯一的区别是特定于FALCON-Unzip的job设置。

五：falcon-Phase

这里是falcon-Phase配置文件的列子