Falcon:三代reads比對組裝工具箱


Falcon
主頁:github: PacificBiosciences/FALCON

簡介

Falcon是一組通過快速比對長reads,從而來consensus和組裝的工具。

Falcon工具包是一組簡單的代碼集合,我使用它們來研究單倍體和二倍體基因組的高效組裝算法。

為了提高計算速度,它有一些后台代碼是使用C來實現的,為了方便一些簡單的前端是用Python編寫的。

Falcon不是一個傻瓜的組裝工具(除了很小的基因組),為了得到最好的結果,你可能需要了解各種分布式計算系統和一些基本的基因組組裝理論。FAQ頁面有一些常見的問題及其回答。


Falcon 基因組組裝工具包 - 使用指南

1. 分層次基因組組裝過程的概述

主要有如下幾個步驟:

  1. Raw sub-reads overlapping for error correction
  2. Pre-assembly and error correction
  3. Overlapping detection of the error corrected reads
  4. Overlap filtering
  5. Constructing graph from overlaps
  6. Constructing contig from graph

簡單翻譯:

  1. 使用Raw sub-reads 構建重疊,准備進行錯誤校正
  2. 預組裝和錯誤校正
  3. 錯誤校正后的reads的重疊檢測
  4. 重疊過濾
  5. 用重疊來構造圖
  6. 用圖來構造contig

每個步驟使用不同的命令行工具,實現不同的算法來完成這項工作。同樣,每一步的計算需求也大不一樣。假定用戶已經擁有合理數量的計算資源。

例如,在合理的時間內裝配100M大小的基因組,可能需要至少32核cpu128 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 選項指向了包含所有輸入數據的文件。來自Dalignerfasta2DBfc_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_dasge_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_optionpa_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_diffmin_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選項

dalignerpa_HPCdaligner_optionovlp_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_filespreads 目錄。
.
這一步的主要輸出存儲在0-rawreads /preads0-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


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM