SingularityでQIIME2の実行環境を作成する

はじめに

前回、解析環境の保守性の高さと、構築とHPCでコンテナファイルに持っていけば解析ができるような移植性の高さに魅力を感じてSingularityをインストールした。

https://zenn.dev/edna_startup/articles/438d5d96342c29

実際にSingularityでソフトウェアを実行する場合、コンテナファイルである.sifが必要となる。

今回は例として、QIIME2と周辺ツールを対象に、Singularityコンテナの定義ファイルである.defの作成から.sifをビルドするまでを記載する。

1. .defファイル

Singularity CEはカスタム色の強いコンテナ構築ツールで、defファイルにはベースとなるOSや、ソフトウェアのインストール、環境変数、メタデータに至る様々な情報を記載できる。

このdefファイルはheaderとsectionの２つのパートに別れている。

1.1 Header

Headerにはコンテナ構築に使用するOSの情報を記載する。

Bootstrapは、ヘッダーの最初に記載する必要がある。

既に構築されているイメージ集に関する情報といったことろ。

Bootstrapエージェントはいくつか種類があり、libraryやdocker、shubが良く使われている印象。

後はマシンに保存されているイメージを指定できるlocalimageやCentOS や Scientific Linux などの yumベースのシステムであるyumなどもある。

今回はBootstrapにdockerを指定して、Dockerレイヤーを使用した基本イメージの作成を行った。

Bootstrapはキーワードの指定が必須で、一番良く使うのはFrom。

簡素に書くなら以下のように記載する。

Bootstrap: docker
From: ubuntu:22.04

これでUbuntu 22.04のSingularityコンテナを作ることができるほか、CentOS, DebianなどのOSも指定することもできる。

これができると、特定のOSでしか稼働しないツール (CentOS7を要求するBCL Convertなど) を、手元の環境で実行できるのでこの状態でも十分使えるツールと言える。

1.2 Section

SectionはHeader以外の情報を定義する。%に続けてセクション名を記載し、以降にセクションにあった情報を記載する。

ちなみにBuild時は/bin/shが利用されこれにはオプションも指定可能。

また、defファイルに記載する各セクションの順序は重要ではないが、Build時には実行される順序がある。

以下の順序で実行される。

1. %setup
2. %files
3. %app*
4. %post
5. %test
6. %environment
7. %startscript
8. %runscript
9. %labels
10. %help

以降で今回使用するsectionについて説明する。

ちなみに、SingularityCE v4以降はテンプレートエンジンのような機能をサポートしており、defファイル内に{{ 変数名 }}と記載しておき、Build時に--build-arg 変数名="入れたい値"とすることで内部で値を受け取ることができる。

バージョンなどを都度変えたい場合などに使用するといい。

テンプレート機能については以下を参照されたい。
https://docs.sylabs.io/guides/4.1/user-guide/definition_files.html#templating-how-to-pass-values-into-definition-files

%post

このセクションでは、wgetを利用したファイルをダウンロードや、ソフトウェアのインストール、設定ファイルの編集、ディレクトリの作成などを実施する。

新たなPCを購入した時にツールをインストールするようなイメージでスクリプトを記載する。ここがうまく書ければ大体意図した形でBuildでき、コンテナファイルを作成できるはず。

注意点として、ホストの環境変数は基本的に使用できない。使用したい場合はSingularityCE v4から導入されたテンプレートを利用することで情報を渡す必要がある。

以下では、必要なソフトウェアをaptでインストールしつつ、miniconda3でcondaを使用できるようにして、qiime2をインストールしている。

