macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

メタゲノムのリードの発生からbinningまで自動でシミュレーションする MAGICIAN

 

 シーケンスリードからメタゲノムアセンブリゲノム(MAGs)を回収することで、微生物群集とその構成員に関するさらなる洞察が可能になり、場合によっては単一分離ゲノム用に設計されたツールでそのような配列を解析することもできる。結果の質は配列の質に依存するため、単一分離ゲノム用のツールの性能をMAG上で事前にテストする必要がある。バイオインフォマティクスレバレッジして、この目的のために、構成が既知の様々な合成テストセットを迅速に作成することができる。

本著者らは、MAGのシミュレーションのための柔軟でユーザーフレンドリーなパイプラインであるMAGICIANを発表する。MAGICIANは、合成メタゲノムシミュレータとメタゲノム解析およびビニングパイプラインを組み合わせ、ユーザーが提供した入力ゲノムに基づいてMAGをシミュレートする。MAGICIANを使用することで、カバレッジの深さがごくわずか(1%)変化するだけでも、ゲノムを回収できるかどうかに大きく影響することが分かった。また、MAGICIANの現在のデフォルトパイプラインで得られたゲノムを抗菌薬耐性遺伝子同定ツールResFinderで解析する際の適合性を評価することで、シミュレートされたMAGの利用を実証した。

MAGICIANを用いることで、一般的に高品質でありながら、実データで発生する問題を反映したMAGをシミュレートすることができ、現実的なベストケースデータを提供することができる。これらのゲノムをResFinderで解析した結果を評価したところ、最も妥当なルッキングのリスクが明らかになった。

 

インストール

Github

https://github.com/KatSteinke/magician

#本体
git clone https://github.com/KatSteinke/magician.git
cd magician/
mamba env create --file requirements.yml
conda activate magician_base

#CAMISIM
#CAMISIMのcustom error profilesが必要。CAMISIMのfolkをクローンしておく
git clone https://github.com/KatSteinke/CAMISIM

default_config.ymlを開いて、クローンしたCAMISIMのパスを指定する

vi default_config.yml

python run_magician.py -h

usage: run_magician.py [-h] [--target TARGET]

                       [--profile_type {mbarc,hi,mi,hi150,own}]

                       [--profile_name PROFILE_NAME]

                       [--profile_readlength PROFILE_READLENGTH]

                       [--insert_size INSERT_SIZE] [--cluster CLUSTER]

                       [--cores CORES] [--config_file CONFIG_FILE]

                       [--snake_flags [SNAKE_FLAGS [SNAKE_FLAGS ...]]]

                       community_file

 

Run MAGICIAN to simulate MAGs for a specified community or set of communities.

Run without arguments for a test run.

 

