Nanopore リードを確認する | iWellness 実践・研修

1. 配布ファイルを確認する

cd /home/bsjinzai
conda activate wgs-long
ls -lh data/nanopore.fastq.gz
gzip -t data/nanopore.fastq.gz
seqkit head -n 2 data/nanopore.fastq.gz

gzip -t は正常なら何も表示せずに終了します。

FASTQ の1リードは通常4行で表されます。1行目は識別子、2行目は塩基配列、3行目は +、4行目は各塩基の品質情報です。今回は電気信号を配列へ変換済みのリードを使います。

2. 長さと量を調べる

seqkit stats -a -T data/nanopore.fastq.gz \
  > results/nanopore.raw.stats.tsv
cat results/nanopore.raw.stats.tsv

num_seqs はリード数、sum_len は総塩基数、avg_len は平均長、N50 は長さで並べたとき累積塩基数が全体の半分に達するリードの長さです。ここでの N50 はリードの N50であり、後で計算するアセンブリの N50 とは対象が違います。

概算カバレッジは 総塩基数 ÷ 推定ゲノムサイズ です。たとえば300,000,000 bpなら、3,200,000 bpで割って約94倍になります。リードがゲノム上に均一に分布することを仮定した目安であり、各領域の実測depthではありません。

3. Filtlong で約300 Mbに絞る

conda install --override-channels -c conda-forge -c bioconda \
  --strict-channel-priority filtlong

set -o pipefail
filtlong --target_bases 300000000 data/nanopore.fastq.gz \
  2> logs/filtlong.log \
  | gzip > qc/nanopore.filtered.fastq.gz

Filtlong は長さと品質を考慮してリードを選びます。--target_bases は残す総塩基数の目標です。元のデータがこれより少なければ、データ量は増えません。リード単位で選ぶため、結果がちょうど300 Mbになるとは限りません。

set -o pipefail は、パイプの前段が失敗した場合も失敗として扱う Bash の設定です。

4. 前処理後を確認する

gzip -t qc/nanopore.filtered.fastq.gz
seqkit stats -a -T data/nanopore.fastq.gz qc/nanopore.filtered.fastq.gz \
  > results/nanopore.qc.stats.tsv
cat results/nanopore.qc.stats.tsv
conda list --explicit > logs/wgs-long.qc.explicit.txt

リード数、総塩基数、平均長、N50を記録します。この後の Flye と Hybrid には、ここで作った同じフィルタ済みリードを使います。元データは data/ に残します。

qc/nanopore.filtered.fastq.gz が読み取れ、総塩基数が0でないことを確認します。長いリードを優先することで、短い配列由来のリードも減る可能性があります。後のアセンブリで小さな配列が見つからない場合の検討材料になります。

参考：SeqKit、Filtlong。