%post
    # Update & Install software 
	apt update && apt upgrade -y 
	apt install -y \
        wget

    # make directory
    mkdir -p /opt/bin

    # Install miniconda
    wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
    bash Miniconda3-latest-Linux-x86_64.sh -b -p /opt/miniconda3

    . /opt/miniconda3/etc/profile.d/conda.sh
    
    conda activate base
    conda update --yes --all

    # set channel
    conda config --append channels conda-forge
    conda config --append channels bioconda
    
    # Install mamba
    conda install conda-forge::mamba

    # Install qiime2 amplicon
    wget https://data.qiime2.org/distro/amplicon/qiime2-amplicon-{{ qiime2_version }}-py38-linux-conda.yml
    mamba env create -n qiime2-amplicon-{{ qiime2_version }} --file qiime2-amplicon-{{ qiime2_version }}-py38-linux-conda.yml

以下では後に利用するyamlファイルなどをダウンロードする時に使用するwgetをインストールしている。

    # Update & Install software 
	apt update && apt upgrade -y 
	apt install -y \
        wget

以下で、ツールの置き場を作成。

    # make directory
    mkdir -p /opt

ここからはminiconda3をインストールして、condaを使えるようにしつつ、高速なインストーラーであるmamabaをインストールする。

    # Install miniconda
    wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
    bash Miniconda3-latest-Linux-x86_64.sh -b -p /opt/miniconda3

    . /opt/miniconda3/etc/profile.d/conda.sh
    
    # Activate
    conda activate base
    conda update --yes --all
    
    # Install mamba
    conda install conda-forge::mamba

続いて、qiime2 2024.2のツールのバージョン情報などが記載されたyamlファイルをダウンロードして、その情報をともに仮想環境を作成後にActivateする。

    # Install qiime2 amplicon 2024.02 
    wget https://data.qiime2.org/distro/amplicon/qiime2-amplicon-{{ qiime2_version }}-py38-linux-conda.yml

    mamba env create -n qiime2-amplicon-{{ qiime2_version }} --file qiime2-amplicon-{{ qiime2_version }}-py38-linux-conda.yml

    conda activate qiime2-amplicon-{{ qiime2_version }}

QIIME2には後発で開発されている様々なプラグインがある。基本的にpipかgithubリポジトリ上からインストールするものが多い模様。今回は、pipでの導入が推奨されていたのでそちらに従ってインストールした。

    # Install plugin
    pip install q2-greengenes2

最後にcacheをクリア

    # cache clear
    conda clean -a

これで%postの処理ステップは終了

%environment

このセクションには、実行時に設定される環境変数を定義する。注意点として、実行時にはコンテナで使用可能となるが、ビルド時には利用できない。ビルド中に同じ変数が必要な場合には、%postセクションに定義する必要がある。

先程%postで作成した仮想環境は常にActivateして使用するわけではないので、Activateした時と同様にPATHに仮想環境内の実行ファイルまでのパスを記載しておく。

また、常に明示的にするのであれば不要な気もするが、ロケール環境変数 LC_ALLにはC.UTF-8を指定した。

LC_ALLが設定されている場合は、環境変数 LANGもそちらに従い、設定がない場合には個別に割り当てられた値が利用されるらしい。

%environment
    # Setting variable for run
    export PATH=/opt:/opt/miniconda3/bin:/opt/miniconda3/envs/qiime2-amplicon-{{ qiime2_version }}/bin:${PATH}
    export LC_ALL=C.UTF-8

%labels

このセクションでは、コンテナ内の%labelファイルにメタデータを追加するために使用する。

%labels
    Author Author
    Build_date 2024.04.21
    Version v1.0

テキストの最初のスペースまでの部分がラベルの名前として解釈され、その後の部分がラベルの値として解釈される。

Build後に以下のコマンドを実行すると、イメージで使用可能なラベルを確認できる。

singularity inspect my_container.sif

Author: Author
Build_date: 2024.04.13
Version: v1.0
org.label-schema.build-arch: amd64
org.label-schema.build-date: Sunday_21_April_2024_15:35:52_JST
org.label-schema.schema-version: 1.0
org.label-schema.usage: /.singularity.d/runscript.help
org.label-schema.usage.singularity.deffile.bootstrap: docker
org.label-schema.usage.singularity.deffile.from: ubuntu:22.04
org.label-schema.usage.singularity.runscript.help: /.singularity.d/runscript.help
org.label-schema.usage.singularity.version: 4.1.2
org.opencontainers.image.ref.name: ubuntu
org.opencontainers.image.version: 22.04

