NECAT: Nanopore數(shù)據(jù)的高效組裝工具

NECAT是肖傳樂(lè)老師團(tuán)隊(duì)開(kāi)發(fā)的一個(gè)針對(duì)Nanopore數(shù)據(jù)組裝的軟件，目前該工具尚未發(fā)表，除了https://github.com/xiaochuanle/NECAT有軟件的介紹外，暫時(shí)沒(méi)有中文資料介紹NECAT的使用。

太長(zhǎng)不看的結(jié)論: Nanopore的組裝推薦用下NECAT。組裝之后是先用MEDAKA做一遍三代polish，然后用NextPolish默認(rèn)參數(shù)做二代polish。

這篇將會(huì)以一篇發(fā)表在Nature Communication上的擬南芥nanopore數(shù)據(jù)介紹如何使用NECAT進(jìn)行組裝，運(yùn)行在CentOS Linux release 7.3.1611 (Core)，64G為內(nèi)存， 20線程(Intel(R) Xeon(R) CPU E5-2640 v4 @ 2.40GHz)，下面是正文。

軟件安裝

NECAT可以在https://github.com/xiaochuanle/NECAT/releases/頁(yè)面獲取最新的軟件下載地址，這里下載的是0.01版本。

wget https://github.com/xiaochuanle/NECAT/releases/download/v0.01/necat_20190307_linux_amd64.tar.gz
tar xzvf necat_20190307_linux_amd64.tar.gz
export PATH=$PATH:$(pwd)/NECAT/Linux-amd64/bin

目前0.01版本不支持gz文件作為輸入，但后續(xù)版本應(yīng)該會(huì)支持。

實(shí)戰(zhàn)

第一步: 新建一個(gè)分析項(xiàng)目

mkdir NECAT && cd NECAT

以發(fā)表在NC上的擬南芥數(shù)據(jù)為例, 下載該數(shù)據(jù)

# 三代測(cè)序
wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR217/003/ERR2173373/ERR2173373.fastq.gz
seqkit seqkit fq2fa ERR2173373.fastq.gz | gzip -c > ERR2173373.fasta

第二步: 創(chuàng)建配置文件

necat.pl config ath_config.txt

配置文件中，主要修改如下幾個(gè)參數(shù)

PROJECT=athaliana #項(xiàng)目名
ONT_READ_LIST=read_list.txt #read所在路徑文件
GENOME_SIZE=120000000 #基因組大小
THREADS=20 # 線程數(shù)
MIN_READ_LENGTH=3000 # 最短的read長(zhǎng)度
CNS_OUTPUT_COVERAGE=45 # 用于組裝的深度

參數(shù)中還有一個(gè)，NUM_ITER=2，它并非是簡(jiǎn)單的重復(fù)2次糾錯(cuò)，它的每一輪的校正目的其實(shí)不同，第一輪的優(yōu)先級(jí)是敏感度(senstitive)， 第二輪之后主要追求速度(fast)。

除了上面的配置參數(shù)外，其他參數(shù)可以不需要修改，使用默認(rèn)的值即可。需要修改的話，參考最后的參數(shù)說(shuō)明部分。

第三步: 序列糾錯(cuò)

necat.pl correct ath_config.txt &

糾錯(cuò)后的reads在athaliana/1-consensus/cns_final.fasta

cns_finla.fasta的統(tǒng)計(jì)信息會(huì)輸出在屏幕中， 或者自己用fsa_rd_stat也能得到同樣的結(jié)果

Count: 206342
Tatal: 3102480870
Max: 112992
Min: 1010
N25: 31940
L25: 18989
N50: 21879
L50: 48506
N75: 13444
L75: 93215

此外我還用time獲取了運(yùn)行時(shí)間，糾錯(cuò)花了大概一個(gè)小時(shí)。

real    55m31.451s
user    815m32.801s
sys     7m55.039s

第四步: contig組裝

necat.pl assemble ath_config.txt &

結(jié)果在athaliana/4-fsa/contigs.fasta

關(guān)于contigs.fata統(tǒng)計(jì)信息會(huì)輸出在屏幕上，同樣用fsa_rd_stat 也可以。

Count: 162
Tatal: 122293198
Max: 14562810
Min: 1214
N25: 13052494
L25: 3
N50: 9503368
L50: 5
N75: 4919866
L75: 10

時(shí)間用了75分鐘

real    74m53.127s
user    1308m29.534s
sys     12m5.032s

第五步: contig搭橋

necat.pl bridge ath_config.txt

結(jié)果在athaliana/6-bridge_contigs/bridged_contigs.fasta

Count: 127
Tatal: 121978724
Max: 14562810
Min: 2217
N25: 13193939
L25: 3
N50: 11146374
L50: 5
N75: 5690371
L75: 9

從N50和N75可以看出這一步會(huì)提高組裝的連續(xù)性。