positional arguments:

  community_file        File with paths to source genomes, sequence type

                        (plasmid/chromosome) and their desired relative copy

                        number in each community to simulate (one column with

                        organisms' copy numbers per community)

 

optional arguments:

  -h, --help            show this help message and exit

  --target TARGET       Desired output file or rule (default: MAGs, statistics

                        and summary files for all communities)

  --profile_type {mbarc,hi,mi,hi150,own}

                        Type of ART error profile to use for CAMISIM: mbarc,

                        hi, mi, hi150, own (default: mbarc)

  --profile_name PROFILE_NAME

                        Base name of custom error profile, if given (name of

                        files without '[1/2].txt'); required with 'own' error

                        profile

  --profile_readlength PROFILE_READLENGTH

                        Read length of custom error profile; required with

                        'own' error profile

  --insert_size INSERT_SIZE

                        Mean insert size for read simulation (default: 270)

  --cluster CLUSTER     For use with snakemake's cluster mode; supply command

                        for submitting jobs as you would with snakemake.

  --cores CORES         Amount of cores Snakemake should use (default: 6)

  --config_file CONFIG_FILE

                        Config file for run

  --snake_flags [SNAKE_FLAGS [SNAKE_FLAGS ...]]

                        Flags to be passed to snakemake, enclosed in quotes

 

 

テストラン

オプションなしで実行すると、テストランできる。"y"を選択。

python run_magician.py community.tsv

snakemakeのランのためにcondaを使って仮想環境が作られるが、依存関係が解決できず失敗した。

 

 

実行方法

MAGICIANの実行には、シミュレーションの対象となる生物のゲノム配列(genbankまたはfasta形式)と、タブ区切りのサンプル分布ファイルが必要。タブ区切りファイルの最初の列はゲノムファイルのパス、2番目は配列の種類(chromosomeまたはプラスミド)、それに続く列はコミュニティの構成名とコミュニティ内の配列の相対的な存在量。

https://github.com/KatSteinke/magician/blob/master/test/data/test_genomes/sample_distributions.tsv

 

TSVファイルを指定して実行する。また、--snake_flagsでSnakemakeに渡すフラグを二重引用符で囲んで指定する必要がある。ドライランには"-n "を使用する。

python run_magician.py community.tsv

最後に以下のようなサマリーレポートもできる。

https://github.com/KatSteinke/magician/tree/master/test/data/sample_summaries

 

 

論文より

  • MAGICIANは、ユーザーが定義した合成コミュニティからのメタゲノムリードとMAGをインシリコでシミュレートし、これらのコミュニティの組成を簡単に変化させ、大量の実験をシミュレートすることができる。
  • 機能を実証するために、関連する生物やプラスミドを導入し、シークエンシングパラメーターを変化させた。その結果、プラスミドは一般的に回収できず、近縁の生物の再構築は困難であった。また、ほとんどのケースで近縁な生物の配列を含むビンは非常に断片化された。

引用

MAGICIAN: MAG simulation for investigating criteria for bioinformatic analysis
Kat Steinke, Sünje J. Pamp & Patrick Munk 
BMC Genomics volume 25, Article number: 55 (2024)

 

関連

https://kazumaxneo.hatenablog.com/entry/2018/08/10/123830

 

 

 

微生物の機能をGO termの形で予測する DeepGOMeta

 

 微生物サンプルの解析は、その多様性と複雑性のために、依然として計算上困難である。ロバストなde novoタンパク質機能予測法の欠如は、これらのサンプルから機能的洞察を導き出すことの難しさを悪化させている。相同性や配列の類似性に依存する従来の予測手法では、新規タンパク質やホモログが知られていないタンパク質の機能を予測できないことが多い。さらに、これらの手法のほとんどは、主に真核生物のデータに対して学習されたものであり、微生物のデータセットに対する評価や適用は行われていない。本研究では、微生物に関連するデータセットで学習させた、ジーオントロジー(GO)の語彙としてのタンパク質機能予測用に設計された深層学習モデルDeepGOMetaを紹介する。このモデルは、新しい評価戦略を用いて検証され、多様な微生物データセットに適用されている。データとコードはhttps://github.com/bio-ontology-research-group/deepgometaで利用できる。

 

DeepGOMetaは微生物種(原核生物古細菌、ウイルス)に属するUniProtKB/Swiss-Prot Knowledgebaseタンパク質でトレーニング、テスト、評価されている。DeepGOMetaは、GO termの形で機能予測を提供する。

インストール

推奨されている通り、dockerイメージを使ってテストした。使用前にdockerがGPUを認識できるようにしておく必要がある(*1)(os: ubuntu20)。

  • The code was developed and tested using python 3.10.
  • You'll need around 30Gb storage and a GPU with >16Gb memory (or you can use CPU)

Github

#ここではdockerを使用(dockerhub)
docker pull coolmaksat/deepgometa

> docker run --rm coolmaksat/deepgometa python predict.py --help

Usage: predict.py [OPTIONS]

 

Options:

  -if, --in-file TEXT        Input FASTA file  [required]

  -dr, --data-root TEXT      Prediction model

  -t, --threshold FLOAT      Prediction threshold

  -bs, --batch-size INTEGER  Batch size for prediction model

  -d, --device TEXT          Device

  --help                     Show this message and exit.

 

 

モデルのダウンロード(800MBほど)

wget https://deepgo.cbrc.kaust.edu.sa/data/deepgometa/data.tar.gz
unzip data.tar.gz

 

実行方法

タンパク質のfastaファイルを指定する。解凍したdata/をdockerと共有するため、解凍したdata/の1つ上に移動して実行する。下のコマンドではexample.faを解凍したdata/に配置している。

docker run --rm --gpus all -v $(pwd)/data:/workspace/deepgometa/data coolmaksat/deepgometa python predict.py -if data/example.fa

実行すると、初めにESM-2のモデルがダウンロードされる。細菌株のproteome 4000配列のアノテーションに40分ほどかかった(GPU: RTX3090)。

 

example.faaを指定した場合、example/にGO termの3つのカテゴリに相当する、example_preds_mf.tsv.gz、example_preds_cc.tsv.gz、example_preds_bp.tsv.gzが出力される。

gless example_preds_mf.tsv.gz

タンパク質ごとに複数のGO termが1行に1つずつ記載されている。多くの配列に複数のGO termがついている。出力について、レポジトリには説明が見当たらなかった。

 

  • DeepGOMetaは、Dockerイメージを使用してNextflowワークフローとして実行することもできる。
  • (論文より)MGnifyタンパク質データベースから2,000個のタンパク質をアノテーションする際に、既存のPfamアノテーションを持つタンパク質は567個に過ぎなかったのに対し、DeepGOMetaは2,000個のタンパク質全てのアノテーションに成功した。一方、このサブセットにはPfamのアノテーションがないため、予測の機能的関連性と精度を直接検証する上では課題がある。

引用

DeepGOMeta: Predicting functions for microbes
Rund Tawfiq,  Kexin Niu,  Robert Hoehndorf,  Maxat Kulmanov

bioRxiv, Posted January 31, 2024.

 

*1

こちらを参考にしました。

https://zenn.dev/derbuihan/articles/a1e636d29e1b51

ロングリードトランスクリプトームの高効率なクラスタリングを行う geluster

 

 ロングリードRNAシーケンス技術の進歩は、トランスクリプトーム解析に明るい未来をもたらした。ロングリードをその起源遺伝子ファミリーにしたがってクラスタリングすることは非常に重要である。しかし、既存のde novoクラスタリングアルゴリズムは、膨大な計算資源を必要とする。

ロングRNA-seqリードをクラスタリングする新しいアルゴリズムGeLusterを開発した。1つのシミュレーションデータセットと9つの実データセットでのテストで、GeLusterは優れた性能を示した。テストしたNanoporeデータセットでは、2番目に高速な手法の2.9~17.5倍の速度で動作し、メモリ消費量は7分の1以下でありながら、高いクラスタリング精度を達成した。また、PacBioデータでもGeLusterは同様の性能を示した。GeLusterは将来の大規模トランスクリプトーム研究の舞台に対応する。

 

インストール

minimap2とsamtoolsが必要。下では2つのツールの導入については省略している。

依存

  • g++ with support for C++11 (e.g. 9.4.0)
  • minimap2 (In the GeLuster paper we tested with version 2.24-r1122)
  • samtools (In the GeLuster paper we tested with version 1.10)

Github

https://github.com/yutingsdu/GeLuster

git clone https://github.com/yutingsdu/GeLuster.git
cd GeLuster/src/
make -j8
cd ../

> ../GeLuster  

[Error] : Reads file is not provided!

    

===========================================================================

 

GeLuster usage:

 

** Required **

 

--reads/-r <string>     : path to the read file

 

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

 

** Options **

 

--help/-h          : Output GeLuster Help Information.

 

--version/-v          : Print current version of GeLuster.

 

--iteration/-i <int>      : Number of GeLuster iterations ([3,9], default: 3).

 

--seqType/-s <string>      : 'cDNA' for ONT cDNA data, 'dRNA' for ONT direct RNA data, or 'PacBio' for pacbio data (default:cDNA).

 

--rform/-f <string>      : 'fq' for FASTQ format reads, 'fa' for FASTA format reads (default: fq).

 

--threads/-t <int>       : Number of threads to be used (default: 10).

 

--multi/-M           : To generate a proxy of gene expression matrix for multiple RNA-seq samples. Input files should be separated by commas.

 

--output_dir/-o <string>  : Output path, default: geluster_outdir.

 

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

 

** Typical commands **

 

A typical GeLuster command might be:

 

  GeLuster -r reads.fastq -f fq -s cDNA -o geluster_outdir 

 

===========================================================================

 

 

 

テストラン

cd GeLuster/sample_test/
./run_me.sh

出力

geluster_outdir/



実行方法

ロングリードを指定する。

GeLuster -r Test.fastq -s cDNA -o out
  • -s    'cDNA' for ONT cDNA data, 'dRNA' for ONT direct RNA data, or 'PacBio' for pacbio data (default:cDNA).
  • -f    'fq' for FASTQ format reads, 'fa' for FASTA format reads (default: fq).

 

出力TSVは1リードにつき1行のテキストファイルとなっている。各行には、リードの名前とそのリードが属するクラスタがプリントされる、

 

サブディレクトリfastq_filesにはクラスタ毎にfastqが保存される。

 

複数のサンプル(細胞)から得られたRNA-seqリードを入力とし、入力されたサンプル(細胞)について遺伝子発現行列を生成することもできる。

GeLuster -r reads1.fastq,reads2.fastq,reads3.fastq -f fq -s cDNA --multi -o geluster_outdir
  • --multi/-M    To generate a proxy of gene expression matrix for multiple RNA-seq samples. Input files should be separated by commas.

出力行列の各項目は、与えられたサンプル(細胞)における特定のクラスタ(遺伝子)のリード数(発現レベル)を表す。

 

レポジトリより

  • GeLusterは、参照情報を使わずにde novoで発現行列を生成する。このアプローチは、参照情報がない生物種でも機能すると期待される。

引用

Highly efficient clustering of long-read transcriptomic data with geluster 
Junchi Ma, Xiaoyu Zhao, Enfeng Qi, Renmin Han, Ting Yu, Guojun Li Author Notes
Bioinformatics, Published: 03 February 2024

ロングリードのハイブリッドエラー訂正を行う HERO

 

一般的に優れているが、次世代シーケンシング(NGS)リードを用いた第3世代シーケンシング(TGS)リードのエラーを修正するハイブリッドアプローチは、ハプロタイプ特異的バリアントを、倍数体サンプルや混合サンプルのエラーと取り違える。HEROは、NGSリードとTGSリードの特定の強みに最適に対応するために、de Bruijnグラフとオーバーラップグラフの両方を利用する最初の「ハイブリッド-ハイブリッド」アプローチとして提案する。広範なベンチマーク実験により、HEROはindelとミスマッチのエラー率を平均65% (27~95%)および20% (4~61%)、ゲノムアセンブリの前にHEROを使用することで、関連するカテゴリーの大部分でアセンブリが大幅に改善される。

 

HEROは2つのステージから構成されている。

 

HEROはHaplotype-aware ERRor cOrrectionを意味する。

インストール

依存

注;HERROはGPUのVRAMをかなり使う。推奨80GBかつ複数GPUとなっている(ヒトゲノムなど)。

 

説明されている通り、レポジトリをcloneして依存するツールを導入(前処理のステップ1と2用)、それからsingularityイメージをビルドした。最後にモデルをダウンロードして実行した(古いsingularityだと動かないので注意。v3.9.5を使用した)。

Github

https://github.com/lbcb-sci/herro

#1
git clone https://github.com/dominikstanojevic/herro.git
cd herro
mamba env create --file scripts/herro-env.yml

#2  singularityイメージのビルド
sudo singularity build herro.sif herro-singularity.def
#もしくはビルド済みのイメージをダウンロード(3.5GB)
wget http://complex.zesoi.fer.hr/data/downloads/herro.sif

> scripts/preprocess.sh

Please place porechop_with_split.sh and no_split.sh in the same directory as this script.

This script requires 4 arguments:

1. The input sequence file. e.g. input.fastq

2. The output prefix. e.g. 'preprocessed' or 'output_dir/preprocessed'

3. The number of threads to be used.

4. The number of parts to split the inputs into for porechop (since RAM usage may be high).

 

> scripts/create_batched_alignments.sh

Please place batch.py in the same directory as this script.

This script requires 4 arguments:

1. The path to the preprocessed reads.

2. The path to the read ids of these reads e.g. from seqkit seq -n -i.

3. The number of threads to be used.

4. The directory to output the batches of alignments.

 

> singularity run --nv herro.sif inference -h

$ singularity run --nv herro.sif inference -h

Subcommand used for error-correcting reads

 

Usage: herro inference [OPTIONS] -m <MODEL> -b <BATCH_SIZE> <READS> <OUTPUT>

 

Arguments:

  <READS>   Path to the fastq reads (can be gzipped)

  <OUTPUT>  Path to the corrected reads

 

Options:

      --read-alns <READ_ALNS>    Path to the folder containing *.oec.zst alignments

      --write-alns <WRITE_ALNS>  Path to the folder where *.oec.zst alignments will be saved

  -w <WINDOW_SIZE>               Size of the window used for target chunking (default 4096) [default: 4096]

  -t <FEAT_GEN_THREADS>          Number of feature generation threads per device (default 1) [default: 1]

  -m <MODEL>                     Path to the model file

  -d <DEVICES>                   List of cuda devices in format d0,d1... (e.g 0,1,3) (default 0) [default: 0]

  -b <BATCH_SIZE>                Batch size per device. B=64 recommended for 40 GB GPU cards.

  -h, --help                     Print help

 

 

 

モデルのダウンロード

wget http://complex.zesoi.fer.hr/data/downloads/model_v0.1.pt

 

実行方法

ステップ1、2はcondaで作った環境で実行する前処理のステップ(HEROは使わない)。ステップ3のherro inferenceコマンドはsingularityイメージを使用する。

 

1、Preprocess reads

PorechopによるONTリードの前処理。10スレッド指定。最後にメモリ使用量を削減するために一過的にリードを分割して扱うための分割数を指定する。分割が不要な場合は、メモリ使用量が増えるが<parts_to_split_job_into>を1に設定する。ONTデータが巨大なら分割が推奨される。

conda activate herro
scripts/preprocess.sh input_fastq out_prefix 10 <parts_to_split_job_into>

Dorado v0.5では、アダプタートリミングが追加されたため、Porechopやduplexツールを使用したアダプタートリミングや分割はおそらく将来不要になる(レポジトリより)。

outprefix.fastq.gzが得られる。

 

2、minimap2 alignment and batching

1の出力のロングリードを指定する。reads.idはseqkitで取得できる。20スレッド指定。

#1 seqkit seq
seqkit seq -ni outprefix.fastq.gz > reads.id

#2 minimap2 alignment and batching
scripts/create_batched_alignments.sh outprefix.fastq.gz reads.id 20 outdir

outdir/を3で指定する。

 

3、Error-correction

singularityイメージを使う。モデルファイル、バッチサイズ、リードと出力を指定する。推奨バッチサイズは、VRAM40GB(32GBでも可能)のGPUでは64、VRAM80GBのGPUでは128となっている。

singularity run --nv --bind <host_path>:<dest_path> herro.sif inference --read-alns outdir -t <feat_gen_threads_per_device> -d <gpus> -m model_v0.1.pt -b <batch_size> <preprocessed_reads> <fasta_output> 

(レポジトリより)GPU は ID を使って指定する。例えば、パラメータ -d の値が 0,1,3 に設定された場合、herro は 1 番目、2 番目、4 番目の GPU カードを使用する。パラメータ -t はデバイスごとに与えられる - 例えば、-t が 8 に設定され、3 つの GPU が使われる場合、herro は合計で 24 の特徴生成パッドを作成する。

 

 

メモ

  • 論文のR-HERO、F-HERO、L-HEROはRatatosk、FMLRC、LoRDECの3つの反復を指す。
  • 論文では、DBGベースのHybrid error correction (HEC) 手法であるLoRDEC、FMLRC、Ratatosk のいずれかを用いてTGSリードを事前補正している。FMLRC、LoRDEC、およびRatatoskとも、最初の3ラウンドの反復で結果が著しく改善している。
  • レポジトリには、HG002のHifi reads/Duplex ONT reads/Corrected UL readsを使ってHifiasmでアセンブルした結果の表が提示されている。未訂正Duplex ONT readsよりHERO Corrected readsの方が著しくアセンブリが改善されていて、エラー率では及ばないものの、連続性ではHifiリードを上回っている。

引用

Hybrid-hybrid correction of errors in long reads with HERO
Xiongbin Kang, Jialu Xu, Xiao Luo & Alexander Schönhuth 
Genome Biology volume 24, Article number: 275 (2023) 

 

関連

https://kazumaxneo.hatenablog.com/entry/2017/06/20/125551

 

 

複数のプロファイルHMMを1つに統合する HMMerge

 

 過去数十年の間に多重配列アライメントのための手法開発が進歩したにもかかわらず、配列の長さが大きく異なるデータセットのアライメントは、特に入力配列に非常に短い配列(シークエンシング技術、または進化の過程で大きく欠失した配列)が含まれる場合、まだ十分に解決されていない問題である。

HMMergeは、配列長の不均一性が高いデータセットのアラインメントを計算する方法であり、与えられた「バックボーン」アラインメントに短い配列を追加する方法である。HMMergeは、その前身であるアラインメント手法UPPとWITCHの技術を基に構築されており、プロファイルHMMのアンサンブルを構築してバックボーンアラインメントを表現し、そのアンサンブルを用いて残りの配列をバックボーンアラインメントに追加する。HMMergeはUPPやWITCHとは異なり、アンサンブルから新しい「マージ」HMMを構築し、そのマージHMMを使ってクエリー配列をアライメントする。HMMergeはWITCHと競合し、非常に短い配列をバックボーンアラインメントに追加する際にWITCHより優位であることを示す。

 

インストール

condaを使っているので、condaで環境を作って依存を導入した。

依存

  • biopython
  • click
  • dendropy (if using a backbone tree)
  • numpy
  • pyhmmer-sepp
  • pytest (for testing)
  • scipy

Github

https://github.com/MinhyukPark/HMMerge

mamba create -n HMMerge python=3 -y
conda activate HMMerge
pip install biopython click dendropy numpy pyhmmer-sepp pytest scipy pathos

#pyhmmer-seppはpyhmmerのフォーク。オリジナルのpyhmmerなら以下の通り導入可能
mamba install -c bioconda pyhmmer

#例ではtrimAIを使用している
mamba install -c bioconda trimal 

#本体
git clone https://github.com/MinhyukPark/HMMerge.git
cd HMMerge/

python main.py --help

Usage: main.py [OPTIONS]

 

Options:

  --input-dir PATH                The input temp root dir of sepp that

                                  contains all the HMMs  [required]

  --backbone-alignment PATH       The input backbone alignment  [required]

  --query-sequence-file PATH      The input query sequence file  [required]

  --output-prefix PATH            Output prefix  [required]

  --input-type [custom|sepp|upp]  The type of input

  --num-processes INTEGER         Number of Processes

  --support-value FLOAT RANGE     the weigt support of Top HMMs to choose for

                                  merge, 1.0 for all HMMs  [0.0<=x<=1.0]

  --equal-probabilities BOOLEAN   Whether to have equal enty/exit

                                  probabilities

  --model [DNA|RNA|amino]         DNA, RNA, or amino acid analysis  [required]

  --output-format [FASTA|A3M]     FASTA or A3M format for the output alignment

                                  [required]

  --debug                         Whether to run in debug mode or not

  --verbose                       Whether to run in verbose mode or not

  --help                          Show this message and exit.

 


テスト

cd HMMerge/
pytest test.py

 

実行方法

ランするには、FASTA形式のバックボーンアライメント、NEWICK形式のバックボーンツリー、アラインメントするFASTA形式の配列(MSAをdecompose(分解)したもの)が必要。

python main.py --input-dir decomposed_alignments_dir/ --backbone-alignment backbone_alignment --query-sequence-file Query_sequences --output-prefix outdir --num-processes 10 --model dna
  • --input-dir     The input temp root dir of sepp that contains all the HMMs  [required]
  • --backbone-alignment     The input backbone alignment  [required]
  • --query-sequence-file      The input query sequence file  [required]
  • --model [DNA|RNA|amino]      DNA, RNA, or amino acid analysis

 

 

引用

HMMerge: an ensemble method for multiple sequence alignment 
Minhyuk Park, Tandy Warnow
Bioinformatics Advances, Volume 3, Issue 1, 2023

 

kraken2のレポートをkrona plotで視覚化する

2024/02/14 誤字修正

 

メタゲノムデータ解析レシピ(ISBN 978-4-7581-2255-9)3章のWEB年度更新で、kraken2のunclassifiledの割合には注意しましょうという説明をしました。その中で、unclassifiledがkrona plotには反映されないと書いたのですが、これはKrakenの出力をBrackenに読み込み、Brackenで推定したレポートをkrona plotにプロットした時の話でした。Krakenの生のレポートをそのまま使った時にはunclassifiledもkronaでプロットされます。ここで修正させていただきます。

ではなぜBrackenをKrakenと共に使うのでしょうか。今日はこれについて書いてみます。まずkrakenのアルゴリズムについて簡潔に説明します。年度更新でも書いてますが、krakenはk-merの完全一致による超高速なシークエンシングリードの分類学的分類のプロファイラーです(Wood DE et al, 2014)。krakenは各リード配列内の全てのk-merを取り出し、それらを全ゲノムのLCA分類群にマッピングしてk-merの数で重み付けします。ルートからリーフ(種)まで、パスに沿ったすべての重みの合計を計算することによってノードの得点を計算します。最大得点のパスが正しい分類とみなし(通常は一意)、リードそれぞれに対して対応する分類学的ラベルを割り当てます。krakenは塩基またはアミノ酸(kraken2のみ)の完全一致を利用しているため、同種であっても、データベースから少し距離のある株のリードの場合、どの分類群とも一致せず、unclassifiedなリードが多数になってしまう可能性があります。数万以上のリファレンスゲノムにスケールできるものの、データベースが充実していないと正しい分類が難しくなるという欠点があるわけです。2019年にメモリ使用量を数分の一と大幅に削減したkraken2の論文も出ていますが、DB依存性が高い点は変わらずです(Wood DE et al, 2019)。

krakenのアルゴリズムのもう1つの問題点として、種間や属間で高度に保存されたゲノム領域があると、最大得点のパスが種や属レベルで2つ以上出てきてしまい、krakenは種や属のレベルまでユニークにリードを分類することができなくなる点が挙げられます。そのような時、krakenはlast common ancestor (LCA)までリードに分類をアサインします。種レベルでユニークな分類が不可能なら、ルートから属まで、属も無理ならルートから科まで、ということです。しかし、このようなリードが多数になると、種レベルや属レベルの定量をしたい時、誤差が無視できないほど大きくなります。2017年に発表されたBracken(Jennifer Lu et al, 2017)はこの点に注目した実装で、リードを確率的に再分布させることによってデータ中の種の存在量を推定します。再分布させたいランクはユーザーが指定します。例えば種を指定したとしたら、Brackenは、論文の式1にあるベイズ推定によって、種レベルより上のノードにアサインされたリードを種ノードの下に分散し、株レベルでアサインされたリードをそれらの親の種上に再分散して確率的なカウント値を算出しようとします。詳しくはそれぞれの原著論文をご確認下さい。

Brackenがこのような性質を持つため、Brackenをkrakenのあとに実行することで、種や属レベル定量結果を改善できます。

 

1、Kraken2

この流れを実際に確認してみましょう。ここではKraken v2を使います。Kraken2とBrackenは前もってインストールされ、データベースも既にローカル環境にダウンロード・ビルドされているとします(詳しくは年度更新を参照して下さい)。

#way1 single sample
kraken2 --db <path>/<to>/kraken2_DB/ --threads 20 --gzip-compressed --paired --report kraken2_report.txt --use-names --classified-out seqs#.fq R1.fq.gz R2.fq.gz > output.txt

#way2 forループ
mkdir report
for file in `\find *R1.fq.gz -maxdepth 1 -type f`; do
 fastq1=$file
fastq2=${file%_R1.fq.gz}_R2.fq.gz
kraken2 --db <path>/<to>/kraken2_DB/ --threads 20 --gzip-compressed --paired --report report/${file%_R1.fq.gz}_kraken2_report.txt --use-names --classified-out seqs#.fq $fastq1 $fastq2 > report/${file%_R1.fq.gz}_output.txt
rm seqs_1.fq seqs_2.fq
done

#way3 たくさんfastqがあるなら、DB読み込みを1度で実行できるフォークを使用(解説)
mkdir report
kraken2 --use-names --report-zero-counts --threads 20 -db <path>/<to>/kraken2_DB/ --output report/ --report report.txt *fastq

-reportで指定するレポートテキスト:kraken2_report.txtの先頭が下になります。

unclassifiedがまずプリントされて、2行目以降にルートから種までの全リードの分類結果の%がプリントされます。レポートは、検出された分類ごとに繰り返しこの結果をプリントします。

 

2,Bracken

続いてBrackenを使い、指定した分類階級のカウント数を推定します。ここでは種を 指定します(-l S)。krakenのレポートテキスト(-i kraken2_report.txt)とリード長(-r 150)を指定します。"-d"でkraken DBを指定します。

#way1 single sample
bracken -d <path>/<to>/kraken2_DB/ -i kraken2_report.txt -r 150 -l S -t 5 -o  bracken_report.txt

#way2 forループ
for file in `\find *R1.fq.gz -maxdepth 1 -type f`; do
bracken -d <path>/<to>/kraken2_DB/ -i ${file%_R1.fq.gz}_kraken2_report.txt -r 150 -l S -t 5 -o ${file%_QT1.fq.gz}_bracken_report.txt
done

通常のレポートなら数秒以内に結果は出力されます。Brackenのレポートテキスト:bracken_report.txtの先頭が下になります。

種レベルを指定したので、先ほどのKrakenレポートのようにルートから種までの全リードの分類がプリントされるのではなく、1行ずつ種レベルの確率的分類がプリントされているのが特徴です。種レベルなので、アサインされていないunclassifiedも消えています。

また、これとは別にkraken形式のレポートファイルも再作成されて保存されます。こちらも未分類のリードはレポートに含まれません。

kraken形式のレポートファイルですが、Brackenのレポートには推定レベル以下のレベルはプリントされず、カウント値が閾値を下回ったレベルは含まれていません。また、Brackenでは扱えない未分類のリードもレポートに含まれていません(特にここがkrakenレポートと異なる部分)。

 

3、krona

最後にkronaを使ってBrackenの結果を階層的な円グラフ形式で視覚化します。Brackenレポート(画像3枚目のkrakenスタイルのレポートの方)をKraken toolsを使ってkrona互換のテキストファイルに変換し、それからkronaでプロットします。

#step2 brackenからのkraken形式brackenレポートをkrona互換形式に変換
git clone https://github.com/jenniferlu717/KrakenTools.git
python KrakenTools/kreport2krona.py -r kraken2_report.txt -o bracken_convert

#step2
#korna plot。複数レポートがあるならワイルドカード指定
ktImportText *bracken_convert -o krona.html

krona plot

 

このようにプロットできますが、年度更新で説明した通り、unclassifiedはKrakenの精度を測るための重要な指標になります。参考資料(*1)のようにレポート先頭を別に保存しておくか、Brackenを介さずにkrona plotを描画しておきましょう。つまり、説明の1のステップ、3のステップの順番です。

#step1のkrakenからのレポートをkrona互換形式に変換
git clone https://github.com/jenniferlu717/KrakenTools.git
python KrakenTools/kreport2krona.py -r kraken2_report.txt -o kraken2_convert

#step2
#korna plot。複数レポートがあるならワイルドカード指定
ktImportText *kraken2_convert -o krona.html

krona plot

Krakenのレポートテキストを使ったため、今度のプロットにはunclassifiedもプロットされています。37%がunclassifiedであることが分かりました(リアルメタゲノム使用)。

 

このように、Krakenのレポートもkronaで可視化しておくとunclassifiedも視覚化されます。unclassifiedはKrakenの信頼性を見積もるために重要な指標になるので、KrakenのレポートもBrackenからのレポートもどちらもプロットしておくと良いのではないかと思います。

引用

Bracken: estimating species abundance in metagenomics data

Jennifer Lu​, Florian P. Breitwieser, Peter Thielen, Steven L. Salzberg

January 2, 2017 PeerJ

https://peerj.com/articles/cs-104/

 

Improved metagenomic analysis with Kraken 2
Derrick E. Wood, Jennifer Lu & Ben Langmead 
Genome Biology volume 20, Article number: 257 (2019) 

 

参考

*1

https://telatin.github.io/microbiome-bioinformatics/Metagenomics-classification-2/

 

注;ここで紹介したフローは一例です。krakenの出力をどのように修正や整形、マージ、視覚化していくかは複数のやり方があります。

 

 

 

 

 

ノイズの多いロングリードからハプロタイプを考慮したde novo二倍体ゲノムアセンブリを行う PECAT

 

 高いシーケンスエラーは、2倍体ゲノムアセンブリへのロングノイズリードの適用を妨げてきた。既存のアセンブラーでは、長ノイズリードに含まれる高シーケンスエラーとヘテロ接合体を区別できず、ハプロタイプスイッチの多いアセンブリーが生成されてしまう。ここでは、長いノイズリードから2倍体ゲノムを再構築するための、段階的エラー修正およびアセンブルツールであるPECATを紹介する。ハプロタイプを考慮したエラー修正法を設計し、シーケンスエラーを修正しながらヘテロ接合体の対立遺伝子を保持することができる。さらに、補正されたリードのSNPエラーを減少させるリードレベルSNPコーラーを開発した。次に、リードを異なるハプロタイプグループに割り当てるために、リードグルーピング法を用いる。アセンブルを高速化するため、PECATは必要な場合にのみローカルアライメントを行う。PECATは、他のアセンブラと比較して、ノイズの多い長いリードのみを用いて二倍体ゲノムを効率的にアセンブルし、より連続したハプロタイプ特異的コンティグを生成する。特に、B. taurus (Bison×Simmental)においては、Nanoporeリードを用いることで、ほぼハプロタイプを分離したアセンブルを実現している。

 

インストール

依存

  • python3 (3.6+)
  • minimap2 (2.17+)
  • racon (v1.4.21+)
  • perl (v5.22.1+)
  • samtools (1.7+)
  • clair3 (v0.1-r12+) (optional)
  • medaka (1.7.2+) (optional)

Github

#依存関係
mamba create -n pecat-env python=3.11 -y
conda activate pecat-env
mamba install minimap2=2.24 racon=1.5 perl=5.32 samtools=1.17 -y

#本体
git clone --recursive https://github.com/lemene/PECAT.git
cd PECAT
make -j 12

#pre-build docker image(dockerhub)
docker pull lemene/pecat:v0.0.3

$ ./pecat.pl 

Usage: pecat.pl correct|assemble|unzip|config cfg_fname

    correct:     correct rawreads

    assemble:    generate haplotype-collapsed assembly

    unzip:       generate diploid assembly

    config:      generate default config file

 

 

テストラン

PECAT/demo/.

cfgfile

 

 

cd PECAT/demo/
#リードを指定してデフォルトパラメータのコンフィグファイル作成
pecat.pl config cfgfile
=> cfgfileができる。

#アセンブリ。configファイルを指定する。
../build/demo/pecat.pl unzip cfgfile

出力

S1/

S1/6_polish/racon/

  • 修正されたリードはS1/1-correct/corrected_reads.fastaファイルにある。
  • primaryコンティグとalternateコンティグはS1/6-polish/racon/にある。
  • dual-formatコンティグはS1/6-polish/racon/にある(haplotype_1.fasta, haplotype_2.fasta)。
  • polish_medaka=1が設定されている場合、PECATはMedakaを使用して上記の結果をさらにポリッシュし、コンティグはS1/6-polish/medakaに置かれる。

 

その他(レポジトリより)

  • cleanup=1を推奨。PECATはテンポラリファイルを削除する。
  • ウシやヒトのような大きなゲノムの場合、corr_rd2rd_optionsとalign_rd2rd_optionsに -f 0.005 または -f 0.002 というパラメータを追加することを推奨する。demo/configs/にあるcfg_cattle_clr, cfg_cattle_ont, cfg_hg002_ontを参照。このパラメータはminimap2に渡され、繰り返されるminimizersの上位0.005または0.002フラクションをフィルタリングすることを意味する。これにより、重複候補が少なくなり、ディスク使用量が減り、エラー修正ステップとアセンブルステップが速くなる。

引用

de novo diploid genome assembly using long noisy reads via haplotype-aware error correction and inconsistent overlap identification
Fan Nie, Neng Huang, Jun Zhang, Peng Ni, Zhenyu Wang, ChuanLe Xiao, Feng Luo, Jianxin Wang

bioRxiv, Posted September 27, 2022. 

bioRxiv, Posted February 17, 2023 (new version).

 

関連

 

このツールは内藤先生のツイートで知りました。ありがとうございます。