🧮

LRGE: 長鎖リードのオーバーラップからゲノムサイズを推定する

2025/11/09に公開

はじめに

ゲノムサイズはオミックス解析のデータ取得量や、適切なデータボリュームに合わせるダウンサンプリング、アセンブリやカバレッジ計算などのパラメータ決定に不可欠です。

一方、真核生物や高いヘテロ接合度のサンプルではショートリード前提のk-merベースの推定ツール(GenomeScope2など)は結果が破綻しやすいです。

著者はロングリード同士のオーバーラップ数に着目し、各リードの期待重複数からゲノムサイズを導く LRGEを開発しました。

minimap2で生成したPAFファイル(Pairwise mApping Format)を利用し、(1)リード長と最小重複長から理論的な重複確率を計算、(2)各リードの推定値の中央値を最終推定値とすることで外れ値の影響を抑えています。

LRGEはRustで実装されており、スタンドアロンのCLI(lrge)や、ライブラリ(liblrge)、Bioconda、コンテナ(Docker/Apptainer)が提供されています。デフォルトは2set(ターゲット:10k、クエリ:5k)のアプローチですが、all-vs-all(ava)のアプローチも選択可能です。

論文では、3370の細菌データセット(ONT: 2468 / PacBio: 902) でLRGEをGenomeScope2、Mash、Ravenと比較し、ONTデータではk-mer法より高精度かつRavenと同等の誤差で10倍以上軽量に動作することを示しました。

加えて、3種類の真核生物 (シロイヌナズナ (Arabidopsis thaliana: 125 Mbp)、キイロショウジョウバエ (Drosophila melanogaster: 143 Mbp)、パン酵母 (Saccharomyces cerevisiae: 12 Mbp))でも検証を実施し、ゲノムサイズが大きくなるにつれて、十分なオーバーラップを得るためのリードが必要にはなるが、実用的な精度を達成しています。

使い方

1. Install LRGE

GitHubに記載のインストール手順は下記のとおりです。

https://github.com/mbhall88/lrge

方法	コマンド例	メモ
Precompiled binary	`curl -sSL lrge.mbh.sh \| sh`	`/usr/local/bin` へ配置されるため権限に注意
Bioconda	`conda install -c bioconda lrge`	依存込みで管理可能
Cargo	`cargo install lrge`	Rustを既に使っている環境向け
コンテナ	`apptainer exec docker://ghcr.io/mbhall88/lrge:latest lrge --help`	依存が多いHPCではApptainer推奨

今回はbioconda経由でインストールします。

# Miniforge3の導入
wget "https://github.com/conda-forge/miniforge/releases/download/25.3.0-3/Miniforge3-25.3.0-3-Linux-x86_64.sh"
bash Miniforge3-25.3.0-3-Linux-x86_64.sh -b -p $PWD/miniforge3

# activate
. ${PWD}/miniforge3/etc/profile.d/conda.sh
MAMBA_ROOT_PREFIX=${PWD}/miniforge3
conda activate

mamba create -n lrge_env -c bioconda lrge
mamba activate lrge_env

USAGE LRGE

Genome size estimation from long read overlaps

Usage: lrge [OPTIONS] <INPUT>

Arguments:
  <INPUT>  Input FASTQ file

Options:
  -o, --output <OUTPUT>      Output file for the estimate [default: -]
  -T, --target <INT>         Target number of reads to use (for two-set strategy; default) [default: 10000]
  -Q, --query <INT>          Query number of reads to use (for two-set strategy; default) [default: 5000]
  -n, --num <INT>            Number of reads to use (for all-vs-all strategy)
  -P, --platform <PLATFORM>  Sequencing platform of the reads [default: ont] [possible values: ont, pb]
  -F, --filter-contained     Exclude overlaps for internal matches
  -t, --threads <INT>        Number of threads to use [default: 1]
  -C, --keep-temp            Don't clean up temporary files
  -D, --temp <DIR>           Temporary directory for storing intermediate files
  -s, --seed <INT>           Random seed to use - making the estimate repeatable
  -q, --quiet...             `-q` only show errors and warnings. `-qq` only show errors. `-qqq` shows nothing
  -v, --verbose...           `-v` show debug output. `-vv` show trace output
  -h, --help                 Print help (see more with '--help')
  -V, --version              Print version

論文で引き合いに出ていたツールもインストールしておく

mamba install -c conda-forge -c bioconda Genomescope2 FastK smudgeplot jellyfish kmc Mash

2. デモデータ

ダウンロードフォルダの作成

mkdir 00dat

Mycobacterium tuberculosis(ゲノムサイズ4.40 Mbp / 4405449 bp)

https://elifesciences.org/articles/98300

wget -O 00dat/Mycobacterium_tuberculosis.fastq.gz "ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR283/049/SRR28370649/SRR28370649_1.fastq.gz"

