macでインフォマティクス

macでインフォマティクス

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

大規模微生物データセットのためのスケーラブルなコアゲノムアライメント Parsnp 2.0

高速なツール 2024 Preprint phylogenetic tree Inference

 

　2016年以降、NCBIでリファレンスゲノムが利用可能な微生物種の数は3倍以上に増えている。Multiple genome alignmentは、共通の祖先を共有する複数のゲノムのヌクレオチドを特定するプロセスであり、多くの下流の比較解析手法の入力として使用される。Parsnpは、現在のゲノムデータ時代に対応できる数少ないマルチプルゲノムアライメント手法の1つであるが、2014年の最初のリリース以来、メジャーリリースは行われていない。このギャップを解決するために、本著者らはParsnp v2を開発した。Parsnp v2では、ユーザーがプログラムの実行をよりコントロールできるようになり、Parsnpをさまざまなユースケースに合わせて調整できるようになった。Parsnpにパーティショニングオプションを導入し、入力を複数の並列アライメントプロセスに分割し、最終的なアライメントにまとめることができるようになった。パーティショニングオプションにより、正確なコアゲノムアライメントを維持しながら、メモリ使用量を4倍以上、実行時間を2倍以上削減することができる。また、アライメントのアンカーは、入力セット全体ではなく、そのパーティション内で保存されていればよいため、パーティショニングワークフローは、アセンブリアーティファクトや小さなばらつきによって引き起こされる複雑さの影響を受けにくい。数千の細菌ゲノムとウイルスゲノムを含むデータセットでのパフォーマンスに注目する。

 

チュートリアル

https://harvest.readthedocs.io/en/latest/content/parsnp/tutorial.html

 

インストール

Github

https://github.com/marbl/parsnp?tab=readme-ov-file

#conda(link)
mamba install -c bioconda parsnp=2.0.3 -y

> parsnp -h

01:07:32 - INFO - |--Parsnp 2.0.3--|

 

usage: parsnp [-h] [-c] -d SEQUENCES [SEQUENCES ...] [-r REFERENCE] [-g GENBANK [GENBANK ...]] [-o OUTPUT_DIR] [-q QUERY] [-U MAX_MUMI_DISTR_DIST | -mmd MAX_MUMI_DISTANCE] [-F] [-M] [--use-ani] [--min-ani MIN_ANI] [--min-ref-cov MIN_REF_COV] [--use-mash]

              [--max-mash-dist MAX_MASH_DIST] [-a MIN_ANCHOR_LENGTH] [-m MUM_LENGTH] [-C MAX_CLUSTER_D] [-z MIN_CLUSTER_SIZE] [-D MAX_DIAG_DIFF] [-n {mafft,muscle,fsa,prank}] [-u] [--partition] [--partition-size PARTITION_SIZE] [--extend-lcbs]

              [--extend-ani-cutoff EXTEND_ANI_CUTOFF] [--extend-indel-cutoff EXTEND_INDEL_CUTOFF] [--match-score MATCH_SCORE] [--mismatch-penalty MISMATCH_PENALTY] [--gap-penalty GAP_PENALTY] [--skip-phylogeny] [--validate-input] [--use-fasttree] [--vcf]

              [-p THREADS] [-P MAX_PARTITION_SIZE] [-v] [-x] [-i INIFILE] [-e] [--no-recruit] [-V]

 

    Parsnp quick start for three example scenarios:

    1) With reference & genbank file:

    python Parsnp.py -g <reference_genbank_file1 reference_genbank_file2 ...> -d <seq_file1 seq_file2 ...>  -p <threads>

 

    2) With reference but without genbank file:

    python Parsnp.py -r <reference_genome> -d <seq_file1 seq_file2 ...> -p <threads>

 

    3) Autorecruit reference to a draft assembly:

    python Parsnp.py -q <draft_assembly> -d <seq_file1 seq_file2 ...> -p <threads>

    

 

optional arguments:

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

 

Input/Output:

  -c, --curated         (c)urated genome directory, use all genomes in dir and ignore MUMi?

  -d SEQUENCES [SEQUENCES ...], --sequences SEQUENCES [SEQUENCES ...]

                        A list of files containing genomes/contigs/scaffolds. If the file ends in .txt, each line in the file corresponds to the path to an input file.

  -r REFERENCE, --reference REFERENCE

                        (r)eference genome (set to ! to pick random one from sequence dir)

  -g GENBANK [GENBANK ...], --genbank GENBANK [GENBANK ...]

                        A list of Genbank file(s) (gbk)

  -o OUTPUT_DIR, --output-dir OUTPUT_DIR

  -q QUERY, --query QUERY

                        Specify (assembled) query genome to use, in addition to genomes found in genome dir

 

MUMi:

  -U MAX_MUMI_DISTR_DIST, --max-mumi-distr-dist MAX_MUMI_DISTR_DIST, --MUMi MAX_MUMI_DISTR_DIST

                        Max MUMi distance value for MUMi distribution

  -mmd MAX_MUMI_DISTANCE, --max-mumi-distance MAX_MUMI_DISTANCE

                        Max MUMi distance (default: autocutoff based on distribution of MUMi values)

  -F, --fastmum         Fast MUMi calculation

  -M, --mumi_only, --onlymumi

                        Calculate MUMi and exit? overrides all other choices!

  --use-ani             Use ani for genome recruitment

  --min-ani MIN_ANI     Min ANI value to allow for genome recruitment.

  --min-ref-cov MIN_REF_COV

                        Minimum percent of reference segments to be covered in FastANI

  --use-mash            Use mash for genome recruitment

  --max-mash-dist MAX_MASH_DIST

                        Max mash distance.

 

