SPAdes で Illumina リードを組み立てる

1. 後半のデータを確認する

Linuxの/home/bsjinzai/data/にある、試料Aのillumina_1.fastq.gzとillumina_2.fastq.gzを使います。R1とR2は同じDNA断片の両端を読んだペアです。片方だけ並べ替えたり、別の試料と組み合わせたりしません。

cd /home/bsjinzai
ls -lh data/illumina_1.fastq.gz data/illumina_2.fastq.gz
gzip -t data/illumina_1.fastq.gz data/illumina_2.fastq.gz
conda activate wgs-long
seqkit stats -a -T data/illumina_1.fastq.gz data/illumina_2.fastq.gz \
  > results/illumina.raw.stats.tsv
cat results/illumina.raw.stats.tsv

Nanoporeとのリード長の違いと、R1・R2のリード数を確認します。件数が同じでもペア対応が正しいとは限りません。

2. fastp で前処理する

conda create -n wgs-short \
  --override-channels -c conda-forge -c bioconda \
  --strict-channel-priority fastp
conda activate wgs-short
fastp --version

fastp \
  -i data/illumina_1.fastq.gz \
  -I data/illumina_2.fastq.gz \
  -o qc/illumina_1.clean.fastq.gz \
  -O qc/illumina_2.clean.fastq.gz \
  --detect_adapter_for_pe \
  --thread 8 \
  --html qc/fastp.html \
  --json qc/fastp.json \
  2> logs/fastp.log

-i と -I がペアの入力、-o と -O が出力です。アダプター検出と品質フィルタを行います。処理後の2ファイルを、SPAdesとUnicyclerの両方に共通で使います。

Windows側のPowerShellで、fastpのレポートをダウンロードします。

New-Item -ItemType Directory -Force "$HOME\genome-workshop\qc"
scp -P PORT USER@HOST:/home/bsjinzai/qc/fastp.html "$HOME\genome-workshop\qc"

Windowsのブラウザーでgenome-workshop\qc\fastp.htmlを開き、処理前後のリード数、塩基数、品質、アダプター除去を確認します。fastp公式説明も参照してください。

3. SPAdes を追加して実行する

conda install --override-channels -c conda-forge -c bioconda \
  --strict-channel-priority spades
spades.py --version

spades.py \
  --isolate \
  -1 qc/illumina_1.clean.fastq.gz \
  -2 qc/illumina_2.clean.fastq.gz \
  -o assembly/illumina \
  -t 8 \
  -m 16 \
  > logs/spades.stdout.log 2> logs/spades.stderr.log
conda list --explicit > logs/wgs-short.explicit.txt

--isolate は単離菌の高カバレッジデータ向けの設定です。-m 16 はメモリ上限16 GBです。CPUの並列数は -t 8 です。

4. FASTA とグラフを確認する

ls -lh assembly/illumina/contigs.fasta
ls assembly/illumina/*gfa assembly/illumina/*fastg
conda activate wgs-long
seqkit stats -a -T assembly/illumina/contigs.fasta \
  > results/illumina.assembly.stats.tsv
cat results/illumina.assembly.stats.tsv

比較には contigs.fasta を使います。scaffolds.fasta はギャップを含む接続情報を使った出力です。

Windows側のPowerShellでグラフをダウンロードし、Bandageでgenome-workshop\assembly\illumina\assembly_graph.fastgを開きます。

New-Item -ItemType Directory -Force "$HOME\genome-workshop\assembly\illumina"
scp -P PORT USER@HOST:/home/bsjinzai/assembly/illumina/assembly_graph.fastg "$HOME\genome-workshop\assembly\illumina"

短いリードでは、その長さやペア間の距離を超える反復配列の経路を決めにくくなります。Nanoporeのグラフと比べ、分岐や断片がどの程度見えるか記録します。

参考：SPAdesの実行オプション、SPAdesの出力。