天天看點

使用Falcon對三代測序進行基因組組裝

Falcon

是PacBio公司開發的用于自家SMRT産出資料的基因組組裝工具。Falcon分為三個部分:

HGAP:PacBio最先開發的工具,用于組裝細菌基因組,名字縮寫自Hierarchical genome-assembly process(層次基因組組裝程序)。 适用于已知複雜度的基因組,且基因組大小不能超過3Gb. 由于是圖形界面,是以用起來會非常友善。

Falcon:和HGAP工作流程相似,可認為是指令行版本的HGAP,能與Falcon-Unzip無縫銜接。

Falcon-Unzip: 适用于雜合度較高或者遠親繁殖或者是多倍體的物種

層次基因組組裝過程(HGAP)分為兩輪.

第一輪是選擇種子序列或者是資料集中最長的序列(通過length_cufoff設定),比較短的序列比對到長序列上用于産生高可信度的一緻性序列。PacBio稱其為預組裝(pre-asembled), 其實和糾錯等價。這一步可能會将種子序列在低覆寫度的區域進行分割(split)或者修整(trim),由falcon_sense_options參數控制,最後得到preads(pre-assembled reads)。

第二輪是将preads互相比對,進而組裝成contigs(contig指的是連續的不間斷的基因組序列, contiguous sequence)

使用Falcon對三代測序進行基因組組裝

HGAP

基因組最後組裝結果是單倍體,但實際上人類、動物和植物大部分的基因組都是二倍體,兩套染色體之間或多或少存在的差異。這種差異在組裝時就是“圖”裡的氣泡(bubble)。PacBio開發的Falcon-Unzip就是用來處理“氣泡”,把不同單倍型的contig分開。

使用Falcon對三代測序進行基因組組裝

Falcon-Unzip

運作參數分類

Falcon的運作非常簡單,就是準備好配置檔案傳給fc_run.py,然後讓fc_run.py排程所有需要的軟體完成基因組組裝即可。隻不過初學者一開始可能會迷失在茫茫的參數中,是以我們要對參數進行劃分,分層了解。

參數從是否直接參與基因組組裝分為任務投遞管理系統相關和實際組裝相關。

任務投遞管理系統相關參數:

任務管理系統類型:

job_type,job_queue

不同階段并發任務數:

default_concurrent_jobs

,

pa_concurrent_jobs

cns_concurrent_jobs

ovlp_concurrent_jobs

不同階段的投遞參數:

sge_option_da,sge_option_la

sge_option_cns

sge_option_pla

sge_option_fc

這些參數和實際的組裝并沒有多大關系,就是控制如何遞交任務,一次遞交多少任務而已。這些你需要根據實際計算機可用資源進行設定,Falcon推薦多計算節點任務方式。實際組裝參數按照不同的任務階段可以繼續劃分

原始序列間重疊檢測和糾錯:

pa_DBsplit_option

pa_HPCdaligner_option,falcon_sense_option

糾錯序列間重疊檢測:

ovlp_DBsplit_option

ovlp_HPCdaligner_option

字元串圖組裝:

overlap_filtering_setting

length_cutoff_pr

除此之外,還有一些全局性的參數,如輸入檔案位置

input_fofn

, 輸入資料類型

input_type

, 基因組預估大小

genome_size

以及隻使用超過一定長度的序列用于組裝

length_cutoff

length_cutoff_pr

。那麼問題來了,面對茫茫多的Falcon參數,我們到底應該如何設定才比較好?

當然是學習前人的經驗,Falcon在

https://pb-falcon.readthedocs.io/en/latest/parameters.html

提供了不同物種組裝的參數設定參考檔案, 我們通過比較不同配置間的參數差異來明确每個參數的意義,最後還需要通過實戰了解。

組裝實戰

這裡用的是 E. coli 資料集。由于三代組裝對資源的消耗是非常厲害的,是以貿然使用較大基因組進行組裝,都不知道什麼時候才能把基因組裝好,一不小心記憶體用盡伺服器卡的隻能重新開機,是以先用大腸杆菌先大緻跑一下。

第一步: 準備資料

建立相應的目錄,然後用

wget

進行資料下載下傳并用

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           

第二步:建立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組裝的各個階段所用的參數,然而一開始我們并不知道哪一個參數才是最優的,通常都需要不斷的調整才行。當然由于目前已經有比較多的物種使用了Falcon進行組裝,是以可以從他們的配置檔案中進行借鑒

wget https://pb-falcon.readthedocs.io/en/latest/_downloads/fc_run_ecoli_local.cfg           

該檔案的大部分内容都不需要修改,除了如下幾個參數

input_fofn

: 這裡的input.fofn就是第二步建立的檔案,input.fofn。建議把該檔案放在cfg檔案的同級目錄下,這樣子就不需要改配置檔案該檔案的路徑了。

# list of fasta files
input_fofn = input.fofn           