seqkitでデータサイズについて把握しておきます。

seqkit stats -T 00dat/*.gz

file    format  type    num_seqs        sum_len min_len avg_len max_len
00dat/Mycobacterium_tuberculosis.fastq.gz       FASTQ   DNA     187189  715157149       3       3820.5  43938

3. 原核生物のゲノムサイズの推定

mkdir 01Results

Mycobacterium tuberculosis ONT リード

lrgeによるゲノムサイズの推定

lrge -s 123 --float-my-boat --q1 0.15 --q3 0.65 -t 32 00dat/Mycobacterium_tuberculosis.fastq.gz -o 01Results/Mycobacterium_tuberculosis_size.txt

-s: シード値の指定
--float-my-boat: 推定値を最も近い整数に丸めない
-q1,q3: パーセンタイルの値を変更
-n: all-vs-allのリードランダムサンプリングのリード数

cat 01Results/Mycobacterium_tuberculosis_size.txt
# 4468119

# PacBio リードを使用する場合 -Pでプラットフォームを指定する
# lrge -P pb -F -T 20000 -Q 10000 pacbio.fq.gz

-t でスレッド数、-P でプラットフォーム補正 (エラー補正モデル) を指定
真核ゲノムやリピートが多いデータでは -F (internal match 除外) と --use-min-ref (小さい集合を minimap2 の reference に回す) を合わせるとメモリ削減になる
lrge は Fastqから内部的にランダムサブサンプリング → minimap2 実行 → PAF 解析 → 推定値と IQR の報告まで自動で完結します。PAFを外部で作る必要はありません。

ゲノムサイズ推定アプローチとパラメーターの調整

2set vs. ava: 2set はメモリ効率重視。ava は高精度ですが計算量が O(n²) なので、数万リード以上の場合は 2setが推奨。
リード数: リード長とクオリティ高いほど、必要なリード数は減ります。ONT で20kbp以上、QV>12を確保できればデフォルトのままでも5%未満の誤差が期待できるようです。
Depthの偏り: プラスミドや rRNA が極端に高コピーなサンプルでは中央値が小さな repliconに引きずられるため、事前にリードを binning するか、ダウンサンプリング時に長い染色体由来リードを優先的に残すと安定します

Genomescope2によるゲノムサイズの推定

kmcでk-mer頻度のヒストグラムを計算

kmc -v -k21 -ci2 -cs10000 -t32 00dat/Mycobacterium_tuberculosis.fastq.gz 01Results/Mycobacterium_tuberculosis_kmc /tmp

k-merカウントのヒストグラムをエクスポート

kmc_tools transform 01Results/Mycobacterium_tuberculosis_kmc histogram 01Results/Mycobacterium_tuberculosis_kmc_reads.histo -cx10000

genomescopeの実行

genomescope2 -i 01Results/Mycobacterium_tuberculosis_kmc_reads.histo -o 01Results/Mycobacterium_tuberculosis_genomescope -k 21 -p 1

-k: 使用するkmerサイズ
-p: 倍数性を指定

Model converged het:0 kcov:133 err:0.00308 model fit:1.18 len:4359196

推定ゲノムサイズが4359196とlrgeや理論値ともかなり近似した値が得られた。

FastKでkmer頻度を算出してSmudgeplotを書いてみる。

FastK -v -t4 -k21 -M16 -T12 00dat/Mycobacterium_tuberculosis.fastq.gz -N01Results/FastK_Table

smudgeplot.py hetmers -L 12 -t 4 -o 01Results/kmerpairs --verbose 01Results/FastK_Table

smudgeplot.py all -o 01Results/Mycobacterium_tuberculosis 01Results/Mycobacterium_tuberculosis_hetmers_text.smu -ylim 200 -t "Mycobacterium_tuberculosis ONT"

しみ(=smudge)として表されるhigh kmerな部分が1箇所ABで存在しているので想定通りの結果になっていると思われます。

Mashによるゲノムサイズの推定

mash sketch -s 100000 -m 10 -r 00dat/Mycobacterium_tuberculosis.fastq.gz -o 01Results/Mycobacterium_tuberculosis

-r: readsモード。総塩基数ではなく k-mer 内容からゲノムサイズを推定して、内部でp値計算などに使用
- 実行ログにEstimated genome sizeが出力される。
-m: 除外するkmer数値の最小値を指定
-k, -s: k-mer長とスケッチサイズ。小さすぎると荒く、大きくすると精度は向上するが計算コストも向上する

Estimated genome size: 4.50523e+06
Estimated coverage:    133.206

こちらも4.50523e+06と理論値に近い値が得られた。また、1コマンドで10秒程度で結果が得られることを踏まえると、簡易的に結果を取得したい場合に良さそう。

おしまい。