🧬

Necatを用いた魚類ゲノムのde novo assembly

2023/08/30に公開

Tips

bioinformatics

ゲノムアセンブリ

tech

Necatとは？

PacBioと比較して少しエラー率の高いナノポアのロングリードを対象にしたアセンブラーです。

他のアセンブラーと比較して実行時間が短く、N50(コンティグ全体の配列長の50%以上となったコンティグ長)も長いといった利点があります。

原著論文では、大腸菌からシロイヌナズナやヒトまで検討しているようなので、Necatは幅広い生物に対応したアセンブラーと思っていいでしょう。

このツールが対象としてるナノポアのデータは、他のシーケンサーと比較して長鎖配列を生成できる一方、エラー分布が広く、エラー率の高いシーケンスが発生します。

そんな配列をNecatはまず、低エラー率サブシーケンス(LERS)の修正と、高エラー率サブシーケンス(HERS)の2段階に分けてエラー修正します(A-C)。

論文 Fig.2

次に、エラー修正されたcontigシーケンスを用いたアセンブル(D)を行い、ナノポアのraw readを用いてbridging(E)する事で連続性の高いcontig配列を生成することができるようです。

配列データの取得

prefetchとfasterq-dumpを使ってNCBI SRAから配列データを取得します。

この論文のデータを使用します。

Nanopore SRR12750269: 130GB程度のデータ量

Data	num seq	sum len	min len	ave len	max len	N50
SRR12750269.fastq.gz	4810089	130432770835	2	27116.5	884012	43600

mkdir 00sra 01fq
prefetch -p SRR12750269 -o 00sra/SRR12750269.sra && fasterq-dump -p -e 32 00sra/SRR12750269.sra -O 01fq/ && pigz 01fq/SRR12750269.fastq

上記コマンドでsraファイルがダウンロードされ、そのファイルからSRR12750269.fastq.gzファイルが生成されます。

仮想環境の作成

仮想環境を作成してactivateしてからツールをインストールすることで、ツールのバージョンを管理します。

python3が必要だったり、2が必要だったり、、これがベースとなる環境で起こると悲惨なので仮想環境は作って作業するに越したことはありません。

仮想環境といっても専用フォルダを作ってそこにツールをポイポイ入れていると思えばイメージはつきやすいです。

# Install pakcage and gen virtual env.
mamba create -n nanopore -y
# Acrivate
mamba activate nanopore

配列データの状態確認

NanoPlotはロングリードシーケンスデータのプロットツールです。webでの利用もできるようですが、今回はLocalで使用します。

https://github.com/wdecoster/NanoPlot

mambaを使ってインストール

mamba install -c bioconda NanoPlot -y

NanoPlotを実行します。

NanoPlot -t 8 --fastq SRR12750269.fastq.gz --loglength -o NanoPlotNormal

NanoPlotNormalフォルダにQscoteと配列長に関する以下のようなFigがpngとhtmlファイルで生成されます。

縦軸はクオリティスコア、横軸は配列長

見てみると、配列長が短くQscoreが低い配列が大半を占めている一方、35kbp付近の比較的長い配列はQscoreが10付近(10bpに1bpの確率でエラー塩基が含まれる)であることが分かります。

上記を踏まえてデータ処理を進めるわけですが、ナノポアの配列は開始点と終点のクオリティが低くなる傾向にあります。

それも別ツールで確認してみましょう。

https://github.com/wdecoster/nanoQC

mambaでインストールします。

mamba install -c bioconda nanoQC -y

QCを実行します。

nanoQC SRR12750269.fastq.gz --outdir SRR12750269.nanoqc

SRR12750269.nanoqcフォルダにnanoQC.htmlというファイルが出力されます。いくつかFigがありますが、Positionごとの平均Qscoreの情報は以下の図になります。

縦軸はQscoreで横軸は各末端からの配列長

一般的に言われているように前後半とも末端40bp程度までは、Qscoreが比較的低い値を示しています。

相対的に見て40bp程度までの配列は意味をなさないのでクオリティフィルタリングのステップで除去してしまいましょう。

Nanopore readのクオリティフィルタリング

chopper

chopperはNanoFiltとNanoLyseをRustという言語で再実装したフィルタリングとトリミングツールになります。

https://github.com/wdecoster/chopper

こちらもmambaでインストールします。

mamba install -c bioconda chopper -y

また配列データが大きいので、gz形式のファイルを並列解凍できるrapidgzipもインストールします。

https://github.com/mxmlnkn/rapidgzip

mamba install -c mxmlnkn rapidgzip -y

では次に、NanoPlotとnanoQCで確認した情報を元に配列長でフィルタリングを進めます。

両末端50bpをtrimします。以下のコマンドを実行すると00fqフォルダにTrim後の配列データであるSRR12750269.chopper.fastq.gzが出力されます。

rapidgzip -c -d -P 32 SRR12750269.fastq.gz | chopper --headcrop 50 --tailcrop 50 --threads 32 | pigz -c > 00fq/SRR12750269.chopper.fastq.gz

Necat実行のための設定

https://github.com/xiaochuanle/