genome_size

seed_coverage

length_cutoff,length-cutoff_pr

這三個參數控制糾錯所用資料量群組裝所用資料量. 如果要讓程式在運作的時候自動确定用于糾錯的資料量,就将

length_cutoff

設定成"-1",同時設定基因組估計大小

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

: 同時運作的任務數。顯然是越多越快,有些配置檔案都寫了192,但是對于大部分人而言是沒有那麼多資源資源的,盲目寫多隻會導緻伺服器當機。

# 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           
參數參考:

第四步:運作程式。

這一步可以寫一個專門的shell腳本,然後用qsub(一種SGE叢集任務投遞指令)進行任務投遞。作為一個單主機使用者,我就直接在指令行中運作。

source ~/opt/biosoft/falcon_unzip/fc_env/bin/activate
fc_run fc_run_ecoli_local.cfg &> fc_run.log &           

在程式運作過程中,可以通過幾行指令來看下程式的進度。

檢查不同階段總共要執行的任務數目

# raw-data
grep '^#' 0-rawreads/run_jobs.sh
# preads
grep '^#' 1-preads_ovl/run_jobs.sh           

目前已經執行的任務數

find 0-rawreads/ -name "job*done" | wc -l
find 0-rawreads/ -name "m_*done" | wc -l           

評估組裝運作結果

E. coli 的 subreads總共有1.01Gb資料,共105.451條reads。

我們可以用

dazzler

裡的指令行來确認dazzler資料庫是否建構正常

DBstats 0-rawreads/raw_reads.db | head -n 17           
使用Falcon對三代測序進行基因組組裝

DBstats

從中我們可以發現,隻有86.1%的read被用于建構dazzler資料庫,這是由于配置檔案裡我們設定的DBsplit參數使用長度大于500bp的read,

pa_DBsplit_option = -x500 -s50

此外,我們還需要檢查 預組裝(preassmebly)中實際用于組裝的資料量

cat 0-rawreads/report/pre_assembly_stats.json           
使用Falcon對三代測序進行基因組組裝

糾錯資料報告

這裡注意一點,這裡的

length_cutoff

是程式根據配置檔案中這一項

length_cutoff = -1

确定,這裡的"-`"意味着程式會自動根據基因組大小和覆寫度确定,但是這不一定是最佳選擇,你應該自己設定。

最後的資料存放在

2-asm-falcon/

檔案夾下,簡單用

ls -lh 2-asm-falcon/

檢視該檔案時,我發現組裝的

p_ctg.fa

大小就隻有1.4M, 而實際的大腸杆菌基因組大小為4.8M左右,差距有點大。

原因就是配置檔案中

length_cutoff

的設定問題,上圖顯示隻有22X用于組裝。

當我将其修改成

length_cutoff=2000

時,糾錯後還有156X資料用于組裝,最後的組裝的大小就成了4.8m。 是以用于糾錯的資料應該盡量多一點,但是門檻值也不要太低,否則增加了不必要的運算量。

由于大腸杆菌基因組小,1Gb的資料相當于深度為200X,深度非常的高,可以适當繼續提高

length_cutfoff

。當然大多數項目都是70X, 100X, 120X的資料, 設定成2000也就差不多了。

此外,建議多用幾個長度門檻值的preads進行組裝,即修改配置檔案中

length_cutoff_pr

, 然後對"2-asm-falcon", "mypwatcher",和 "all.log" 重命名, 這樣子重新運作腳本就會從糾錯後序列開始。

我嘗試了5k,10k,15k這三個參數,其中5k效果最差,裝出了20條序列,基因組隻有560k。10k裝出了3條序列,基因組大小是4.5M,N501.7M。 15k效果最好,裝出了一條4.6M的序列。

在length_cutoff分别是2k和15k情況下,length_cutoff_pr的長度都設定為15k,最後組裝結果是length_cutoff 設定為2k時,最後序列長度更長

官方文檔的坑

如果你在組裝過程中用

ps aux | grep 你的使用者名 | wc -l

檢查自己用了多少線程時。你會發現有一個時期會突然出現會出現好幾百個

python -m falcon_kit.mains.consensus

程序。

我在組裝自己物種的時候出現了700多個,最後256G記憶體完全被消耗完,導緻無法登陸伺服器隻好重新開機。

據我猜測應該是和

falcon_sense_option

中的

--n_cores

cns_concurrent_jobs

有關,我當時這兩個參數分别設定了20和32,那麼結果就要用到600多個程序。而每個程序大概會消耗1G記憶體,那麼峰值就會是600G。

而這裡組裝 E. coli 基因組這兩個參數分别是6和10,用

free -h

看記憶體消耗時也差不多時用了60多G記憶體。

是以一定要根據實際情況調整,配置檔案中和并行運算的參數,量力而行,不然重新開機等着你。

本文作者:徐洲更

繼續閱讀