主頁:github: PacificBiosciences/FALCON
簡介
Falcon是一組通過快速比對長reads,從而來consensus和組裝的工具。
Falcon工具包是一組簡單的代碼集合,我使用它們來研究單倍體和二倍體基因組的高效組裝算法。
為了提高計算速度,它有一些后台代碼是使用C來實現的,為了方便一些簡單的前端是用Python編寫的。
Falcon不是一個傻瓜的組裝工具(除了很小的基因組),為了得到最好的結果,你可能需要了解各種分布式計算系統和一些基本的基因組組裝理論。FAQ頁面有一些常見的問題及其回答。
Falcon 基因組組裝工具包 - 使用指南
1. 分層次基因組組裝過程的概述
主要有如下幾個步驟:
- Raw sub-reads overlapping for error correction
- Pre-assembly and error correction
- Overlapping detection of the error corrected reads
- Overlap filtering
- Constructing graph from overlaps
- Constructing contig from graph
簡單翻譯:
- 使用Raw sub-reads 構建重疊,准備進行錯誤校正
- 預組裝和錯誤校正
- 錯誤校正后的reads的重疊檢測
- 重疊過濾
- 用重疊來構造圖
- 用圖來構造contig
每個步驟使用不同的命令行工具,實現不同的算法來完成這項工作。同樣,每一步的計算需求也大不一樣。假定用戶已經擁有合理數量的計算資源。
例如,在合理的時間內裝配100M大小的基因組,可能需要至少32核cpu 和 128 Gb RAM。代碼是按照集群計算環境來編寫的。需要一個作業隊列來進行長時間的腳本工作,以及cpu-rich的計算工作隊列。
一個包含了一組sub-reads的文件,fc_run.py
腳本可以驅動工作流來管理和檢查數據依賴和提交每一步的工作,從給定的數據中生成一個組裝的草圖。
fc_run.py
是工作流驅動腳本,需要運行在一台計算機上,允許持續長時間地工作一段時間來完成整個裝配過程。它需要一個配置文件fc_local.cfg作為單個輸入。原始序列數據的輸入文件包含在配置文件中。
配置文件可用於控制計算資源的使用和重要算法參數,根據輸入數據集達到最有組裝的效果。然而,沒有萬能的方式來猜測各種具體項目的最合適的選項和參數。參數調優的最好方法是了解一些組裝理論,一點點實現方法,這樣才能正確的了解更改某些選項對組裝結果的影響。快速瀏覽reads長度分布、整體覆蓋和調整相應的選項也是非常重要。
這里我們會詳細介紹基因組分層組裝策略 以及 fc_run.py
各個參數的含義
快速入門
2. 使用Raw sub-reads構建重疊(為了錯誤校正)
這個版本的 Falcon 的 overlapping 是用 Gene Myers' Daligner
的修改版本來做的(參考 github 的 forked version)。主要的修改是adapting some I/O for the downstream bioinformatics process。可以通過git diff
來查看有哪些修改。
git diff #查看版本之間差異
input_fofn
選項指向了包含所有輸入數據的文件。來自Daligner
的fasta2DB
在 fc_run.py
會被調用(I/O集中器,會在你執行 fc_run.py
的電腦節點出運行,如果這在你的服務器集群中是個問題,你必須修改代碼來打包相關的操作到一個腳本中,以便提交到你的任務管理系統中)。
input_fofn #指向了包含所有輸入數據的文件
這個版本的 fc_run.py
支持從錯誤校正的reads中運行組裝程序(可以跳過錯誤校正環節,直接開始組裝)
input_type = preads #輸入為 all error-corrected reads,會跳過錯誤校正,直接開始最終的重疊群組裝步驟
input_type = raw #原始數據
你需要決定cutoff的長度,一般,選擇一定的閾值,使得你能得到15x 到 20x的覆蓋度是比較好的。 然而,如果計算資源充足,校正reads你有其他用途,你可以設置lower length cutoff,來獲得更多的error corrected reads。
最終組裝的時候,更多的reads不意味着更好的組裝結果,因為會有噪音。
有時,嘗試不同的length_cutoff_pr
是有意義的,較初次overlapping step for error correction而言,后期計算更cheap。
我們的策略是選擇一個短的length_cutoff 進行一次計算。之后, 我們使用不同的length_cutoff_pr來獲得更好的組裝結果。
length_cutoff #控制錯誤校正過程中的cutoff
length_cutoff_pr #控制之后組裝重疊群步驟的cutoff
pa_concurrent_jobs
選項控制由fc_run.py
提交的 number of concurrent jobs。sge_option_da
和 sge_option_la
控制作業隊列和daligner jobs 的 number of slots。daligner
使用的線程數默認是4。然而,依據集群設置、計算節點的內存數量,你可以使用超過4的slots。為了選擇最好的數量,你最好咨詢你們實驗室的HPC專家,做一些前期的測試。
pa_concurrent_jobs #fc_run.py提交的當前作業數
sge_option_da #作業隊列
sge_option_la #daligner作業的slots數
最終運行作業的數量取決於你怎么"splits"序列數據庫,仔細閱讀Gene Myers's blog ( http://dazzlerblog.wordpress.com ) ,了解如何設置 pa_DBsplit_option
和 pa_HPCdaligner_option
選項。一般,對大型基因組而言,在pa_DBsplit_option
中,你應該使用-s400
(400Mb sequence per block)。這樣作業的數量會變少,但是持續時間會延長。如果,作業調度系統限制了一個作業可以運行的時間,你應該設置a smaller number for the -s
option。
另一個影響作業總數的參數是pa_HPCdaligner_option
中的 -dal
option。它決定how many blocks are compared to each in single jobs。大的 -dal
給出了大的作業,但是數量少;反之。小的,作業小,但是提交的數量多。
這個工作流中,daligner
生成的trace point沒有被用到(為了高效,應該使用trace point,但是首先你的知道如何正確的pull out它們)。pa_HPCdaligner_option
中的 -s1000
使得trace point變稀疏,以此節約磁盤空間。我們使用-l1000
忽視短於1kb的reads。
pa_DBsplit_option #-s400 for large genomes smaller number for the -s
pa_HPCdaligner_option #
-dal # determines how many blocks are compared to each in single jobs;
-s1000 #
-l1000 #
3. 預組裝和錯誤校正
daligner
的輸出時一組.las
文件,它包含了reads間相互比對的信息。Such information is dumped as sequences for error correction by a binary executable LA4Falcon to fc_consensus.py.(什么意思?) fc_consensus.py
生成了consensus(這部分是用C寫的)。
fc_consensus.py
有許多選項,可以使用falcon_sense_option
來控制它。大部分情況下,--min_cov
和 --max_n_read
是最重要的選項。--min_cov
控制了什么時候a seed read gets trimmed or broken due to low coverage(什么意思?),--max_n_read
puts a cap on the number of reads used for error correction(什么意思?)。
在高度重復的基因組匯總,設置小的 --max_n_read
,來確保 consensus code does not waste time aligning repeats。最長的合適的overlaps 被用來進行correction ,來減少the probability of collapsed repeats。
你可以使用cns_concurrent_jobs
來設置提交給任務管理系統的當前作業的最大數量。
falcon_sense_option #
--min_cov #when a seed read gets trimmed or broken due to low coverage???
--max_n_read #
cns_concurrent_jobs #
4. 錯誤校正后的reads的overlaps檢測
這部分與最初的overlapping 十分相似,只有部分修改,daligner
只以本地的raw reads作為默認值。
fc_run.py generates a fasta file of error-corrected reads where the fasta header is parse-able by daligner.(什么意思?)
下面的參數控制這個步驟的計算過程:
sge_option_pda = -pe smp 8 -q jobqueue
sge_option_pla = -pe smp 2 -q jobqueue
ovlp_concurrent_jobs = 32
ovlp_DBsplit_option = -x500 -s50
ovlp_HPCdaligner_option = -v -dal4 -t32 -h60 -e.96 -l500 -s1000
這個設置與第一步的overlapping 相平行,主要的區別是ovlp_HPCdaligner_option
的-e
選項。錯誤率小了很多,因此我們期待 much higher correlation between the p-reads。
5. overlaps過濾
不是所有的overlaps 都是獨立的,所以利用某些filtering step來減少計算和組裝圖的復雜性是可能的。例如:如果一個read被全部包含在另一個read里,那么它們的overlap不會提供額外的構建基因組的信息。由於重疊關系的傳遞屬性(transitive property),許多overlap 信息可以被簡單的推測出來。
事實上,構建contigs 的第一步就是刪除 transitive reducible edges(???),這意味着,我們只需要best n overlaps
在5' or 3' 端。overlap_filtering_setting
option的--bestn
parameter參數可以被用來控制每個read報告的maximum overlap。
overlap_filtering_setting
--bestn #control the maximum overlap reported for each read
另一個有用的啟發式,僅保留有平均 5' and 3' 覆蓋的reads。這是因為,如果一個read以重復結束,它會有更高的重復。這種reads不提供價值,不能解決相關的重復問題。我們可以過濾掉它們,希望能有reads(跨越重復,在兩端有正常的獨特的anchors)。同樣,如果一個reads一端的覆蓋度太低,那它可能有太多的測序錯誤。這種 reads產生了 "spurs"在組裝圖中,通常都會將它們過濾掉。--max_cov
和 --min_cov
就是用來過濾那些有太高或者太低的reads。
--max_cov #過濾太高的reads
--min_cov #過濾太低的reads
過濾的腳本也允許過濾一些 "split" reads。如果一個read在5' 和 3'端有非常不相等的覆蓋度,這也是某一段是重復的信號。 --max_diff
可以用來過濾掉一端比另一端的覆蓋哦高很多的reads。
--max_diff #過濾掉一端比另一端的覆蓋哦高很多的reads
使用這些參數的正確的數量是多少?真的很難設置正確,如果 error corrected reads的總覆蓋度比知道的length cut off 並且合理的高 (e.g. greater than 20x),那么將平均的覆蓋度min_cov
設置為5,max_cov
設置為3倍將會很安全, max_diff
設置為2倍。然而,在低覆蓋度的情況下,將min_cov
設為1或2將會更好。一個名為fc_ovlp_stats.py
的幫助腳本可以幫助dump(丟棄)。 A helper script called fc_ovlp_stats.py can help to dump the number of the 3' and 5' overlap of a given length cutoff, you can plot the distribution of the number of overlaps to make a better decision.(???)
你也可以設置max_diff
和min_cov
為很高的值,來避免過濾,在某些特殊的場合。
這個過濾的過程會顯著過濾掉高度拷貝重復的信息。也就是,這些重復會被過濾掉,並且不會在最終的組裝結果中顯示。如果你對這些重復感興趣,總之,單獨處理這些重復會更加有效和實用。如何解決重復序列這個問題具有非常重要的意義。
6. 用overlaps構造圖
鑒於overlapping 數據,string graph 由 fc_ovlp_to_graph.py
來創建,使用默認的參數。 fc_ovlp_to_graph.py
生成一些文件代表組裝后的最終的 string graph。最終的ctg_path
包含了每一個contig的圖的信息。一個 contig是一個simple paths 和 compound paths的線性的path。Compound paths are those subgraph that is not simple but have unique inlet and outlet after graph edge reduction.
They can be induced by genome polymorphism or sequence errors. By explicitly encoding such information in the graph output, we can examine the sequences again to classify them later.
7. 用圖來構造contig
create draft contigs的最后一步就是find a single path of each contig graph,以及generate sequences accordingly。在a contig graph is not a simple path的情況下,我們尋找 end-to-end path that has the most overlapped bases。這被稱為 primary contig
。圖中的每一個path,如果它與 primary one
不同,我們會construct the associated contig。在the associated contigs are induced by sequencing error的情況下,the identity of the alternative contig and the primary contig will be high ( > 99% identity most of time)。在there are true structure polymorphism的情況下,there are typically bigger difference between the associated contigs and the primary contigs。
fc_graph_to_contig.py
腳本,takes the sequence data and graph output to construct contigs。它會生成所有的 associated contigs at this moment。將來,一些后處理程序,來刪除errors中的some of the associated contigs將會被開發出來。(正在開發)
8. 通用的daligner選項
daligner由pa_HPCdaligner_option
和ovlp_HPCdaligner_option
控制。
為了限制內存,我們可以使用-M
選項。對於人類基因組組裝,我們測試 -M 32,每個daligner使用32G RAM。其他可能性正在調查中。
更多的daligner的細節,參見: dazzlerblog。
工作目錄結構、工作恢復和故障排除
代碼設計工作在單一目錄。工作目錄的典型布局是這樣的:
$ ls -l
total 56
drwxr-xr-x 84 jchin Domain Users 8192 Nov 30 12:30 0-rawreads
drwxr-xr-x 18 jchin Domain Users 4096 Nov 30 12:33 1-preads_ovl
drwxr-xr-x 2 jchin Domain Users 4096 Nov 30 12:44 2-asm-falcon
-rwxr-xr-x 1 jchin Domain Users 1041 Nov 30 12:13 fc_run.cfg
-rw-r--r-- 1 jchin Domain Users 378 Nov 29 23:20 input.fofn
drwxr-xr-x 2 jchin Domain Users 4096 Nov 30 12:13 scripts
drwxr-xr-x 2 jchin Domain Users 24576 Nov 30 12:33 sge_log
典型的input.fofn
是這樣的:
/mnt/data/Analysis_Results/m140913_050931_42139_c100713652400000001823152404301535_s1_p0.1.subreads.fasta
/mnt/data/Analysis_Results/m140913_050931_42139_c100713652400000001823152404301535_s1_p0.2.subreads.fasta
/mnt/data/Analysis_Results/m140913_050931_42139_c100713652400000001823152404301535_s1_p0.3.subreads.fasta
1. 內部0-rawreads目錄
目錄0-rawreads
包括:overlapping原始序列的所有的腳本和數據,它包含不同的job_ *
和m_ *
目錄。
.
例如,如果我們將E. coli的數據分成20個塊,目錄就會如下:
cns_done job_00011 job_00024 job_00037 job_00050 m_00003 m_00016
da_done job_00012 job_00025 job_00038 job_00051 m_00004 m_00017
job_00000 job_00013 job_00026 job_00039 job_00052 m_00005 m_00018
job_00001 job_00014 job_00027 job_00040 job_00053 m_00006 m_00019
job_00002 job_00015 job_00028 job_00041 job_00054 m_00007 m_00020
job_00003 job_00016 job_00029 job_00042 job_00055 m_00008 preads
job_00004 job_00017 job_00030 job_00043 job_00056 m_00009 prepare_db.sh
job_00005 job_00018 job_00031 job_00044 job_00057 m_00010 raw_reads.db
job_00006 job_00019 job_00032 job_00045 job_00058 m_00011 rdb_build_done
job_00007 job_00020 job_00033 job_00046 job_00059 m_00012 run_jobs.sh
job_00008 job_00021 job_00034 job_00047 las_files m_00013
job_00009 job_00022 job_00035 job_00048 m_00001 m_00014
job_00010 job_00023 job_00036 job_00049 m_00002 m_00015
job_ *
目錄存儲了每個daligner的輸出工作。
m_ *
目錄是用於合並工作。
有一些空文件,它們的名字是done
結尾。這些文件的時間戳是用來跟蹤工作流階段。您可以修改時間戳和重啟fc_run.py
來觸發為工作流程的一部分做某些做計算。然而,不推薦這樣,除非你有一些時間來讀fc_run.py
的源代碼,並且知道如何跟蹤工作流內部的依賴性。(例如,如果你在成功組裝后touch rdb_build_done
,並且重新運行fc_run.py
。因為所有中間過程取決於文件,rdb_build_done
比任何中間文件的都要新,就會觸發fc_run.py
,重復整個裝配過程。
. las_files
存儲了最后比對信息。
如果你不打算重新運行工作,但是想知道如何比對的樣子,你可以刪除所有job_ *
和m_ *
目錄但保留las_files
和preads
目錄。
.
這一步的主要輸出存儲在0-rawreads /preads
。0-rawreads/preads
內部的out.%04d.fa
的是輸出的reads的fasta文件。如果這一步成功完成, 哨兵文件cns_done
將被創建。
2. 1-preads_ovl目錄
該目錄存儲着p-read與p-read overlapping的數據。它與0-rawreads大體一致,但是卻沒有consensus 這一步。主要輸出是1-preads_ovl/
目錄里的 *.las
。
3. 2-asm-falcon目錄
這是最終的輸出目錄。
它包含了組裝圖和draft contigs的所有信息,細節描述參照Graph output format
一節。
fc_ovlp_to_graph.py 選項 和 組裝圖輸出格式
下面是fc_ovlp_to_graph.py
的使用信息:
usage: fc_ovlp_to_graph.py [-h] [--min_len MIN_LEN] [--min_idt MIN_IDT]
[--lfc]
overlap_file
a example string graph assembler that is desinged for handling diploid genomes
positional arguments:
overlap_file a file that contains the overlap information.
optional arguments:
-h, --help show this help message and exit
--min_len MIN_LEN minimum length of the reads to be considered for
assembling
--min_idt MIN_IDT minimum alignment identity of the reads to be considered
for assembling
--lfc use local flow constraint method rather than best overlap
method to resolve knots in string graph
De-dup example and strategy
TODO
Low coverage assembly
TODO
Assembly layout rule
CPU使用
Applications of the assembly graph
TODO
Reproducibility and replicability
序列 alignment 和 consensus 的C代碼
Other TODOs
Incremental overlapping
Pre-processing repeat for overlapping
術語
subread:子read,
full-pass subread:
pre-assembly:
error correction:
proper overlap:
string graph:
contig:
primary contig:
associated contig:
p-read:
compound path:
simple path:
Quiver:
科普時間:
什么是virtualenv?
每個應用可能需要各自擁有一套“獨立”的Python運行環境。virtualenv就是用來為一個應用創建一套“隔離”的Python運行環境。
資料:
virtualenv - 廖雪峰
virtualenv -- python虛擬沙盒
Virtualenv入門基礎教程
幫助文檔
現在默認分支是“master”,其包含Falcon最新的源碼。
當前最新的綜合版本是v0.3.0,安裝指南參照:v0.3+ Integration Installation Guide
上一個版本是v0.2.2,具體參照:v0.2.2 release github repository