%help

このセクションでは、singularity run-help <container>.sifを実行時に表示されるhelpメッセージを記載する。

%help
	This container for 16s amplicon sequence analysis using qiime2 version 2024.02.

次はコンテナをbuildする。

2. .defファイルからSingularityコンテナをBuildする

ここまで説明してきたdefファイルをまとめると以下の通り。

qiime2-2024-02.def

Bootstrap: docker
From: ubuntu:22.04

%post
    # Install 
	apt update && apt upgrade -y 
	apt install -y \
        wget

    # make directory
    mkdir -p /opt

    # Install miniconda
    wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
    bash Miniconda3-latest-Linux-x86_64.sh -b -p /opt/miniconda3

    . /opt/miniconda3/etc/profile.d/conda.sh
    
    # Activate
    conda activate base
    conda update --yes --all
    
    # Install mamba
    conda install conda-forge::mamba

    # Install qiime2 amplicon 2024.02 
    wget https://data.qiime2.org/distro/amplicon/qiime2-amplicon-{{ qiime2_version }}-py38-linux-conda.yml
    mamba env create -n qiime2-amplicon-{{ qiime2_version }} --file qiime2-amplicon-{{ qiime2_version }}-py38-linux-conda.yml
    conda activate qiime2-amplicon-{{ qiime2_version }}

    # Install plugin
    pip install q2-greengenes2

    # cache clear
    conda clean -a

%environment
    # Setting variable for run
    export PATH=/opt:/opt/miniconda3/bin:/opt/miniconda3/envs/qiime2-amplicon-{{ qiime2_version }}/bin:${PATH}
    export LC_ALL=C.UTF-8

%labels
    Author Author
    Build_date 2024.04.21
    Version v1.0

%help
	This container for 16s amplicon sequence analysis using qiime2 version 2024.02.

singularityのbuildコマンドでここまで作成したdefファイルからSingularityコンテナをビルドする。

singularity build --fakeroot --build-arg qiime2_version="2024.2" qiime2-2024-02.sif qiime2-2024-02.def

またよく使われる周辺ツールのdefファイルは以下の通り。後で使うのでこちらもBuildしておく。

qiime2_othertools.def

Bootstrap: docker
From: ubuntu:22.04

%post
    # Install 
	apt update && apt upgrade -y 
	apt install -y \
        wget

    # make directory
    mkdir -p /opt

    # Install miniconda
    wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
    bash Miniconda3-latest-Linux-x86_64.sh -b -p /opt/miniconda3

    . /opt/miniconda3/etc/profile.d/conda.sh

    # Activate    
    conda activate base
    conda update --yes --all

    # set channel
    conda config --append channels conda-forge
    conda config --append channels bioconda
    
    # Install mamba
    conda install conda-forge::mamba

    # Install informatics tools
    mamba create -n info-tools \
    python=3.10 \
    fastp=0.23.4 \
    multiqc=1.21 \
    seqkit=2.8.1 \
    parallel \
    aria2 \
    pigz \
    sra-tools=3.0 \
    parallel-fastq-dump \
    seqfu

%environment
    # Setting variable for run
    export PATH=/opt:/opt/miniconda3/bin:/opt/miniconda3/envs/info-tools/bin:${PATH}
    export LC_ALL=C.UTF-8

%labels
    Author Author
    Build_date 2024.04.21
    Version v1.0

%help
	This container for informatics tools with 16s amplicon sequence analysis using qiime2 version 2024.2.

singularity build --fakeroot qiime2_othertools.sif qiime2_othertools.def

3. せっかくなので作ったコンテナとQiime2を使って簡単なデータ解析

使用するのは以下の論文の16s rRNAのシーケンスデータ。

