這裡的新版指的是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, 這也就意味著我們再也不需要用到他們原本笨拙的安裝指令碼，浪費時間在安裝軟體上。

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 

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]按照如下的方式設定，那麼同時會出現 $ 8 \times 50 = 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視覺化，可以發現組裝效果不錯。