🧬

SqueezeMeta v1.7.0: 多機能なメタゲノム解析パイプライン

2025/04/26に公開

はじめに

こちらで紹介したSqueezeMetaが1.7.0にバージョンアップされました。

またチュートリアルも更新されてsphinxで作成されています。

https://squeezemeta.readthedocs.io/en/v1.7.0/index.html

全体的にアップデート箇所が多いようなので再度インストールしてテスト実行してみます。

SqueezeMeta

SqueezeMetaは、分析の複数のステップでカスタムスクリプトと外部ソフトウェアパッケージを組み合わせて使用します。

01.run_assembly.pl: アセンブル
02.run_barrnap.pl: RNAの予測と分類
03.run_prodigal.pl: ORF (CDS) 予測
04.rundiamond.pl: 系統および遺伝子機能データベースに対する相同性検索
05.run_hmmer.pl: Pfamデータベースを使ったHmmerによる検索
06.lca.pl: 遺伝子の系統分類
07.fun3assign.pl: 遺伝子の機能割り当て（オプション）
08.blastx.pl: 遺伝子予測やヒットがないコンティグの部分に対するBlastx
09.summarycontigs3.pl: コンティグの系統分類結果と分類学的差異のチェック
10.mapsamples.pl: 遺伝子とコンティグのカバレッジと存在量推定
11.mcount.pl: 分類群の存在量の推定
12.funcover.pl: 機能の存在量の推定
13.mergeannot2.pl: 以前の結果を統合してORFテーブルを取得
14.runbinning.pl: ビニング
15.dastool.pl: DAS Toolとによるビンの統合
16.addtax2.pl: ビンの系統分類と分類上の相違のチェック
17.checkbins.pl: CheckM2によるビンのクオリティチェック (オプションでGTDB-Tkによる分類)
18.getbins.pl: 以前の結果をマージしてビンテーブルを取得
19.getcontigs.pl: 以前の結果をマージしてコンティグテーブルを取得
20.minpath.pl: 各ビンのKEGGおよびMetaCyc pathwayの予測
21.stats.pl: 最終的な統計量などの取得
sqm2tables.py: 系統分類および機能的プロファイルを含むテーブルの生成

SqueezeMeta v1.6からの主な変更点

GitHubのリリースノートより一部抜粋

SqueezeMeta 互換性
- CheckM2のデータベース導入必須:
  - 既存ユーザーはCheckM2用のファイルを追加でダウンロード
SqueezeMeta 新機能
- -extbins: 既存のゲノムやBinのアノテーションおよび定量解析が可能に (-extassemblyのMulti Fasta対応版のような位置づけ)
- -onlybins: アノテーションをスキップすることでビン取得までの処理を高速化することが可能
- --gtdbtk: GTDB-tkによるビンの分類が可能 (データベースは別途ダウンロード)
- CheckMの導入: ビンの完全性・コンタミ率の評価 (CheckM1の旧分類問題を解消)
- --fastnr: DIAMONDの--fastモードで高速アノテーション (精度と速度はトレードオフ)
- コンティグ・ビンのばらつき計算が簡略化
- sqm2tables.py: 自動的に最終ステップで実行されるように
- 依存パッケージはCondaベースに
SqueezeMeta 変更・バグFix
- プロジェクト名がコンティグ名・ビン名の接頭語として追加されるように
- Step10で生成されるBAMにサンプル情報を追加
- ~~DIAMONDのバージョンを2.1.9に更新~~
- Step13実行時にORFが重複して出力されるバグの修正
- MetaBAT2の新バージョンと互換性のない問題を修正
- SPAdes実行時、従来のcontigs.fastaではなくscaffolds.fastaを使用（rnaspadesの場合はtranscripts.fasta）
SQMtools
- loadSQM: 複数のSqueezeMetaのプロジェクトを同時読み込み可能に
- ビンの定義・修正・手動キュレーション機能を追加。手動またはsubset関数でビンの再編集と品質再評価が可能
- exportContigs, exportORFs, exportBinsで配列データを簡単にエクスポート可能
- コピー数の計算基準をRecA遺伝子から、10種類の普遍的なSingleコピー遺伝子の中央値へ変更(single_copy_genesパラメータで変更可)
SQMtools 変更・バグFix
- load_sequences=FALSE に設定することで、低メモリでの読み込みが可能
- exportPathway関数で出力先ディレクトリを指定できる output_dir 引数を追加
- ORF の開始・終了位置を明示的にトラッキング
- plotFunctionsやexportPathwayのデフォルト定量法をcopy_numberに変更
- combineSQMlite実行後に一部IDが欠損する問題を修正
- data.tableパッケージが自動的に読み込まれない不具合を修正
- ORFやコンティグが1つしかない場合のsubset実行時のエラーを修正

SqueezeMetaのインストール