https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP367264&o=acc_s%3Aa

Metadataをダウンロードして、SraRunTable.txtのAssay TypeでAMPLICONのAccession Noを抽出する。

フィルタリングをかけた後のファイルはこちらからもダウンロードできるようにしておきました(SraRunTable.txt)。

https://github.com/NaokiShibata/Demodata/tree/main/singularity_qiime2

3.1 配列データの取得

SraRunTable.txtのうち必要な情報をフィルタリングし、ファイルの1列目の2行目以降を抽出して、arrayとして`data_arrayに代入します。

Accession Numberは全部で298サンプルありそのまま実行するとsraファイルだけで2GB程度になる。もしStrageを節約したい場合は、shufでランダムに配列情報を取得して新しいarrayを作成する。

# Accession Numnberのarray化
mapfile -t data_array < <(awk -F '[,]' 'NR > 1 {print $1}' "SraRunTable.txt")

# 元の配列からランダムに10個の要素を選択して新しい配列に代入する場合は以下を使う
# random_data_array=($(shuf -e "${data_array[@]}" -n 10))

コンテナ内にダウンロードしておいたaria2を使用して、sraファイルをダウンロードする。

# SRA用のフォルダ作成
mkdir 00sra

# sraファイルのダウンロード
for i in ${data_array[@]}; do
    singularity exec qiime2_othertools.sif aria2c \
    -x 16 \
    -ctrue \
    --dir=00sra \
    https://sra-pub-run-odp.s3.amazonaws.com/sra/${i}/${i}
done

続いて、こちらもコンテナ内にダウンロードしておいたparallel-fastq-dumpとpigzを使用して、sra形式からPaierd Endのfastq.gzファイルを生成する。

# フォルダの作成
mkdir 00fq

# sra形式 -> fastq -> fastq.gz
printf "%s\n" "${data_array[@]}" | xargs -I {} -P 4 singularity exec qiime2_othertools.sif sh -c "
parallel-fastq-dump \
--sra-id 00sra/{} \
--threads 8 \
--outdir 00fq/ \
--split-files;

pigz -p 8 00fq/{}*fastq
"

Accession Numberをxargsに渡し、並列でコマンドを実行している。

• printf "%s\n" "${data_array[@]}"
  • 配列の要素を改行で区切って出力します。これにより、各行が xargs によって個別のコマンドとして扱われる
• xargs -I {} -P 4
  • -I {} は置換文字列 {} を使用して、入力行をコマンドの任意の位置に挿入できるようにする
  • -P 4 は同時に実行するプロセスの最大数を4に設定します。この数はシステムのリソースに応じて調整可能
• singularity exec qiime2_othertools.sif sh -c "...":
  • singularity exec でコンテナ内部でシェルコマンドを実行します。sh -c は複数のコマンドを含む一連のシェルコマンドを扱うために使用される

これでfastq.gzファイルが生成できた。

3.2 Manifestファイルの作成とQiime2へのデータのインポート

以降ちょっと古いですが、以下の記事を参考にすすめる。

https://edna-blog.com/analyze/qiime2/

まずはじめにseqfuでManifestファイルを作成。

singularity exec qiime2_othertools.sif \
seqfu metadata \
--format manifest \
-1 "_1" \
-2 "_2" \
--threads 32 \
00fq > Manifest.tsv

続いて、qiime tools importで解析対象のfastqファイル群をqza形式に変換。

--input-pathにはさっき作成した、Pathが記載されているManifest.tsvを指定。--typeと--iput-formatはこちらを参考にして選択。

mkdir -p 01output/qza 01output/qzv log

singularity exec qiime2-2024-02.sif \
tools import \
--type 'SampleData[PairedEndSequencesWithQuality]' \
--input-format PairedEndFastqManifestPhred33V2 \
--input-path Manifest.tsv \
--output-path 01output/qza/00demux.qza \
2>&1 | tee log/qiime_demux.log

3.3 DADA2によるデノイジング、エラー配列の修正

論文にはMiseqで250PEのシーケンスを下との記載があるが、編集前のSraRunTable.txtには600bpのアンプリコンをシーケンス(300PE)したとの情報があります。seqkitで確認しても、R1, R2ともに300bpあったので300PEで間違いないでしょう。

DADA2の処理についても論文に記載のパラメーターに準拠して解析を実行しており、DADA2のpooling methodにはpseudo-poolingを使用したと記載があったので、--p-pooling-method puseudoとしています。

--p-pooling-methodに選択できる、'independent', 'pseudo'については、実行時間が増えるもののpuseudoのほうが感度は高い模様。

また、過去のissueで議論されているように、pooling-methodをpuseudoにした場合、chimera-methodもpooldにすることが望ましいようで、将来的に実装されるとのことみたいだが今の所実装されておらず、デフォルトのconsensusが指定されている。今回はそのままとした。

singularity exec qiime2-2024-02.sif \
qiime dada2 denoise-paired \
--i-demultiplexed-seqs ./01output/qza/00demux.qza \
--p-trim-left-f 0 \
--p-trim-left-r 0 \
--p-trunc-len-f 250 \
--p-trunc-len-r 180 \
--p-pooling-method pseudo \
--o-representative-sequences ./01output/qza/02rep-seqs-dada2_f250_r180.qza \
--o-table ./01output/qza/02table-dada2_f250_r180.qza \
--o-denoising-stats ./01output/qza/02stats-dada2_f250_r180.qza \
--p-n-threads 0 \
--verbose 2>&1 | tee log/qiime_dada2.log

3.4 read filtering

論文の情報に従いリード数でテーブルデータに対してフィルタリングを実施する。

singularity exec qiime2-2024-02.sif \
qiime feature-table filter-features \
--i-table 01output/qza/02table-dada2_f250_r180.qza \
--p-min-frequency 2500 \
--o-filtered-table 01output/qza/03frequency-filtered-table.qza

次に、フィルタリング後に残った情報に対応する配列データを抽出する。この作業をしないとテーブルデータと齟齬が出る。

singularity exec qiime2-2024-02.sif \
qiime feature-table filter-seqs \
--i-data 01output/qza/02rep-seqs-dada2_f250_r180.qza \
--i-table 01output/qza/03frequency-filtered-table.qza \
--o-filtered-data 01output/qza/03filtered-rep-seqs.qza

フィルタリング後のデータをqiimeのアーティファクトファイルから抽出したい場合は、tools exportで取得することができる。

singularity exec qiime2-2024-02.sif \
qiime tools export \
--input-path 01output/qza/02rep-seqs-dada2_f250_r180.qza \
--output-path 01output/other1/

singularity exec qiime2-2024-02.sif \
qiime tools export \
--input-path 01output/qza/03filtered-rep-seqs.qza \
--output-path 01output/other2/

seqkit statsで配列データのサマリーを取得

singularity exec qiime2_othertools.sif \
seqkit stats -a -T \
01output/other1/dna-sequences.fasta > 01output/other1/dna-sequences.tsv

singularity exec qiime2_othertools.sif \
seqkit stats -a -T \
01output/other2/dna-sequences.fasta > 01output/other2/dna-sequences.tsv

	format	type	num_seqs	sum_len	min_len	avg_len	max_len	Q1	Q2	Q3	sum_gap	N50	N50_num	Q20(%)	Q30(%)	AvgQual	GC(%)
フィルタリング前	FASTA	DNA	1921	485669	250	252.8	255	253.0	253.0	253.0	0	253	3	0.00	0.00	0.00	52.74
フィルタリング後	FASTA	DNA	225	56894	252	252.9	254	253.0	253.0	253.0	0	253	2	0.00	0.00	0.00	52.83

だいぶ配列数が減って変わりましたね。

3.5 Greengenes2による系統分類

少し変わり種として、2023年に更新されたGreengenes2のデータベースを使って系統分類してみる。

https://github.com/biocore/q2-greengenes2/

Greengenes2は初期のGreengenesから設計が大きく更新されており、全ゲノムのデータからマーカー遺伝子に基づく系統樹を作成して分類を決めるというデータベースの構築方法をとっている。これにより、16s rRNAとショットガンの両方で直接統合可能な結果を得ることができるらしい。

https://www.nature.com/articles/s41587-023-01845-1

QIIME2での使用においては、一般的に用いられるq2-feature-classifierプラグインを使用した系統分類の他にq2-greengenes2が利用できるので、以下の流れで実行しようと思ったがエラーが出てうまく実行できなかった。

Error q2-greengenes2

データベースはこちらからダウンロードした : 2022.10.taxonomy.asv.nwk.qza

mkdir 00ref
wget https://ftp.microbio.me/greengenes_release/current/2022.10.taxonomy.asv.nwk.qza \
--no-check-certificate \
-o 00ref/2022.10.taxonomy.asv.nwk.qza

singularity exec qiime2-2024-02.sif \
qiime greengenes2 filter-features \
--i-reference 00ref/2022.10.taxonomy.asv.nwk.qza \
--i-feature-table 01output/qza/02table-dada2_f250_r180.qza \
--o-filtered-feature-table 01output/qza/04table-dada2_f250_r180_gg2.qza \
--verbose 2>&1 | tee log/qiime_filter-features_gg2.log

singularity exec qiime2-2024-02.sif \
qiime greengenes2 taxonomy-from-table \
--i-reference-taxonomy 00ref/2022.10.taxonomy.asv.nwk.qza \
--i-table 01output/qza/04table-dada2_f250_r180_gg2.qza \
--o-classification 01output/qza/04gg2_taxonomy.qza \
--verbose 2>&1 | tee log/qiime_taxonomy-from-table_gg2.log

Issueもあって解決法が提案されているがあまり時間をかけたくないので、以下のデータベースを使って、q2-feature-classifierで系統分類を実施した。

データベースはこちらからダウンロードした : 2022.10.backbone.full-length.nb.qza

データベースのダウンロード

# Reference用のフォルダの作成と移動
mkdir 00ref && cd 00ref
# ダウンロード
wget https://ftp.microbio.me/greengenes_release/current/2022.10.backbone.full-length.nb.qza 
# もとのフォルダへ移動
cd ..

データベースが用意できればfeature-classifier classify-sklearnで以下のように系統分類を実行する。

singularity exec qiime2-2024-02.sif \
qiime feature-classifier classify-sklearn \
--i-classifier 00ref/2022.10.backbone.full-length.nb.qza \
--i-reads 01output/qza/03filtered-rep-seqs.qza \
--o-classification 01output/qza/04taxonomy_gg2.qza \
--p-n-jobs 10 \
--verbose 2>&1 | tee log/qiime_taxonomy-from-table_gg2-sklearn.log 

10分程度で終了

3.6 Barplotの描画

SraRunTable.txtの情報を整理して適切な形にしたmetadata.tsvを使って、Barplotを描画。

metadata.tsvは以下を参照。

https://github.com/NaokiShibata/Demodata/tree/main/singularity_qiime2

singularity exec qiime2-2024-02.sif \
qiime taxa barplot \
--i-table 01output/qza/03frequency-filtered-table.qza \
--i-taxonomy 01output/qza/04taxonomy_gg2.qza \
--m-metadata-file metadata.tsv \
--o-visualization 01output/qzv/gg2-taxa-bar-plots.qzv

Level4: Orderでsample typeごとに並び替えたbarplot

3.7 α多様性の視覚化 : rarefaction curveの描画

いくつかの多様度指数を算出後、まとめてrarefactionとして描画します。

まずは、phylogeny align-to-tree-mafft-fasttreeで系統樹を作成。

singularity exec qiime2-2024-02.sif \
qiime phylogeny align-to-tree-mafft-fasttree \
--i-sequences 01output/qza/03filtered-rep-seqs.qza \
--output-dir 01output/other/phylogeny-align-to-tree-mafft-fasttree

次に、有根系統樹を使っていくつかの多様度指数のrarefraction curveを作成・描画する。

singularity exec qiime2-2024-02.sif \
qiime diversity alpha-rarefaction \
--i-table 01output/qza/03frequency-filtered-table.qza \
--i-phylogeny 01output/other/phylogeny-align-to-tree-mafft-fasttree/rooted_tree.qza \
--m-metadata-file metadata.tsv \
--p-max-depth 20000 \
--o-visualization  01output/qzv/05alpha-rarefaction-plot.qzv \
--verbose 2>&1 | tee log/qiime_rarefaction.log

observed_features, sample_setの表示

3.8 α多様性の視覚化 : カテゴリごとのα多様性の比較

Metadata.tsvで規定しているカテゴリごとのα多様性を見比べることができる。

singularity exec qiime2-2024-02.sif \
qiime diversity core-metrics-phylogenetic \
--i-phylogeny 01output/other/phylogeny-align-to-tree-mafft-fasttree/rooted_tree.qza \
--i-table 01output/qza/03frequency-filtered-table.qza \
--p-sampling-depth 10000 \
--p-n-jobs-or-threads 32 \
--m-metadata-file metadata.tsv \
--output-dir 01output/other/diversity-core-metrics-phylogenetic
--verbose 2>&1 | tee log/qiime_core-metrics-phylogenetic.log

singularity exec qiime2-2024-02.sif \
qiime diversity alpha-group-significance \
--i-alpha-diversity 01output/other/diversity-core-metrics-phylogenetic/observed_features_vector.qza \
--m-metadata-file metadata.tsv \
--o-visualization 01output/qzv/06alpha-group-sig-obs-feats.qzv \
--verbose 2>&1  | tee tee log/qiime_alpha-group-significance.log

カテゴリごとに確認されたASV数のbox plot

Kruskal-Wallisによる差の検定の結果

他にもたくさんプラグインがあり多様な解析ができる。

3.9 β多様性の視覚化 : PCoA (主座標分析)

PCoAはデータセット内のオブジェクト間の距離を視覚化するために用いられる。一般的にBray-Cutis, Jaccardなどの距離指標を用いた行列から算出・描画する。

以下のスクリプトでは、--p-metricで指定できる値をloopで指定してすべての結果が出力されるようにした。

# 指定できるmetricsについてすべてのplotを生成
metrics_array=(aitchison braycurtis canberra canberra_adkins chebyshev cityblock correlation cosine dice euclidean generalized_unifrac hamming jaccard jensenshannon kulsinski matching minkowski rogerstanimoto russellrao seuclidean sokalmichener sokalsneath sqeuclidean unweighted_unifrac weighted_normalized_unifrac weighted_unifrac yule)

for met in ${metrics_array[@]};do
    singularity exec qiime2-2024-02.sif \
    qiime diversity beta-rarefaction \
    --i-table 01output/qza/03frequency-filtered-table.qza \
    --p-metric ${met} \
    --p-clustering-method nj \
    --p-sampling-depth 15000 \
    --m-metadata-file metadata.tsv \
    --o-visualization 01output/qzv/07${met}-rarefaction-plot.qzv \
    --verbose 2>&1 | tee log/qiime_beta-rarefaction_${met}.log
done

各距離指標に基づく3Dのこんな図が得られる

3.10 β多様性解析 : PERMANOVA

PERMANOVAは1つのグループ内のサンプル間の距離 (グループ内距離) が、別のグループ内のサンプルまでの距離 (グループ間距離)と異なるという仮説をテストする。

comp_set=(host_sex sample_type)

for comp in ${comp_set[@]}; do
    singularity exec qiime2-2024-02.sif \
    qiime diversity beta-group-significance \
    --i-distance-matrix 01output/other/diversity-core-metrics-phylogenetic/unweighted_unifrac_distance_matrix.qza \
    --m-metadata-file metadata.tsv \
    --m-metadata-column ${comp} \
    --p-pairwise \
    --p-method 'permanova' \
    --o-visualization 01output/qzv/09uw_unifrac-${comp}-significance.qzv \
    --verbose
done

for comp in ${comp_set[@]}; do
    singularity exec qiime2-2024-02.sif \
    qiime diversity beta-group-significance \
    --i-distance-matrix 01output/other/diversity-core-metrics-phylogenetic/weighted_unifrac_distance_matrix.qza \
    --m-metadata-file metadata.tsv \
    --m-metadata-column ${comp} \
    --p-pairwise \
    --p-method 'permanova' \
    --o-visualization 01output/qzv/09w_unifrac-${comp}-significance.qzv \
    --verbose
done

PERMANOVAとセットでよく解析されるのがPERMDISP。こちらはグループ間のβ多様性の等分散性を比較検定できる。

comp_set=(host_sex sample_type)

for comp in ${comp_set[@]}; do
    singularity exec qiime2-2024-02.sif \
    qiime diversity beta-group-significance \
    --i-distance-matrix 01output/other/diversity-core-metrics-phylogenetic/unweighted_unifrac_distance_matrix.qza \
    --m-metadata-file metadata.tsv \
    --m-metadata-column ${comp} \
    --p-pairwise \
    --p-method 'permdisp' \
    --o-visualization 01output/qzv/09uw_unifrac-${comp}-significance_permdisp.qzv \
    --verbose
done

for comp in ${comp_set[@]}; do
    singularity exec qiime2-2024-02.sif \
    qiime diversity beta-group-significance \
    --i-distance-matrix 01output/other/diversity-core-metrics-phylogenetic/weighted_unifrac_distance_matrix.qza \
    --m-metadata-file metadata.tsv \
    --m-metadata-column ${comp} \
    --p-pairwise \
    --p-method 'permdisp' \
    --o-visualization 01output/qzv/09w_unifrac-${comp}-significance_permdisp.qzv \
    --verbose
done

左: PERMDISP, 右: PERMANOVA, カテゴリ: Sample_type

PERMANOVAとPERMDISPはセットで使われることが多い。その必要性と考察の仕方については、以下の和文を参考に参考にされたい。

https://www.jstage.jst.go.jp/article/seitai/66/1/66_1/_pdf

さいごに

これでなんとなくだがsingularity コンテナが作成できるようになった。defファイルをgithubで管理したりすることで保守的なversion管理ができ、作成したsifファイルを別のPCやスパコンへ持って行くことですぐに解析ができるようになる。はず。

qiime2も年々新しくなって知らないpluginが沢山追加されているので取りこぼさないようにしないといけない。とりあえずはじめは以下のプロトコル論文などを参考に進めるのがいいと思う。

https://currentprotocols.onlinelibrary.wiley.com/doi/10.1002/cpbi.100

次は、workflow言語としてNextflowでも触ろうかな。

出典

• https://docs.sylabs.io/guides/4.1/user-guide/definition_files.html#templating-how-to-pass-values-into-definition-files
• https://geniac.readthedocs.io/en/latest/conda.html
• https://telatin.github.io/singularities/containers/qiime.html
• https://github.com/apptainer/singularity
• https://github.com/sylabs/singularity
• https://github.com/jianhong/16S_pipeline/blob/main/singularity.def
• https://bioinformatics.ccr.cancer.gov/docs/qiime2/Lesson5/
• https://bioinformatics.ccr.cancer.gov/docs/qiime2/Lesson6/
• https://currentprotocols.onlinelibrary.wiley.com/doi/10.1002/cpbi.100
• https://www.jstage.jst.go.jp/article/seitai/66/1/66_1/_pdf