mambaを使ってインストールします。

mamba install -c bioconda necat -y

もしNecatの実行でperl5libのエラーが出たら、自分の環境に合わせて自身の仮想環境のperl5にpathを通しておくと良いです。

こんな感じ↓

PERL5LIB=$HOME/miniconda3/envs/nanopore/lib/perl5/5.32/

USAGE

Usage: necat.pl correct|assemble|bridge|config cfg_fname
    correct:     correct rawreads
    assemble:    generate contigs
    bridge:      bridge contigs
    config:      generate default config file

解析パラメーターの設定

configファイルを作成して、解析対象データや生物種の情報について記載します。

config fileの作成

necat config Opo_config.txt

テンプレートファイルが作成されるので、以下のようにPROJECT=, ONT_READ_LIST=, GENOME_SIZE=, THREADS=の部分を編集します。

PROJECT=NecatRun
ONT_READ_LIST=FqList.txt
GENOME_SIZE=1156123904
THREADS=32

カレントディレクトリにFqList.txt作成してSRR12750269.chopper.fastq.gzの相対パスを記載する。

./00fq/SRR12750269.chopper.fastq.gz

ここまでこれば実行。

Necatの実行

Nanopore raw readのError修正

necat correct Opo_config.txt

パイプラインは40xの最長Raw read(PREP_OUTPUT_COVERAGE)を修正します。

修正されたReadは./NecatRun/1-consensus/cns_iter${NUM_ITER}/cns.fastaに出力され、30Xの最長修正済Read(CNS_OUTPUT_COVERAGE)が次のAssembleに利用されます。配列データは、./NecatRun/1-consensus/cns_final.fastaとして出力されています。

Error correctされたReadのAssembly

修正済みRaw readを以下のコマンドでアッセンブルします。
Error修正のステップが完了していない場合でも、コマンドは最初に修正のステップを自動的に実行します。

necat assemble Opo_config.txt

Assembly済みのcontigは./NecatRun/4-fsa/contigs.fastaとして出力されます。

Bridge contig

アセンブルされたコンティグを更にブリッジング(論文 Fig.2 E)します。

necat bridge Opo_config.txt

Bridging後のコンティグファイルは./NecatRun/6-bridge_contigs/bridged_contigs.fastaとして出力されます。

また、config.txtを明記した部分以外編集していなければ、POLISH_CONTIGSがtrueとなっているはずなので、同時に修正済みReadを使用してブリッジングコンティグをポリッシングします。

ポリッシュされたコンティグ配列は、./NecatRun/6-bridge_contigs/polished_contigs.fastaとして出力されます。

確認

seqkitを使ってプロファイルを確認します。

seqkit stats -T -a ./NecatRun/6-bridge_contigs/polished_contigs.fasta > necat_asm.tsv

tsvファイルが出力されます。抜粋すると以下のとおりです。

Tool	num seq	sum len	min len	ave len	max len	N50
Necat chopper	691	1,237,921,001	611	1,791,492	25,192,510	5,730,779

コンティグ配列は691本で、N50は5.6M baseでした。

ちなみにchopperによるトリミングを実施していないデータを並べてみると以下の通りです。

Tool	num seq	sum len	min len	ave len	max len	N50
Necat chopper	691	1,237,921,001	611	1,791,492	25,192,510	5,730,779
Necat	686	1,238,210,351	577	1,804,971.4	22,948,987	5,518,680

chopperによる処理を加えることで、N50のコンティグ長が約0.2Mbp、最大配列長も25Mbpと約2Mbp長くなりました。

Low Qualityなデータのtrimはアセンブリ結果の改善には有効なようです。

論文のN50は22.25Mbpだったので、Nanopore単独のアセンブリであればこんなもんなのかもしれませんが、もしかするとやり方間違ってるかもしれません。

ちなみにFlyeでアセンブルした結果をさっきの表に加えると以下のとおりです。

Tool	num seq	sum len	min len	ave len	max len	N50
Necat	686	1,238,210,351	577	1,804,971.4	22,948,987	5,518,680
Necat chopper	691	1,238,128,675	756	1,791,792.6	25,196,898	5,610,204
Flye	2,113	1,101,710,698	582	521,396	27,517,816	3,370,922

FlyeはNecatと同じポリッシング1回での実行ですが、--scaffoldオプションをつけました。

初めての実行だったので正しいかどうか自信ないです。bridge contig vs scaffoldなので同じ土台で話をするのはちょっと違うかな？

実行時間は圧倒的にNecatの方が短かったです。

Necat run time

real	745m12.084s
user	19680m55.369s
sys	69m42.640s

最後に今回生成された配列データファイルは大きいので、再帰的に探索して圧縮します。

find ./ -name "*.fasta" -type f -execdir pigz {} \;

次は追加でpolishingをかけてみます。

Information

自作PC1
necat : 0.0.1_update20200803
seqkit : 2.5.1
rapidgzip : 0.9.0
NanoPlot : 1.41.6
De Coster, W., & Rademakers, R. (2022). NanoPack2: population-scale evaluation of long-read sequencing data. Bioinformatics, 39.