作業ディレクトリを作成します。

mkdir squeezemeta

今回はConda/Mambaでインストールします。

Miniforge3の導入

カレントディレクトリにminiforge3のインストーラーをダウンロードします。その後、インストーラーを実行してminiforge3をインストールします。既にconda環境がある場合でも問題ないです。

wget "https://github.com/conda-forge/miniforge/releases/lasqm_results/download/Miniforge3-$(uname)-$(uname -m).sh"
bash Miniforge3-$(uname)-$(uname -m).sh -b -p $PWD/miniforge3

導入したminiforge3のデフォルト環境にactivateしてから、SqueezeMeta用の仮想環境を作成します。

. ${PWD}/miniforge3/etc/profile.d/conda.sh
. ${PWD}/miniforge3/etc/profile.d/mamba.sh

conda activate

導入したminiforge3のcondaにパスが通っているか確認します。

conda activate
conda info -e

*マークが指定したパスのbase環境を指していればOKです。

SqueezeMeta用の仮想環境の作成

mamba createで仮想環境を作成します。

インストール時の注意点
v1.7.0ではexport CONDA_OVERRIDE_CUDA=12.4と、conda config --set remote_connect_timeout_secs 60の実行が必要な場合があるようです。

以下のIssueのようなエラーが出たら試してみてください。

https://github.com/jtamames/SqueezeMeta/issues/962

export CONDA_OVERRIDE_CUDA=12.4
conda config --set remote_connect_timeout_secs 60

仮想環境の作成と同時にSqueezeMetaをインストールします。また、Host除去用にkneaddata, 並列圧縮用にpigzをインストールします。

mamba create -n SqueezeMeta -c conda-forge -c bioconda -c fpusan squeezemeta=1.7 kneaddata pigz --no-channel-priority --override-channels

仮想環境のActivate

mamba activate SqueezeMeta

USAGE SqueezeMeta.pl

Usage: SqueezeMeta.pl -m <mode> -p <project name> -s <samples file> -f <sequence dir> [options]

