Falcon是PacBio公司開發(fā)的用于自家SMRT產(chǎn)出數(shù)據(jù)的基因組組裝工具。Falcon分為三個部分：
- HGAP：PacBio最先開發(fā)的工具，用于組裝細菌基因組，名字縮寫自Hierarchical genome-assembly process(層次基因組組裝進程)。 適用于已知復(fù)雜度的基因組，且基因組大小不能超過3Gb. 由于是圖形界面，所以用起來會非常方便。
- Falcon：和HGAP工作流程相似，可認為是命令行版本的HGAP，能與Falcon-Unzip無縫銜接。
- Falcon-Unzip: 適用于雜合度較高或者遠親繁殖或者是多倍體的物種
層次基因組組裝過程(HGAP)分為兩輪.
第一輪是選擇種子序列或者是數(shù)據(jù)集中最長的序列(通過length_cufoff設(shè)置)，比較短的序列比對到長序列上用于產(chǎn)生高可信度的一致性序列。PacBio稱其為預(yù)組裝(pre-asembled), 其實和糾錯等價。這一步可能會將種子序列在低覆蓋度的區(qū)域進行分割(split)或者修整(trim)，由falcon_sense_options參數(shù)控制，最后得到preads(pre-assembled reads)。
第二輪是將preads相互比對，從而組裝成contigs(contig指的是連續(xù)的不間斷的基因組序列, contiguous sequence)
基因組最后組裝結(jié)果是單倍體，但實際上人類、動物和植物大部分的基因組都是二倍體，兩套染色體之間或多或少存在的差異。這種差異在組裝時就是“圖”里的氣泡(bubble)。PacBio開發(fā)的Falcon-Unzip就是用來處理“氣泡”，把不同單倍型的contig分開。
運行參數(shù)分類
Falcon的運行非常簡單，就是準備好配置文件傳給fc_run.py，然后讓fc_run.py調(diào)度所有需要的軟件完成基因組組裝即可。只不過初學(xué)者一開始可能會迷失在茫茫的參數(shù)中，所以我們要對參數(shù)進行劃分，分層理解。
參數(shù)從是否直接參與基因組組裝分為任務(wù)投遞管理系統(tǒng)相關(guān)和實際組裝相關(guān)。
任務(wù)投遞管理系統(tǒng)相關(guān)參數(shù)：
- 任務(wù)管理系統(tǒng)類型：
job_type,job_queue - 不同階段并發(fā)任務(wù)數(shù)：
default_concurrent_jobs,pa_concurrent_jobs,cns_concurrent_jobs,ovlp_concurrent_jobs - 不同階段的投遞參數(shù)：
sge_option_da,sge_option_la,sge_option_cns,sge_option_pla,sge_option_fc
這些參數(shù)和實際的組裝并沒有多大關(guān)系，就是控制如何遞交任務(wù)，一次遞交多少任務(wù)而已。這些你需要根據(jù)實際計算機可用資源進行設(shè)置，F(xiàn)alcon推薦多計算節(jié)點任務(wù)方式。實際組裝參數(shù)按照不同的任務(wù)階段可以繼續(xù)劃分
- 原始序列間重疊檢測和糾錯:
pa_DBsplit_option,pa_HPCdaligner_option,falcon_sense_option - 糾錯序列間重疊檢測：
ovlp_DBsplit_option,ovlp_HPCdaligner_option - 字符串圖組裝:
overlap_filtering_setting,length_cutoff_pr
除此之外，還有一些全局性的參數(shù)，如輸入文件位置input_fofn, 輸入數(shù)據(jù)類型input_type, 基因組預(yù)估大小genome_size以及只使用超過一定長度的序列用于組裝 length_cutoff和length_cutoff_pr。那么問題來了，面對茫茫多的Falcon參數(shù)，我們到底應(yīng)該如何設(shè)置才比較好？
當然是學(xué)習(xí)前人的經(jīng)驗，F(xiàn)alcon在https://pb-falcon.readthedocs.io/en/latest/parameters.html 提供了不同物種組裝的參數(shù)設(shè)置參考文件， 我們通過比較不同配置間的參數(shù)差異來明確每個參數(shù)的意義，最后還需要通過實戰(zhàn)了解。
組裝實戰(zhàn)
這里用的是 E. coli 數(shù)據(jù)集。由于三代組裝對資源的消耗是非常厲害的，因此貿(mào)然使用較大基因組進行組裝，都不知道什么時候才能把基因組裝好，一不小心內(nèi)存用盡服務(wù)器卡的只能重啟，所以先用大腸桿菌先大致跑一下。
第一步： 準備數(shù)據(jù)
創(chuàng)建相應(yīng)的目錄，然后用wget進行數(shù)據(jù)下載并用tar進行解壓縮。
mkdir -p assembly_test/pb-data && cd assembly_test/pb-data
wget https://downloads.pacbcloud.com/public/data/git-sym/ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz
tar -zxvf ecoli.m140913_050931_42139_c100713652400000001823152404301535_s1_p0.subreads.tar.gz
第二步：創(chuàng)建FOFN
FOFN指的是包含文件名的文件(file-of-file-names), 每一行里面都要有fasta文件的全路徑.
/path/to/data/ecoli.1.subreads.fasta
/path/to/data/ecoli.2.subreads.fasta
/path/to/data/ecoli.3.subreads.fasta
第三步：準備配置文件
配置文件控制著Falcon組裝的各個階段所用的參數(shù)，然而一開始我們并不知道哪一個參數(shù)才是最優(yōu)的，通常都需要不斷的調(diào)整才行。當然由于目前已經(jīng)有比較多的物種使用了Falcon進行組裝，所以可以從他們的配置文件中進行借鑒
wget https://pb-falcon.readthedocs.io/en/latest/_downloads/fc_run_ecoli_local.cfg
該文件的大部分內(nèi)容都不需要修改，除了如下幾個參數(shù)
input_fofn: 這里的input.fofn就是第二步創(chuàng)建的文件，input.fofn。建議把該文件放在cfg文件的同級目錄下，這樣子就不需要改配置文件該文件的路徑了。
# list of fasta files
input_fofn = input.fofn
genome_size, seed_coverage,length_cutoff,length-cutoff_pr 這三個參數(shù)控制糾錯所用數(shù)據(jù)量和組裝所用數(shù)據(jù)量. 如果要讓程序在運行的時候自動確定用于糾錯的數(shù)據(jù)量，就將length_cutoff設(shè)置成"-1"，同時設(shè)置基因組估計大小genome_size和用于糾錯的深度seed_coverage。
# otherwise, user can specify length cut off in bp (e.g. 2000)
length_cutoff = -1
genome_size = 4652500
seed_coverage = 30
# The length cutoff used for overalpping the preassembled reads
length_cutoff_pr = 12000
jobqueue: 這里用的是單主機而不是集群，所以其實隨便取一個名字就行，但是對于SGE則要選擇能夠提交的隊列名。
jobqueue = your_queue
xxx_concurrent_jobs: 同時運行的任務(wù)數(shù)。顯然是越多越快，有些配置文件都寫了192，但是對于大部分人而言是沒有那么多資源資源的，盲目寫多只會導(dǎo)致服務(wù)器宕機。
# job concurrency settings for...
# all jobs
default_concurrent_jobs = 30
# preassembly
pa_concurrent_jobs = 30
# consensus calling of preads
cns_concurrent_jobs = 30
# overlap detection
ovlp_concurrent_jobs = 30
參數(shù)參考: https://pb-falcon.readthedocs.io/en/latest/parameters.html
第四步：運行程序。
這一步可以寫一個專門的shell腳本，然后用qsub(一種SGE集群任務(wù)投遞命令）進行任務(wù)投遞。作為一個單主機用戶，我就直接在命令行中運行。
source ~/opt/biosoft/falcon_unzip/fc_env/bin/activate
fc_run fc_run_ecoli_local.cfg &> fc_run.log &
在程序運行過程中，可以通過幾行命令來看下程序的進度。
檢查不同階段總共要執(zhí)行的任務(wù)數(shù)目
# raw-data
grep '^#' 0-rawreads/run_jobs.sh
# preads
grep '^#' 1-preads_ovl/run_jobs.sh
當前已經(jīng)執(zhí)行的任務(wù)數(shù)
find 0-rawreads/ -name "job*done" | wc -l
find 0-rawreads/ -name "m_*done" | wc -l
評估組裝運行結(jié)果
E. coli 的 subreads總共有1.01Gb數(shù)據(jù)，共105.451條reads。
我們可以用dazzler里的命令行來確認dazzler數(shù)據(jù)庫是否構(gòu)建正常
DBstats 0-rawreads/raw_reads.db | head -n 17
從中我們可以發(fā)現(xiàn)，只有86.1%的read被用于構(gòu)建dazzler數(shù)據(jù)庫，這是由于配置文件里我們設(shè)置的DBsplit參數(shù)使用長度大于500bp的read，pa_DBsplit_option = -x500 -s50
此外，我們還需要檢查 預(yù)組裝(preassmebly)中實際用于組裝的數(shù)據(jù)量
cat 0-rawreads/report/pre_assembly_stats.json
這里注意一點，這里的length_cutoff是程序根據(jù)配置文件中這一項length_cutoff = -1確定，這里的"-`"意味著程序會自動根據(jù)基因組大小和覆蓋度確定，但是這不一定是最佳選擇，你應(yīng)該自己設(shè)置。
最后的數(shù)據(jù)存放在2-asm-falcon/文件夾下，簡單用ls -lh 2-asm-falcon/查看該文件時，我發(fā)現(xiàn)組裝的p_ctg.fa大小就只有1.4M， 而實際的大腸桿菌基因組大小為4.8M左右，差距有點大。
原因就是配置文件中length_cutoff的設(shè)置問題，上圖顯示只有22X用于組裝。
當我將其修改成length_cutoff=2000時，糾錯后還有156X數(shù)據(jù)用于組裝，最后的組裝的大小就成了4.8m。 因此用于糾錯的數(shù)據(jù)應(yīng)該盡量多一點，但是閾值也不要太低，否則增加了不必要的運算量。
由于大腸桿菌基因組小，1Gb的數(shù)據(jù)相當于深度為200X，深度非常的高，可以適當繼續(xù)提高
length_cutfoff。當然大多數(shù)項目都是70X， 100X, 120X的數(shù)據(jù)， 設(shè)置成2000也就差不多了。
此外，建議多用幾個長度閾值的preads進行組裝，即修改配置文件中 length_cutoff_pr， 然后對"2-asm-falcon", "mypwatcher",和 "all.log" 重命名, 這樣子重新運行腳本就會從糾錯后序列開始。
我嘗試了5k，10k，15k這三個參數(shù)，其中5k效果最差，裝出了20條序列，基因組只有560k。10k裝出了3條序列，基因組大小是4.5M，N501.7M。 15k效果最好，裝出了一條4.6M的序列。
在length_cutoff分別是2k和15k情況下，length_cutoff_pr的長度都設(shè)置為15k，最后組裝結(jié)果是length_cutoff 設(shè)置為2k時，最后序列長度更長
官方文檔的坑
如果你在組裝過程中用ps aux | grep 你的用戶名 | wc -l 檢查自己用了多少線程時。你會發(fā)現(xiàn)有一個時期會突然出現(xiàn)會出現(xiàn)好幾百個python -m falcon_kit.mains.consensus 進程。
我在組裝自己物種的時候出現(xiàn)了700多個，最后256G內(nèi)存完全被消耗完，導(dǎo)致無法登陸服務(wù)器只好重啟。
據(jù)我猜測應(yīng)該是和falcon_sense_option中的--n_cores和cns_concurrent_jobs有關(guān)，我當時這兩個參數(shù)分別設(shè)置了20和32，那么結(jié)果就要用到600多個進程。而每個進程大概會消耗1G內(nèi)存，那么峰值就會是600G。
而這里組裝 E. coli 基因組這兩個參數(shù)分別是6和10，用free -h看內(nèi)存消耗時也差不多時用了60多G內(nèi)存。
所以一定要根據(jù)實際情況調(diào)整，配置文件中和并行運算的參數(shù)，量力而行，不然重啟等著你。