組裝結(jié)果polish

對(duì)Nanopore組裝結(jié)果進(jìn)行polish的常用軟件有下面3個(gè)

由于擬南芥的基因組比較小，我分別用了Medaka和racon對(duì)輸出結(jié)果進(jìn)行polish(因?yàn)闆](méi)有原始信號(hào)數(shù)據(jù)，因此nanopolish用不了)，代碼如下

Medaka

NPROC=20
BASECALLS=ERR2173373.fasta
DRAFT=athaliana/6-bridge_contigs/bridged_contigs.fasta
OUTDIR=medaka_consensus
medaka_consensus -i ${BASECALLS} -d ${DRAFT} -o ${OUTDIR} -t ${NPROC} -m r941_min_high

參數(shù)目前發(fā)生了變化，具體Medaka的使用參考https://github.com/nanoporetech/medaka

三輪Racon:

gzip -dc ERR2173373.fastq.gz  > ERR2173373.fastq
minimap2 -t 20 ${DRAFT} ERR2173373.fastq > round_1.paf 
racon -t 20 ERR2173373.fastq round_1.paf ${DRAFT} > racon_round1.fasta
minimap2 -t 20 racon_round1.fasta ERR2173373.fastq > round_2.paf 
racon -t 20 ERR2173373.fastq round_2.paf racon_round1.fasta> racon_round2.fasta
minimap2 -t 20 racon_round2.fasta ERR2173373.fastq > round_3.paf
racon -t 20 ERR2173373.fastq round_3.paf racon_round2.fasta> racon_round3.fasta

在后續(xù)評(píng)估質(zhì)量的時(shí)候，我發(fā)現(xiàn)單純用三代polish的結(jié)果還不是很好，因此我用他們提供的二代測(cè)序，用NextPolish對(duì)NECAT的結(jié)果進(jìn)行polish。

# 二代測(cè)序
prefetch ERR2173372
fasterq-dump -O  . ERR2173372

run.cfg內(nèi)容如下, 其中sgs.fofn記錄的就是解壓后的ERR2173372_1.fastq和ERR2173372_2.fastq的路徑

[General]                
job_type = local  
job_prefix = nextPolish  
task = 1212 
rewrite = no 
rerun = 3
parallel_jobs =  8  
multithread_jobs = 20 
genome = input.fasta 
genome_size = auto 
workdir = ./nextpolish 
polish_options = -p {multithread_jobs}
[sgs_option]             
sgs_fofn = ./sgs.fofn 
sgs_options = -max_depth 100 -bwa

我考慮了兩種情況，一種是直接用二代polish，另一種是三代polish之后接二代polish。

結(jié)果評(píng)估

在計(jì)算時(shí)間上，我之前用Canu跑了相同的數(shù)據(jù)，設(shè)置原始錯(cuò)誤率0.5，糾錯(cuò)后錯(cuò)誤率為0.144，用3個(gè)節(jié)點(diǎn)（每個(gè)節(jié)點(diǎn)12個(gè)線程），運(yùn)行了3天時(shí)間，但是NECAT只需要3個(gè)小時(shí)左右就能完成相同的分析，這個(gè)速度差異實(shí)在是太明顯了。

用Minimap2 + dotPlotly繪制CANU，NECAT和擬南芥參考基因組的共線性圖

minimap2 -t 20 -x asm5 Athaliana.fa NECAT.fa > NECAT.paf
pafCoordsDotPlotly.R  -i NECAT.paf -o NECAT  -l -p 10 -k 5
minimap2 -t 20 -x asm5 Athaliana.fa CANU.fa > CANU.paf
pafCoordsDotPlotly.R  -i CANU.paf -o CANU  -l -p 10 -k 5

NECAT的結(jié)果

NECAT

CANU的結(jié)果

CANU

NECAT和CANU都和參考基因組有著良好的共線性，但是NECAT的連續(xù)性更好，幾乎成一條直線。

之后，我使用了QUAST來(lái)評(píng)估Canu，NECAT初步組裝，NECAT用Medaka， nanopolish和racon糾錯(cuò)的結(jié)果（MD: MEDAKA, RC: RACON, NP:NextPolish）。

quast.py -t 100 --output-dir athaliana --circos \
    CANU.fa \
    NECAT.fa \
    NECAT_MD.fa \
    NECAT_MD_NP.fa \
    NECAT_NP.fa \
    NECAT_RC.fa \
    NECAT_RC_NP.fa \
    -r Athaliana.fa  \
    -g TAIR10_GFF3_genes.gff &

一些描述基本信息

