Racon で配列を補正する

1. minimap2 と Racon を追加する

cd /home/bsjinzai
conda activate wgs-long
conda install --override-channels -c conda-forge -c bioconda \
  --strict-channel-priority minimap2 racon
minimap2 --version
racon --version

minimap2 で、各リードが組み立てた配列のどの場所に対応するかを調べます。Racon はその対応とリードを使い、合意配列を計算します。

2. リードをマッピングする

minimap2 -x map-ont -t 8 \
  assembly/nanopore/assembly.fasta \
  qc/nanopore.filtered.fastq.gz \
  > results/nanopore.paf 2> logs/minimap2.log
wc -l results/nanopore.paf

-x map-ont はNanoporeリードのマッピング用設定です。引数の順は対応付け先のアセンブリ、リードです。PAFに対応関係を保存します。行数はマッピングの記録数で、重複や複数の対応を含むためリード総数とは限りません。

3. Racon を1回実行する

racon -t 8 -u \
  qc/nanopore.filtered.fastq.gz \
  results/nanopore.paf \
  assembly/nanopore/assembly.fasta \
  > results/nanopore.racon.fasta 2> logs/racon.log

Racon の3つの入力は、リード、対応関係、補正対象の配列の順です。-u は、補正できなかった配列も出力に残す指定です。配列が出力に残っていても、その全体が補正されたとは限りません。

seqkit stats -a -T \
  assembly/nanopore/assembly.fasta \
  results/nanopore.racon.fasta \
  > results/racon.comparison.tsv
cat results/racon.comparison.tsv
conda list --explicit > logs/wgs-long.racon.explicit.txt

contig数や長さが変わったか記録します。長さが変わらなくても塩基の置換が起きている場合があります。長さやN50だけでは塩基精度の向上を判定できません。

4. 補正結果を確認する

Raconの出力はFASTAです。Flyeの assembly_graph.gfa は更新されないため、BandageではFlyeのグラフを観察します。

Flyeはすでにpolishingしているため、Raconを追加して必ず結果が良くなるとは限りません。必要なら後半の参照配列との比較で、不一致やindelの変化も調べます。さらにRaconを反復する場合は、更新したFASTAへ再度マッピングしてから実行します。

元のFASTAと results/nanopore.racon.fasta が両方残り、後者が空でなくSeqKitで読めることを確認します。Racon後のFASTAだけから環状化を再判定したことにはなりません。

参考：minimap2、Raconの入出力とオプション。