QUAST でアセンブリを比べる

1. 比較対象を整理する

2. QUAST を実行する

cd /home/bsjinzai
conda create -n wgs-eval \
  --override-channels -c conda-forge -c bioconda \
  --strict-channel-priority quast
conda activate wgs-eval
quast.py --version

quast.py \
  assembly/nanopore/assembly.fasta \
  results/nanopore.racon.fasta \
  assembly/illumina/contigs.fasta \
  assembly/hybrid/assembly.fasta \
  -l Flye,Flye_Racon,Illumina,Hybrid \
  --min-contig 500 \
  -t 8 \
  -o results/quast \
  > logs/quast.stdout.log 2> logs/quast.stderr.log
conda list --explicit > logs/wgs-eval.explicit.txt

--min-contig 500 で、基本指標の対象を500 bp以上に揃えます。これより短い配列は主な比較から外れるため、短い配列の回収状況は別途確認します。レポートには別の長さ閾値の行もあるので、同じ行どうしを比較してください。

3. レポートを読む

Windows側のPowerShellで、QUASTの結果をフォルダごとダウンロードします。

New-Item -ItemType Directory -Force "$HOME\genome-workshop\results"
scp -P PORT -r USER@HOST:/home/bsjinzai/results/quast "$HOME\genome-workshop\results"

-rでレポートに必要な関連ファイルもコピーされます。Windowsのブラウザーでgenome-workshop\results\quast\report.htmlを開きます。表形式の記録にはreport.tsvを使えます。

N50はつながりの指標です。大きくなっても塩基が正しくなったとは言えません。Racon後の変化が小さい場合も、その値だけで補正の有無や効果を結論づけないでください。

4. グラフと並べて考える

比較の記録シートへ値と観察を転記します。環状化の情報はFlyeやUnicyclerの出力から記録し、QUASTのN50から推定しません。Racon後については、Flyeの構造情報を参照していることを明記します。

• Illumina単独とNanopore単独では、連続性とグラフにどんな違いがあったか。
• Hybridで変化した配列は、長い配列か短い配列か。
• 総長が増減した理由として、どんな可能性があるか。
• いまの結果から分かることと、追加の検証が必要なことは何か。

5. 発展：菌種の開示後に参照と比べる

菌種の開示後に行う追加課題です。参照配列はdata/reference.fastaを使います。対応する染色体と付随配列を含む完成ゲノムです。

quast.py \
  assembly/nanopore/assembly.fasta \
  results/nanopore.racon.fasta \
  assembly/illumina/contigs.fasta \
  assembly/hybrid/assembly.fasta \
  -l Flye,Flye_Racon,Illumina,Hybrid \
  -r data/reference.fasta \
  --min-contig 500 \
  -t 8 \
  -o results/quast-reference \
  > logs/quast-reference.stdout.log 2> logs/quast-reference.stderr.log

参照ありではgenome fraction、mismatch、indel、misassemblyなどを確認できます。参照なしの結果からこれらを推測しません。参照との差にはアセンブリエラーだけでなく、実際の試料差や参照側の問題が含まれる可能性もあります。

DFASTへ渡すFASTAを選び、採用理由を記録します。次章はHybridのFASTAを使う例です。

参考：QUASTマニュアル。