CANU        N50 = 4875070,  L50 = 7, Total length = 114689024, GC % = 36.09 
NECAT       N50 = 11146374, L50 = 5, Total length = 121978724, GC % = 36.50
NECAT_MD    N50 = 11216803, L50 = 5, Total length = 122101599, GC % = 36.54
NECAT_MD_NP N50 = 11405151, L50 = 5, Total length = 124142955, GC % = 36.30
NECAT_NP    N50 = 11399084, L50 = 5, Total length = 124735066, GC % = 36.36
NECAT_RC    N50 = 11212098, L50 = 5, Total length = 122519370, GC % = 36.4
NECAT_RC_NP N50 = 11406553, L50 = 5, Total length = 124618502, GC % = 36.34

在BUSCO完整度上， 以embryophyta_odb10作為物種數(shù)據(jù)庫(kù), 其中ONTmin_IT4是發(fā)表的文章里的結(jié)果, Athalina則是擬南芥的參考基因組，我們以它們的BUSCO值作為參照。

Athalina     : C:98.6%[S:98.0%,D:0.6%],F:0.4%, M:1.0%, n:1375
ONTmin_IT4   : C:98.4%[S:97.7%,D:0.7%],F:0.7%, M:0.9%, n:1375
CANU         : C:22.9%[S:22.8%,D:0.1%],F:20.2%,M:56.9%,n:1375
NECAT        : C:36.6%[S:36.6%,D:0.0%],F:22.9%,M:40.5%,n:1375
NECAT_MEDAKA : C:53.6%[S:53.2%,D:0.4%],F:21.0%,M:25.4%,n:1375
NECAT_RACON  : C:45.3%[S:45.2%,D:0.1%],F:23.1%,M:31.6%,n:1375

二代Polish后的BUSCO結(jié)果如下(MD: MEDAKA, RC: RACON, NP:NextPolish)：

Athalina   : C:98.6%[S:98.0%,D:0.6%],F:0.4%,M:1.0%,n:1375
ONTmin_IT4 : C:98.4%[S:97.7%,D:0.7%],F:0.7%,M:0.9%,n:1375
NECAT_NP   : C:98.6%[S:97.9%,D:0.7%],F:0.4%,M:1.0%,n:1375   
NECAT_MD_NP: C:98.7%[S:98.0%,D:0.7%],F:0.4%,M:0.9%,n:1375
NECAT_RC_NP: C:98.5%[S:97.8%,D:0.7%],F:0.4%,M:1.1%,n:1375

從以上這些數(shù)據(jù)，你可以得到以下幾個(gè)洞見(jiàn):

  • 在Nanopore的組裝上，NECAT效果優(yōu)于Canu，無(wú)論是連續(xù)性還是N50上
  • MEDAKA三代polish效果好于RACON。在速度上，MEDAKA比三遍RACON都慢，并且MEDAKA會(huì)將一些可能的錯(cuò)誤組裝給打斷
  • Nanopore的數(shù)據(jù)用NECAT組裝后似乎用NextPolish進(jìn)行polish后就行，但是由于物種比較小，可能不具有代表性。

結(jié)論: Nanopore的組裝建議用NECAT。組裝之后是先用MEDAKA做一遍三代polish，然后用NextPolish默認(rèn)參數(shù)做二代polish。

配置文件補(bǔ)充

這部分對(duì)配置文件做一點(diǎn)簡(jiǎn)單補(bǔ)充。

下面這些參數(shù)相對(duì)簡(jiǎn)單，不需要過(guò)多解釋，按照自己需求修改

  • CLEANUP: 運(yùn)行完是否清理臨時(shí)文件，默認(rèn)是0，表示不清理
  • USE_GRID: 是否使用多節(jié)點(diǎn), 默認(rèn)是false
  • GRID_NODE: 使用多少個(gè)節(jié)點(diǎn)，默認(rèn)是0，當(dāng)USE_GRID為true時(shí)，按照自己實(shí)際情況設(shè)置

以下的參數(shù)則是需要根據(jù)到具體的軟件中去查看具體含義，需要和軟件開(kāi)發(fā)者討論

  • OVLP_FAST_OPTIONS: 第二輪糾錯(cuò)時(shí), 傳給oc2pmov
  • OVLP_SENSITIVE_OPTIONS: 第一輪糾錯(cuò)時(shí), 傳給oc2pmov
  • CNS_FAST_OPTIONS: 第二輪糾錯(cuò)時(shí)，傳給oc2cns
  • CNS_SENSITIVE_OPTIONS: 第一輪糾錯(cuò)時(shí)，傳給oc2cns
  • TRIM_OVLP_OPTIONS: 傳給oc2asmpm
  • ASM_OVLP_OPTIONS: 傳給oc2asmpm
  • FSA_OL_FILTER_OPTIONS: 參數(shù)傳給fsa_ol_filter
  • FSA_ASSEMBLE_OPTIONS: 參數(shù)傳給fsa_assemble
  • FSA_CTG_BRIDGE_OPTIONS: 參數(shù)傳給fsa_ctg_bridge

參考資料

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容