Arguments:

 Mandatory parameters:
   -m <mode>: Mode (sequential, coassembly, merged, seqmerge) (REQUIRED)
   -s|-samples <samples file>: Samples file (REQUIRED)
   -f|-seq <sequence dir>: fastq/fasta read files' directory (REQUIRED)
   -p <project name>: Project name (REQUIRED in coassembly and merged modes)
  
 Restarting
   --restart: Restarts the given project where it stopped (project must be speciefied with -p option) (will NOT overwite previous results, unless --force-overwrite is also provided)
   -step <step number>: In combination with --restart, restarts the project starting in the given step number (combine with --force_overwrite to regenerate results)
   --force_overwrite: Do not check for previous results, and overwrite existing ones
   
 Filtering: 
   --cleaning: Filters with Trimmomatic (Default: No)
   -cleaning_options [options]: Options for Trimmomatic (Default:LEADING:8 TRAILING:8 SLIDINGWINDOW:10:15 MINLEN:30)
   
 Assembly: 
   -a: assembler <megahit, spades, rnaspades, spades-base, canu, flye> (Default: megahit)
   -assembly_options [options]: Extra options to be passed when calling the mapper
   -c|-contiglen <size>: Minimum length of contigs (Default: 200)
   -extassembly <file>: External assembly, path to a fasta file with contigs (overrides all assembly steps).
   --sg|--singletons: Add unassembled reads to the contig file, as if they were contigs  
   -contigid <string>: Nomenclature for contigs (Default: assembler´s name)
   --norename: Don't rename contigs (Use at your own risk, characters like '_' in contig names will make it crash)
   
 Mapping: 
   -map: mapping software <bowtie, bwa, minimap2-ont, minimap2-pb, minimap2-sr> (Default: bowtie) 
   -mapping_options [options]: Extra options to be passed when calling the mapper

 ONT support: 
   --minion: Run on MinION reads (assembler: canu; mapper: minimap2-ont; consensus: 20) (Default: no)

 Annotation:
   -g <int>: Number of targets for DIAMOND global ranking during taxonomic assignment (Default: 100) 
   -db <file>: Specify a new taxonomic database
   --nodiamond: Check if Diamond results are already in place, and just in that case skips the Diamond run (Default: no)
   --nocog: Skip COG assignment (Default: no)
   --nokegg: Skip KEGG assignment (Default: no)
   --nopfam: Skip Pfam assignment  (Default: no)
   --fastnr: Run DIAMOND in --fast mode for taxonomic assignment (Default: no)
   --euk: Drop identity filters for eukaryotic annotation (Default: no)
   -consensus <value>: Minimum percentage of genes for a taxon needed for contig consensus (Default: 50)
   --D|--doublepass: First pass looking for genes using gene prediction, second pass using Diamond BlastX  (Default: no)
   -extdb <database file>: List of user-provided databases
   -b|-block-size <block size>: block size for diamond against the nr database (Default: calculate automatically)
   
 Binning:
   --nobins: Skip all binning  (Default: no). Overrides -binners
   --onlybins: Run only assembly, binning and bin statistics (including GTDB-Tk if requested) (Default: no)
   -binners: Comma-separated list with the binning programs to be used (available: maxbin, metabat2, concoct)  (Default: concoct,metabat2)
   -taxbinmode <s,c,s+c,c+s>: Source of taxonomy annotation of bins (s: SqueezeMeta; c: CheckM; s+c: SqueezeMeta+CheckM;  c+s: CheckM+SqueezeMeta; (Default: s). THIS HAS BEEN DEPRECATED, USE --gtdbtk INSTEAD IF YOU NEED A MORE PRECISE BIN TAXONOMY.
   --nomarkers: Skip retrieval of universal marker genes from bins (you will still get completeness/contamination estimates, but won't be able to do bin refining in SQMtools)
   --gtdbtk: Run GTDB-Tk to classify the bins. Requires a working GTDB-Tk installation in available in your environment
   -gtdbtk_data_path: Path to the GTDB database, by default it is assumed to be present in /path/to/SqueezeMeta/db/gtdb
   -extbins: Path to a directory containing external genomes/bins provided by the user. There must be one file per genome/bin, containing each contigs in the fasta format. This overrides the assembly and binning steps

 
 Performance:
   -t <threads>: Number of threads (Default: 12)
   -canumem <mem>: memory for canu in Gb (Default: 32)
   --lowmem: run on less than 16Gb of memory (Default:no)

 Other:
   -test <step>: Running in test mode, stops AFTER the given step number
   --empty: Creates a empty directory structure and conf files, does not run the pipeline
   
 Information:
   -v: Version number  
   -h: This help

kneaddata用のTrimmomaticはReleaseページからダウンロードして解凍しておきます。

wget https://github.com/usadellab/Trimmomatic/files/5854859/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip 
rm Trimmomatic-0.39.zip

# 確認
java -jar Trimmomatic-0.39/trimmomatic-0.39.jar -version

テストスクリプトの実行を実行します。

test_install.pl

Checking the OS
        linux OK

Checking that tree is installed
        tree --help OK

Checking that ruby is installed
        ruby -h OK

Checking that java is installed
        java -h OK

Checking that all the required perl libraries are available in this environment
        perl -e 'use Term::ANSIColor' OK
        perl -e 'use DBI' OK
        perl -e 'use DBD::SQLite::Constants' OK
        perl -e 'use Time::Seconds' OK
        perl -e 'use Tie::IxHash' OK
        perl -e 'use Linux::MemInfo' OK
        perl -e 'use Getopt::Long' OK
        perl -e 'use File::Basename' OK
        perl -e 'use DBD::SQLite' OK
        perl -e 'use Data::Dumper' OK
        perl -e 'use Cwd' OK
        perl -e 'use XML::LibXML' OK
        perl -e 'use XML::Parser' OK
        perl -e 'use Term::ANSIColor' OK

Checking that all the required python libraries are available in this environment
        python3 -h OK
        python3 -c 'import numpy' OK
        python3 -c 'import scipy' OK
        python3 -c 'import matplotlib' OK
        python3 -c 'import dendropy' OK
        python3 -c 'import pysam' OK
        python3 -c 'import Bio.Seq' OK
        python3 -c 'import pandas' OK
        python3 -c 'import sklearn' OK
        python3 -c 'import nose' OK
        python3 -c 'import cython' OK
        python3 -c 'import future' OK

Checking that all the required R libraries are available in this environment
        R -h OK
        R -e 'library(doMC)' OK
        R -e 'library(ggplot2)' OK
        R -e 'library(data.table)' OK
        R -e 'library(reshape2)' OK
        R -e 'library(pathview)' OK
        R -e 'library(DASTool)' OK
        R -e 'library(SQMtools)' OK

Checking binaries
        spades.py OK
        metabat2 OK
        jgi_summarize_bam_contig_depths OK
        samtools OK
        bwa OK
        minimap2 OK
        diamond OK
        hmmsearch OK
        cd-hit-est OK
        kmer-db OK
        aragorn OK
        mothur OK

Checking that SqueezeMeta is properly configured...

        SqueezeMeta doesn't know where the databases are located!

        If you didn't download them yet please run:

        * perl <Path to dir>/squeezemeta/miniforge3/envs/SqueezeMeta/SqueezeMeta/utils/install_utils/make_databases.pl /download/path/
                (To download the latest source data and compile the databases in your server)

        * perl <Path to dir>/squeezemeta/miniforge3/envs/SqueezeMeta/SqueezeMeta/utils/install_utils/download_databases.pl /download/path/
                (To download a pre-compiled version of the database, which is much quicker)

        If the databases are already present in your server, you can configure this installation of SqueezeMeta with:

        * perl <Path to dir>/squeezemeta/miniforge3/envs/SqueezeMeta/SqueezeMeta/utils/install_utils/configure_nodb.pl /path/to/db
                (Where /path/to/db is the path to the "db" folder generated when downloading the databases)


------------------------------------------------------------------------------

WARNING: Some SqueezeMeta dependencies could not be found in your environment or failed to execute!
        - Databases were not installed, or are not configured in this installation of SqueezeMeta

データベースが無いためWARNINGが出ています。次の章でデータベースの構築を実施します。

データベースのダウンロードと構築

SqueezeMetaの処理過程で系統分類にはGenBank nr、機能分類にはeggnogやKEGG、Pfamといったように複数のデータベースを使用します。

これらを取得するためにdownload_databases.plを実行します。SqueezeMeta に必要なすべてのデータベースをフォーマット済みの状態で取得可能です。

# フォルダの作成
mkdir sqm_db
# ダウンロードスクリプトの実行
download_databases.pl ./sqm_db

私は上記コマンドでダウンロードが終了しなかったので、下記FTPサイトからSqueezeMetaDB.tar.gzを直接ダウンロードしました。

http://andes.cnb.csic.es/SqueezeMeta/

mkdir sqm_db

array=(SqueezeMetaDB.md5 SqueezeMetaDB.tar.gz)

# Download
for i in ${array[@]}; do
    wget http://andes.cnb.csic.es/SqueezeMeta/${i} -o sqm_db/${i}
done

# MD5値の確認
md5sum -c SqueezeMetaDB.md5

tarファイルの解凍

tar -xvf sqm_db/SqueezeMetaDB.tar.gz -C sqm_db

configure_nodb.pl <Path to db>/db

足りないファイル？が1つ仮想環境内にダウンロードされましたが最終的にエラーが出ていなければOKです。

Checking that SqueezeMeta is properly configured... checking database in <Path to db>/db
        nr.db OK
        CheckM manifest OK
        CheckM database OK
        LCA_tax DB OK
        Database was built on Fri Jun  7 18:05:02 2024

GTDB

次にオプションで使用可能なGTDBのデータベースファイルをダウンロードします。こちらも140GB程度になるのでストレージの容量には注意してください。

まず、データベースとgtdb-tkで互換性があるためインストールされているGTDB-tkのバージョンを確認します。

gtdbtk -v
# gtdbtk: version 2.4.0 Copyright 2017 Pierre-Alain Chaumeil, Aaron Mussig and Donovan Parks

v2.4.0は2025/04時点で最新のR226まで対応していたのでR226をダウンロードします。

mkdir sqm_db/db/gtdb
wget https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/auxillary_files/gtdbtk_package/full_package/gtdbtk_data.tar.gz -o sqm_db/gtdb/gtdbtk_data.tar.gz

# ミラーサイトを使う場合は以下を実行
# wget https://data.gtdb.ecogenomic.org/releases/latest/auxillary_files/gtdbtk_package/full_package/gtdbtk_data.tar.gz (mirror for Australia)

気長に待ちます。ダウンロードが完了したら解凍しておきます。

tar xvzf sqm_db/gtdb/gtdbtk_data.tar.gz -C sqm_db/db/gtdb -I pigz

デモデータのダウンロード

以下の論文のデータを使用します。

Relating the metatranscriptome and metagenome of the human gut

mkdir 00RawFastq
array=(SRR769395 SRR769396 SRR769397 SRR769398 SRR769399 SRR769400 SRR769401 SRR769402 SRR769403 SRR769404 SRR769405 SRR769406 SRR769407 SRR769408 SRR769409 SRR769410 SRR769411 SRR769412 SRR769413 SRR769414 SRR769415 SRR769416 SRR769417 SRR769418 SRR769419 SRR769420 SRR769421 SRR769422 SRR769423 SRR769424 SRR769425 SRR769426 SRR769427 SRR769428 SRR769429 SRR769430 SRR769431 SRR769432 SRR769433 SRR769434 SRR769435 SRR769436 SRR769437 SRR769438 SRR769439 SRR769440 SRR769441 SRR769442 SRR769509 SRR769510 SRR769511 SRR769512 SRR769513 SRR769514 SRR769515 SRR769516 SRR769518 SRR769519 SRR769520 SRR769521 SRR769522 SRR769523 SRR769524 SRR769525 SRR769526 SRR769527 SRR769528 SRR769529 SRR769530 SRR769531 SRR769532 SRR769533 SRR769534 SRR769535 SRR769536 SRR769537 SRR769538 SRR769539 SRR769540)

for i in ${array[@]}; do
  curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR769/${i}/${i}_1.fastq.gz -o 00RawFastq/${i}_1.fastq.gz 
  curl -L ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR769/${i}/${i}_2.fastq.gz -o 00RawFastq/${i}_2.fastq.gz 
done

データQC

メタゲノムのデータに対して上記の記事と同様の処理を実施します。

mkdir 01Ref
wget https://huttenhower.sph.harvard.edu/kneadData_databases/Homo_sapiens_hg39_T2T_Bowtie2_v0.1.tar.gz -o 01Ref/Homo_sapiens_hg39_T2T_Bowtie2_v0.1.tar.gz

# 解凍
tar xvf 01Ref/Homo_sapiens_hg39_T2T_Bowtie2_v0.1.tar.gz -C 01Ref/ 
rm 01Ref/Homo_sapiens_hg39_T2T_Bowtie2_v0.1.tar.gz

以下のコマンドでkneaddataを実行します。

mkdir 02kneaddata
array=(SRR769509 SRR769510 SRR769511 SRR769512 SRR769513 SRR769514 SRR769515 SRR769516 SRR769518 SRR769519 SRR769520 SRR769521 SRR769522 SRR769523 SRR769524 SRR769525 SRR769526 SRR769527 SRR769528 SRR769529 SRR769530 SRR769531 SRR769532 SRR769533 SRR769534 SRR769535 SRR769536 SRR769537 SRR769538 SRR769539 SRR769540)

# Set the number of parallel jobs (adjust as needed)
for acc in ${array[@]};do
    kneaddata \
    --trimmomatic Trimmomatic-0.39/ \
    -t 32 \
    -i1 "00RawFastq/${acc}_1.fastq.gz" \
    -i2 "00RawFastq/${acc}_2.fastq.gz" \
    --reference-db "01Ref/bowtie2-index" \
    --output-prefix "${acc}" \
    --trimmomatic-options "ILLUMINACLIP:$PWD/Trimmomatic-0.39/adapters/NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:4:20 MINLEN:50" \
    -o "02kneaddata/${acc}"

    pigz -p 32 ${PWD}/02kneaddata/${acc}/*_paired_1.fastq && \
    pigz -p 32 ${PWD}/02kneaddata/${acc}/*_paired_2.fastq && \
    rm -rv ${PWD}/02kneaddata/${acc}/*.fastq # 中間ファイルを消す
done

RNAのデータに対してはfastpでQCします。fastpはsingularityコンテナにて取得しました。

singularity build --fakeroot def/fastp:0.24.0.sif docker://quay.io/biocontainers/fastp:0.24.0--heae3180_1

array=(SRR769395 SRR769396 SRR769397 SRR769398 SRR769399 SRR769400 SRR769401 SRR769402 SRR769403 SRR769404 SRR769405 SRR769406 SRR769407 SRR769408 SRR769409 SRR769410 SRR769411 SRR769412 SRR769413 SRR769414 SRR769415 SRR769416 SRR769417 SRR769418 SRR769419 SRR769420 SRR769421 SRR769422 SRR769423 SRR769424 SRR769425 SRR769426 SRR769427 SRR769428 SRR769429 SRR769430 SRR769431 SRR769432 SRR769433 SRR769434 SRR769435 SRR769436 SRR769437 SRR769438 SRR769439 SRR769440 SRR769441 SRR769442)

# Set the number of parallel jobs (adjust as needed)
num_jobs=2
path1=00RawFastq
path2=02kneaddata

mkdir ${path2}

process_sample() {
    fq_name=$1
    path1=$2
    path2=$3
    echo "Processing sample: ${fq_name}"

    mkdir ${path2}/${fq_name}

    ./def/fastp:0.24.0.sif fastp \
    --in1 ${path1}/${fq_name}_1.fastq.gz \
    --in2 ${path1}/${fq_name}_2.fastq.gz \
    --out1 ${path2}/${fq_name}/${fq_name}_paired_1.fastq.gz \
    --out2 ${path2}/${fq_name}/${fq_name}_paired_2.fastq.gz \
    --html ${path2}/${fq_name}/${fq_name}_report.html \
    --json ${path2}/${fq_name}/${fq_name}_report.json \
    --qualified_quality_phred 30 \
    --length_required 50 \
    --detect_adapter_for_pe \
    --trim_poly_g \
    --cut_front \
    --thread 16

    echo "Finished processing sample: ${fq_name}"
}

# Export the function for parallel execution
export -f process_sample

# Run fastp for each sample in parallel
parallel -j $num_jobs process_sample \
::: "${array[@]}" \
::: "${path1}" \
::: "${path2}"

fastqだけ移動させておきます。

mv 02kneaddata/*/*.fastq.gz 02kneaddata/

サンプルリストファイル

headerなしのタブ区切り形式でサンプル情報が記載されたファイルを作成します。

解析の実行

SqueezeMetaを実行します。

SqueezeMeta.pl -m coassembly -p sqm_results -s samplelist.tsv -f 02kneaddata -b 9 -t 32 --gtdbtk -gtdbtk_data_path <Path to gtdb database dir>

必須パラメーターは以下の通り。

Option	Argument
-m	sequencial, coassembly, merged, seqmerge>
-p	プロジェクト名 (コアセンブリおよびマージモードでは必須)
-s/-samples	サンプルファイルパス
-f/-seq	Fastqファイルのディレクトリ
--gtdbtk	gtdbtkの実行フラグ
-gtdbtk_data_path	GTDBのデータベースフォルダパス

上記の他に-bと-tオプションを指定しています。-tは使用するスレッド数です。-bはHomology searchで使用するDIAMONDのnrデータベースに対するブロックサイズです。ブロックサイズは最大で16で自動計算されますが、128GBのメモリでは足りなかったため意図的に小さくなるよう指定しました。

結果

実行から約1日と2時間で終了しました。結果を含めたLogや中間ファイルはsqm_resultsフォルダに出力されます。その中にあるresultsがSqueezeMetaの各解析結果がまとめられたフォルダです。

tree -L 1 sqm_results/results

sqm_results/results
├── 01.sqm_results.fasta
├── 02.sqm_results.16S.txt
├── 02.sqm_results.rnas
├── 02.sqm_results.trnas
├── 02.sqm_results.trnas.fasta
├── 03.sqm_results.faa
├── 03.sqm_results.fna
├── 03.sqm_results.gff
├── 06.sqm_results.fun3.tax.noidfilter.wranks
├── 06.sqm_results.fun3.tax.wranks
├── 07.sqm_results.fun3.cog
├── 07.sqm_results.fun3.kegg
├── 07.sqm_results.fun3.pfam
├── 10.sqm_results.mappingstat
├── 11.sqm_results.mcount
├── 12.sqm_results.cog.funcover
├── 12.sqm_results.kegg.funcover
├── 13.sqm_results.orftable
├── 18.sqm_results.bintable
├── 19.sqm_results.contigtable
├── 20.sqm_results.kegg.pathways
├── 20.sqm_results.metacyc.pathways
├── 21.sqm_results.stats
├── bins
└── tables

3 directories, 23 files

基本的な出力は公式Documentや前回の記事を参考にしてください。今回は-nokeggのオプションをつけていないのでkeggやその関連結果も出力されています。

Stats

論文 Fig1

SRAのmetadataから抽出した各サンプルのグループ情報を使ってコンティグへのマッピングに関連するStatsを確認します。グループにはシーケンス法、サンプルタイプ、保存法を設定しています (SalivaのFrozenのデータSRAになかったような。。)

上記のグループ情報をcompsheet.csvとしてカレントディレクトリに保存します。

データを読み込んでグループ情報を結合

import pandas as pd
impoort seaborn

# リードからORFへのマッピングStats
mappingstat_df = pd.read_csv("sqm_results/results/10.sqm_results.mappingstat", skiprows=1, sep='\t')
mappingstat_df = mappingstat_df[~mappingstat_df["# Sample"].str.contains("#")]
mappingstat_df = mappingstat_df.rename(columns={"# Sample":"SampleID"})
# compsheetの読み込み
compsheet_df = pd.read_csv("compsheet.csv", header=0)
# 結合
merged_df = pd.merge(mappingstat_df, compsheet_df, on="SampleID")

StoolサンプルのWGSにおけるMapping率とMapped readsの関係性

plotスクリプト

sns.pairplot(
  data=merged_df.query("Preserve != 'Saliva' & Type=='stool' & Assay=='WGS'"),
  vars=['Mapping perc', 'Mapped reads'],
  hue='Preserve',
  kind='reg',
  diag_kind='kde',
  size=3,
  plot_kws={
      'scatter_kws': {'alpha': 0.4},
      'line_kws': {'linestyle': '--', 'linewidth': .7}
      })

WGSであれば保存法はどれでであってもコンティグへのMapping率は高い値が得られそうですね。

StoolサンプルのRNA-SeqにおけるMapping率とMapped readsの関係性

plotスクリプト

sns.pairplot(
  data=merged_df.query("Preserve != 'Saliva' & Type=='stool' & Assay=='RNA-Seq'"),
  vars=['Mapping perc', 'Mapped reads'],
  hue='Preserve',
  kind='reg',
  diag_kind='kde',
  size=3,
  plot_kws={
      'scatter_kws': {'alpha': 0.4},
      'line_kws': {'linestyle': '--', 'linewidth': .7}
      })

Mapping率が高ければそれだけ各菌種の発現情報が得られているということになります。

結果を見てみるとトランスクリプトームはRNAlaterが保存性も高くユニークなリードを保持しやすい一方、冷凍サンプル(Whole)は収量が得られにくいサンプルがあったのかMapping率も2峰性の分布を示しマッピングされ他リードもばらつきがありデータが安定しなさそうな雰囲気を感じます。凍結融解が結果に影響を及ぼしているしているんでしょうか。

簡単なBinの確認

評価を行うにあたって何らかの基準があると望ましいです。微生物のゲノムの品質評価ツールとして広く使われているCheckMはMIMAG/MISAGの規格を採用しています(CheckM2はモデルに基づく予測を行うので系統に依存しません)。

Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea
(細菌および古細菌の単一増幅ゲノム（MISAG）およびメタゲノムアセンブルゲノム（MIMAG）に関する最小限の情報)

MIMAGの品質カテゴリは複数の要素で構成されていますが、CompletenessとContaminationに着目した場合以下のように定義されています。

カテゴリ	Completeness	Contamination
High-quality draft	≥90%	≤5%(Purity ≥95%)
Medium-quality draft	≥50%	≤10%(Purity ≥90%)
Low-quality draft	<50%	任意

Completenessは、既知の「必須マーカー遺伝子（single-copy marker genes）」の検出率に基づいて推定されます。そのBinに属するはずのゲノム情報がどれだけ揃っているかを示すイメージです。

Contaminationは本来 1 回だけ存在するはずのマーカー遺伝子が複数回見つかった場合、それは別のゲノムの断片が混じっている（＝汚染されている）可能性があるとみなします。その Bin に、他の生物種や系統の断片がどれだけ混入しているかを示すイメージです。

SqueezeMetaも内部でCheckMが動いているので採用することに特に問題は無いでしょう。

今回の結果をseabornで視覚化してみました。

作図用pythonスクリプト

import pandas as pd
import seaborn as sns
import numpy as np
import matplotlib.pyplot as plt

# 使用する列名
use_col=["Bin ID","Completeness", "Contamination"]
# Dataframeの作成
df = pd.read_table("sqm_results/results/18.sqm_results.bintable", skiprows=1, sep="\t", index_col=False, usecols=use_col)
# 純度に変換
df["Purity"] = 100 - df["Contamination"]

# 条件
conditions = [
    (df['Completeness'] >= 90 ) & (df['Purity'] >= 95),
    (df['Completeness'] >= 50 ) & (df['Purity'] >= 90),
]
choices = ['High','Mid']
df['Category'] = np.select(conditions, choices, default = 'Low')

# plot
g=sns.jointplot(data=df, x= "Completeness", y="Purity", kind = "scatter", hue="Category", hue_order =['High','Mid','Low'])
g.ax_joint.set_xlim(-5,105)
g.ax_joint.set_ylim(-5,105)

df["Category"].value_counts()

High: 15
Mid: 101
Low: 182

データ薄めのショートリードの結果でもあるので、High-qualityに属するものはかなり少ないですね。

ここから数値的基準や分類情報などを包括的に評価したものがMetagenome Assembled Genome (=MAG) として扱われます(たぶん)。

sqm2ipath.pl

iPath(Interactive Pathways Explorer)は、KEGGやMetaCycなどの代謝経路データベースをもとにした経路マップの視覚化ツールです。検出された酵素や経路を視覚的に表示するために用いられます。

https://pathways.embl.de/

sqm2ipath.plはSqueezeMetaの出力をiPathでプロットできる形式に変換します。

mkdir ipath_output
sqm2ipath.pl -p sqm_results -classification kegg -o ipath_output/sqm_results.kegg.ipath.out
sqm2ipath.pl -p sqm_results -classification cog -o ipath_output/sqm_results.cog.ipath.out

-taxonで対象を絞ったり出来ますが今回は実施していません。

sqm_results.cog.ipath.outとsqm_results.kegg.ipath.outファイルが出力されればOKです。

iPathを用いた視覚化

iPathは4つの経路マップをwebツールで提供しています。

Metaboilsm: 細胞内の主要な代謝経路（解糖系、TCA回路、アミノ酸・脂質代謝など）を表示
Secondary metabolites: 二次代謝産物の合成経路に特化したマップ
Antibiotics: さまざまな抗生物質の生合成に関与する経路の概要
Microbial metabolism: 多様な環境における微生物代謝経路

Antibiotics: COG

Antibiotics: KEGG

Microbial metabolism: COG

Microbial metabolism: KEGG

ORFの高感度検出オプション

SqueezeMetaの実行にはProdigalによるORF予測の精度を補完・向上させる-D(doublepass)オプションが利用できます。通常のORF予測で見逃された遺伝子をBlastXで補う機能です。

処理配下の流れで行われます

通常のアノテーション後、アノテーションが付かなかった領域やヒットしなかった ORF を抽出
それらの領域に対して DIAMOND BlastX を実行
新たにヒットした配列を追加ORFとして扱い、パイプラインを継続

以下のような利点があります

欠損ORF誤予測の補正が可能(例: フレームシフトによる予測ミス)
特に真核生物やウイルスが含まれるメタゲノムで有効
原核メタゲノムではORF数が約2–3%増加する見込み

デメリットはDIAMONDによるBlastXの追加実行に伴い解析時間が増えることです。

-DをつけてSqueezeMetaを再実行しました。約1日と5時間で終了しました。

SqueezeMeta.pl -m coassembly -p sqm_results_D -s samplelist.tsv -f 02kneaddata -b 9 -t 32 --gtdbtk -gtdbtk_data_path <Path to gtdb database dir>

データの読み込み

script

import polars as pl

df_D = pl.read_csv(
    "sqm_results_D/results/13.sqm_results_D.orftable",
    skip_rows=1,
    separator='\t',
    null_values=["NA", "N/A", "null", ""])

df = pl.read_csv(
    "sqm_results/results/13.sqm_results.orftable",
    skip_rows=1,
    separator='\t',
    null_values=["NA", "N/A", "null", ""])

ORF数の確認

# ORF数
print(f"Number of ORFs: {df.height}")
print(f"Number of ORFs with -D: {df_D.height}")

Number of ORFs: 1812245
Number of ORFs with -D: 1848985

2%程度、ORF数が向上しています。

# Methodごと
print(df['Method'].sort().value_counts())
print(df_D['Method'].sort().value_counts())

Method	count	count (with -D)
Aragorn	7662	5604
Prodigal	1802171	1791275
barrnap	2412	2407
blastx	-	49699

AragornはtRNAやtmRNAの推定ツールです。barrnapはrRNAの探索ツールです。ProdigalとblastxはCDSの探索を目的に使用されています。

blastxと合わせると-Dオプションの付与で2%程度のORF数の向上が確認できました。

print(f"Number of ORFs (with assigned gene): {df['ORF ID'].filter(df['Gene name'].is_not_null()).len()}")
print(f"Number of ORFs with -D (with assigned gene): {df_D['ORF ID'].filter(df_D['Gene name'].is_not_null()).len()}")

Number of ORFs (with assigned gene): 933670 (51%)
Number of ORFs with -D (with assigned gene): 937656 (50%)

遺伝子が割り当てられているORFは半数程度で-Dなしと比べるとほぼ同等数(4000程度向上)。

KEGG IDが割り振られているもの

print(f"Number of ORFs (with assigned KEGG ID): {df['ORF ID'].filter(df['KEGG ID'].is_not_null()).len()}")
print(f"Number of ORFs with -D (with assigned KEGG ID): {df_D['ORF ID'].filter(df_D['KEGG ID'].is_not_null()).len()}")

Number of ORFs (with assigned KEGG ID): 944801 (52%)
Number of ORFs with -D (with assigned KEGG ID): 950982 (51%)

COG IDが割り振られているもの

print(f"Number of ORFs (with assigned COG ID): {df['ORF ID'].filter(df['COG ID'].is_not_null()).len()}")
print(f"Number of ORFs with -D (with assigned COG ID): {df_D['ORF ID'].filter(df_D['COG ID'].is_not_null()).len()}")

Number of ORFs (with assigned COG ID): 1421565 (78%)
Number of ORFs with -D (with assigned COG ID): 1437558 (77%)

Pfamの情報が割り振られているもの。

print(f"Number of ORFs (with assigned Pfam): {df['ORF ID'].filter(df['PFAM'].is_not_null()).len()}")
print(f"Number of ORFs with -D (with assigned Pfam): {df_D['ORF ID'].filter(df_D['PFAM'].is_not_null()).len()}")

Number of ORFs (with assigned Pfam): 796121 (43%)
Number of ORFs with -D (with assigned Pfam): 796122 (43%)

機能情報が1つでも得られそうなデータの数 (COGFUN, KEGG, PFAM)

df_FUN=df.filter(df['COGFUN'].is_not_null() | df['KEGGFUN'].is_not_null() | df['PFAM'].is_not_null())
print(f"Number of ORFs: {df_FUN['ORF ID'].len()}")
df_D_FUN=df_D.filter(df_D['COGFUN'].is_not_null() | df_D['KEGGFUN'].is_not_null() | df_D['PFAM'].is_not_null())
print(f"Number of ORFs with -D: {df_D_FUN['ORF ID'].len()}")

Number of ORFs: 1333635 (73%)
Number of ORFs with -D: 1345873 (72%)

全ORFのうち73%程度はCOG,KEGG,PFAMの情報を持ち、データベースごとに見ると、全ORFに対してKEGGは52%程度、COGは78%程度、PFAMは43%程度の機能情報が期待できそうです。

Binのクオリティについても確認してみます。