MUM search:

  -a MIN_ANCHOR_LENGTH, --min-anchor-length MIN_ANCHOR_LENGTH, --anchorlength MIN_ANCHOR_LENGTH

                        Min (a)NCHOR length (default = 1.1*(Log(S)))

  -m MUM_LENGTH, --mum-length MUM_LENGTH, --mumlength MUM_LENGTH

                        Mum length

  -C MAX_CLUSTER_D, --max-cluster-d MAX_CLUSTER_D, --clusterD MAX_CLUSTER_D

                        Maximal cluster D value

  -z MIN_CLUSTER_SIZE, --min-cluster-size MIN_CLUSTER_SIZE, --minclustersize MIN_CLUSTER_SIZE

                        Minimum cluster size

 

LCB alignment:

  -D MAX_DIAG_DIFF, --max-diagonal-difference MAX_DIAG_DIFF, --DiagonalDiff MAX_DIAG_DIFF

                        Maximal diagonal difference. Either percentage (e.g. 0.2) or bp (e.g. 100bp)

  -n {mafft,muscle,fsa,prank}, --alignment-program {mafft,muscle,fsa,prank}, --alignmentprog {mafft,muscle,fsa,prank}

                        Alignment program to use

  -u, --unaligned       Output unaligned regions

 

Sequence Partitioning:

  --partition           Evenly split input sequences across separate runs of parsnp, then merge results

  --partition-size PARTITION_SIZE

                        Number of sequences in a partitioned run of parsnp

 

LCB Extension:

  --extend-lcbs         Extend the boundaries of LCBs with an ungapped alignment

  --extend-ani-cutoff EXTEND_ANI_CUTOFF

                        Cutoff ANI for lcb extension

  --extend-indel-cutoff EXTEND_INDEL_CUTOFF

                        Cutoff for indels in LCB extension region. LCB extension will be at most min(seqs) + cutoff bases

  --match-score MATCH_SCORE

                        Value of match score for extension

  --mismatch-penalty MISMATCH_PENALTY

                        Value of mismatch score for extension (should be negative)

  --gap-penalty GAP_PENALTY

                        Value of gap penalty for extension (should be negative)

 

Misc:

  --skip-phylogeny      Do not generate phylogeny from core SNPs

  --validate-input      Use Biopython to validate input files

  --use-fasttree        Use fasttree instead of RaxML

  --vcf                 Generate VCF file.

  -p THREADS, --threads THREADS

                        Number of threads to use

  -P MAX_PARTITION_SIZE, --max-partition-size MAX_PARTITION_SIZE

                        Max partition size (limits memory usage)

  -v, --verbose         Verbose output

  -x, --xtrafast

  -i INIFILE, --inifile INIFILE, --ini-file INIFILE

  -e, --extend

  --no-recruit

  -V, --version         show program's version number and exit

 

 

 

実行方法

-gでリファレンスゲノムのgenbankファイル、もしくは-rでリファレンスゲノムのfastaファイルを指定する。-dで比較したいバクテリアのfastaゲノムを含むディレクトリを指定する。

parsnp -g ref.gbk -d genome_dir
  • -g    GENBANK
  • -d   A list of files containing genomes/contigs/scaffolds. If the file ends in .txt, each line in the file corresponds to the path to an input file.
  • -r    reference genome (set to ! to pick random one from sequence dir)

 

（レポジトリより）Parsnp2には、-partitionで起動する新しいモードがある。このモードでは、入力ゲノムをランダムにp個ずつのゲノムグループに分割する。pはデフォルトで50で、-partition-size=pで変更できる。Parsnpは各グループに対して独立に実行され、各グループのアライメント結果は全入力ゲノムの1つのアライメントにマージされる。このモードは、計算量を減らすことができ、大規模なデータセットを対象としている。

parsnp -r examples/mers_virus/ref/England1.fna -d examples/mers_virus/genomes/*.fna --partition --partition-size 10 -o examples-out-partitioned

 

出力

  • parsnp.xmfa - コアゲノムアライメント。
  • parsnp.ggr - harvest-toolkit によって生成されたアラインメントを圧縮したもの。このファイルは Gingr でアラインメントを可視化するために使用できる。
  • parsnp.snps.mblocks - 各配列のコア SNP シグネチャを fasta 形式で表したファイル。parsnp.treeを生成するために使用される。
  • parsnp.tree - 結果の系統樹。
  • パーティションモードで実行した場合、Parsnpは出力ディレクトリにパーティションフォルダを作成する。

 

引用

Parsnp 2.0: Scalable Core-Genome Alignment for Massive Microbial Datasets
Bryce Kille, Michael G Nute, Victor Huang, Eddie Kim, Adam M. Phillippy,  Todd J. Treangen

bioRxiv, Posted January 31, 2024.

 

関連

https://kazumaxneo.hatenablog.com/entry/2017/11/26/205234