macでインフォマティクス (original) (raw)

メタゲノミクスは、環境やヒトに関連するマイクロバイオーム研究に革命をもたらした。しかし、生物学的プロセスや分子機能が知られているタンパク質の数は限られており、これが大きなボトルネックとなっている。原核生物やウイルスでは、同じ生物学的プロセスに関与する遺伝子を、保存された遺伝子クラスターとして共局在させておくことが進化上好ましい。逆に、遺伝子近傍の保存は機能的関連を示す。Spacedustは、保存された遺伝子クラスターを系統的にde novoで発見するためのツールである。相同タンパク質のマッチを見つけるために、Foldseekとの高速で高感度な構造比較を使用する。部分的に保存されたクラスターは、新しいクラスタリングと次数保存のP値を用いて検出される。1,308の細菌ゲノムの全対全分析でSpacedustの感度を実証し、420万遺伝子の58%を含む72,843の保存遺伝子クラスターを同定した。また、特殊なツールによってアノテーションされた抗ウイルス防御系クラスターの95%を回復した。Spacedustの高い感度とスピードは、配列決定された膨大な数の細菌、古細菌、ウイルスゲノムの大規模アノテーションを促進するだろう。

インストール

ubuntuu22.04LTSで公開されているstatic binaryを使用した（CPU: TR3990X）。

Github

static Linux AVX2 build

wget https://mmseqs.com/spacedust/spacedust-linux-avx2.tar.gz; tar xvzf spacedust-linux-avx2.tar.gz; export PATH=$(pwd)/spacedust/bin/:$PATH

static Linux SSE4.1 build

wget https://mmseqs.com/spacedust/spacedust-linux-sse41.tar.gz; tar xvzf spacedust-linux-sse41.tar.gz; export PATH=$(pwd)/spacedust/bin/:$PATH

static macOS build (universal binary with SSE4.1/AVX2/M1 NEON)

wget https://mmseqs.com/spacedust/spacedust-osx-universal.tar.gz; tar xvzf spacedust-osx-universal.tar.gz; export PATH=$(pwd)/spacedust/bin/:$PATH

> spacedust -h

$ spacedust -h

Spacedust is a tool to discover conserved gene clusters between any pairs of contig/genomes

spacedust Version: 16b020301be952232d6eb2eaa2cd2ad0933d68b0

© Ruoshi Zhang ruoshi.zhang@mpinat.mpg.de

usage: spacedust []

Main workflows for database input/output

createsetdb Create sequence set database from FASTA (and GFF3) input of contigs/genomes

aa2foldseek Map a sequence DB to reference foldseek DB

clusterdb Build a searchable cluster database from sequence DB or foldseek structure DB

clustersearch Find clusters of colocalized hits between any query-target sequence/profile set database

Special-purpose utilities

besthitbyset For each set of sequences compute the best element and update p-value

combinehits Group hits and compute a combined E-value for each query-target set pair

summarizeresults Summarize results on clustered hits

clusterhits Find clusters of hits by agglomerative hierarchical clustering and compute their clustering and ordering P-values

>spacedust createsetdb -h

$ spacedust createsetdb -h

usage: spacedust createsetdb <i:fastaFile1[.gz|bz2]> ... <i:fastaFileN[.gz|bz2]> <o:setDB> [options]

By Ruoshi Zhang ruoshi.zhang@mpinat.mpg.de & Milot Mirdita milot@mirdita.de

options: misc:

--dbtype INT Database type 0: auto, 1: amino acid 2: nucleotides [0]

--shuffle BOOL Shuffle input database [0]

--createdb-mode INT Createdb mode 0: copy data, 1: soft link data and write new index (works only with single line fasta/q) [0]

--id-offset INT Numeric ids in index file are offset by this value [0]

--min-length INT Minimum codon number in open reading frames [30]

--max-length INT Maximum codon number in open reading frames [32734]

--max-gaps INT Maximum number of codons with gaps or unknown residues before an open reading frame is rejected [2147483647]

--contig-start-mode INT Contig start can be 0: incomplete, 1: complete, 2: both [2]

--contig-end-mode INT Contig end can be 0: incomplete, 1: complete, 2: both [2]

--orf-start-mode INT Orf fragment can be 0: from start to stop, 1: from any to stop, 2: from last encountered start to stop (no start in the middle) [1]

--forward-frames STR Comma-separated list of frames on the forward strand to be extracted [1,2,3]

--reverse-frames STR Comma-separated list of frames on the reverse strand to be extracted [1,2,3]

--translation-table INT 1) CANONICAL, 2) VERT_MITOCHONDRIAL, 3) YEAST_MITOCHONDRIAL, 4) MOLD_MITOCHONDRIAL, 5) INVERT_MITOCHONDRIAL, 6) CILIATE

  1. FLATWORM_MITOCHONDRIAL, 10) EUPLOTID, 11) PROKARYOTE, 12) ALT_YEAST, 13) ASCIDIAN_MITOCHONDRIAL, 14) ALT_FLATWORM_MITOCHONDRIAL

  2. BLEPHARISMA, 16) CHLOROPHYCEAN_MITOCHONDRIAL, 21) TREMATODE_MITOCHONDRIAL, 22) SCENEDESMUS_MITOCHONDRIAL

  3. THRAUSTOCHYTRIUM_MITOCHONDRIAL, 24) PTEROBRANCHIA_MITOCHONDRIAL, 25) GRACILIBACTERIA, 26) PACHYSOLEN, 27) KARYORELICT, 28) CONDYLOSTOMA

  4. MESODINIUM, 30) PERTRICH, 31) BLASTOCRITHIDIA [1]

--translate INT Translate ORF to amino acid [0]

--use-all-table-starts BOOL Use all alternatives for a start codon in the genetic table, if false - only ATG (AUG) [0]

--add-orf-stop BOOL Add stop codon '*' at complete start and end [0]

--gff-type STR Comma separated list of feature types in the GFF file to select

--stat STR One of: linecount, mean, min, max, doolittle, charges, seqlen, firstline

--tsv BOOL Return output in TSV format [0]

--gff-dir STR Path to gff dir file

common:

--compressed INT Write compressed output [0]

-v INT Verbosity level: 0: quiet, 1: +errors, 2: +warnings, 3: +info [3]

--threads INT Number of CPU-cores used (all by default) [128]

expert:

--write-lookup INT write .lookup file containing mapping from internal id, fasta id and file number [1]

--create-lookup INT Create database lookup file (can be very large) [0]

references:

- Steinegger M, Soding J: MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology, 35(11), 1026-1028 (2017)

> spacedust aa2foldseek -h

$ spacedust aa2foldseek -h

usage: spacedust aa2foldseek <i:inputDB> <i:targetDB> [options]

By Ruoshi Zhang ruoshi.zhang@mpinat.mpg.de & Milot Mirdita milot@mirdita.de

options: prefilter:

--seed-sub-mat TWIN Substitution matrix file for k-mer generation [aa:VTML80.out,nucl:nucleotide.out]

-s FLOAT Sensitivity: 1.0 faster; 4.0 fast; 7.5 sensitive [4.000]

-k INT k-mer length (0: automatically set to optimum) [0]

--k-score TWIN k-mer threshold for generating similar k-mer lists [seq:2147483647,prof:2147483647]

--alph-size TWIN Alphabet size (range 2-21) [aa:21,nucl:5]

--max-seqs INT Maximum results per query sequence allowed to pass the prefilter (affects sensitivity) [10]

--split INT Split input into N equally distributed chunks. 0: set the best split automatically [0]

--split-mode INT 0: split target db; 1: split query db; 2: auto, depending on main memory [2]

--split-memory-limit BYTE Set max memory per split. E.g. 800B, 5K, 10M, 1G. Default (0) to all available system memory [0]

--comp-bias-corr INT Correct for locally biased amino acid composition (range 0-1) [1]

--comp-bias-corr-scale FLOAT Correct for locally biased amino acid composition (range 0-1) [1.000]

--diag-score BOOL Use ungapped diagonal scoring during prefilter [1]

--exact-kmer-matching INT Extract only exact k-mers for matching (range 0-1) [1]

--mask INT Mask sequences in k-mer stage: 0: w/o low complexity masking, 1: with low complexity masking [1]

--mask-prob FLOAT Mask sequences is probablity is above threshold [0.900]

--mask-lower-case INT Lowercase letters will be excluded from k-mer search 0: include region, 1: exclude region [0]

--min-ungapped-score INT Accept only matches with ungapped alignment score above threshold [15]

--add-self-matches BOOL Artificially add entries of queries with themselves (for clustering) [0]

--spaced-kmer-mode INT 0: use consecutive positions in k-mers; 1: use spaced k-mers [1]

--spaced-kmer-pattern STR User-specified spaced k-mer pattern

--local-tmp STR Path where some of the temporary files will be created

align:

-c FLOAT List matches above this fraction of aligned (covered) residues (see --cov-mode) [0.900]

--cov-mode INT 0: coverage of query and target

1: coverage of target

2: coverage of query

3: target seq. length has to be at least x% of query length

4: query seq. length has to be at least x% of target length

5: short seq. needs to be at least x% of the other seq. length [0]

-a BOOL Add backtrace string (convert to alignments with mmseqs convertalis module) [0]

--alignment-mode INT How to compute the alignment:

0: automatic

1: only score and end_pos

2: also start_pos and cov

3: also seq.id

4: only ungapped alignment [0]

--alignment-output-mode INT How to compute the alignment:

0: automatic

1: only score and end_pos

2: also start_pos and cov

3: also seq.id

4: only ungapped alignment

5: score only (output) cluster format [0]

--wrapped-scoring BOOL Double the (nucleotide) query sequence during the scoring process to allow wrapped diagonal scoring around end and start [0]

-e DOUBLE List matches below this E-value (range 0.0-inf) [1.000E-03]

--min-seq-id FLOAT List matches above this sequence identity (for clustering) (range 0.0-1.0) [0.900]

--min-aln-len INT Minimum alignment length (range 0-INT_MAX) [0]

--seq-id-mode INT 0: alignment length 1: shorter, 2: longer sequence [0]

--alt-ali INT Show up to this many alternative alignments [0]

--max-rejected INT Maximum rejected alignments before alignment calculation for a query is stopped [2147483647]

--max-accept INT Maximum accepted alignments before alignment calculation for a query is stopped [2147483647]

--score-bias FLOAT Score bias when computing SW alignment (in bits) [0.000]

--realign BOOL Compute more conservative, shorter alignments (scores and E-values not changed) [0]

--realign-score-bias FLOAT Additional bias when computing realignment [-0.200]

--realign-max-seqs INT Maximum number of results to return in realignment [2147483647]

--corr-score-weight FLOAT Weight of backtrace correlation score that is added to the alignment score [0.000]

--gap-open TWIN Gap open cost [aa:11,nucl:5]

--gap-extend TWIN Gap extension cost [aa:1,nucl:2]

--zdrop INT Maximal allowed difference between score values before alignment is truncated (nucleotide alignment only) [40]

profile:

--pca Pseudo count admixture strength

--pcb Pseudo counts: Neff at half of maximum admixture (range 0.0-inf)

misc:

--taxon-list STR Taxonomy ID, possibly multiple values separated by ','

--stat STR One of: linecount, mean, min, max, doolittle, charges, seqlen, firstline

--tsv BOOL Return output in TSV format [0]

common:

--sub-mat TWIN Substitution matrix file [aa:blosum62.out,nucl:nucleotide.out]

--max-seq-len INT Maximum sequence length [65535]

--db-load-mode INT Database preload mode 0: auto, 1: fread, 2: mmap, 3: mmap+touch [0]

--threads INT Number of CPU-cores used (all by default) [128]

--compressed INT Write compressed output [0]

-v INT Verbosity level: 0: quiet, 1: +errors, 2: +warnings, 3: +info [3]

--remove-tmp-files BOOL Delete temporary files [0]

references:

- Steinegger M, Soding J: MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology, 35(11), 1026-1028 (2017)

> spacedust clusterdb -h

$ spacedust clusterdb -h

usage: spacedust clusterdb <i:inputDB> [options]

By Ruoshi Zhang ruoshi.zhang@mpinat.mpg.de & Milot Mirdita milot@mirdita.de

options: prefilter:

--seed-sub-mat TWIN Substitution matrix file for k-mer generation [aa:VTML80.out,nucl:nucleotide.out]

-s FLOAT Sensitivity: 1.0 faster; 4.0 fast; 7.5 sensitive [4.000]

-k INT k-mer length (0: automatically set to optimum) [0]

--k-score TWIN k-mer threshold for generating similar k-mer lists [seq:2147483647,prof:2147483647]

--alph-size TWIN Alphabet size (range 2-21) [aa:21,nucl:5]

--max-seqs INT Maximum results per query sequence allowed to pass the prefilter (affects sensitivity) [300]

--split INT Split input into N equally distributed chunks. 0: set the best split automatically [0]

--split-mode INT 0: split target db; 1: split query db; 2: auto, depending on main memory [2]

--split-memory-limit BYTE Set max memory per split. E.g. 800B, 5K, 10M, 1G. Default (0) to all available system memory [0]

--comp-bias-corr INT Correct for locally biased amino acid composition (range 0-1) [1]

--comp-bias-corr-scale FLOAT Correct for locally biased amino acid composition (range 0-1) [1.000]

--diag-score BOOL Use ungapped diagonal scoring during prefilter [1]

--exact-kmer-matching INT Extract only exact k-mers for matching (range 0-1) [0]

--mask INT Mask sequences in k-mer stage: 0: w/o low complexity masking, 1: with low complexity masking [1]

--mask-prob FLOAT Mask sequences is probablity is above threshold [0.900]

--mask-lower-case INT Lowercase letters will be excluded from k-mer search 0: include region, 1: exclude region [0]

--min-ungapped-score INT Accept only matches with ungapped alignment score above threshold [15]

--add-self-matches BOOL Artificially add entries of queries with themselves (for clustering) [0]

--spaced-kmer-mode INT 0: use consecutive positions in k-mers; 1: use spaced k-mers [1]

--spaced-kmer-pattern STR User-specified spaced k-mer pattern

--local-tmp STR Path where some of the temporary files will be created

align:

-c FLOAT List matches above this fraction of aligned (covered) residues (see --cov-mode) [0.800]

--cov-mode INT 0: coverage of query and target

1: coverage of target

2: coverage of query

3: target seq. length has to be at least x% of query length

4: query seq. length has to be at least x% of target length

5: short seq. needs to be at least x% of the other seq. length [0]

-a BOOL Add backtrace string (convert to alignments with mmseqs convertalis module) [0]

--alignment-mode INT How to compute the alignment:

0: automatic

1: only score and end_pos

2: also start_pos and cov

3: also seq.id

4: only ungapped alignment [0]

--alignment-output-mode INT How to compute the alignment:

0: automatic

1: only score and end_pos

2: also start_pos and cov

3: also seq.id

4: only ungapped alignment

5: score only (output) cluster format [0]

--wrapped-scoring BOOL Double the (nucleotide) query sequence during the scoring process to allow wrapped diagonal scoring around end and start [0]

-e DOUBLE List matches below this E-value (range 0.0-inf) [1.000E-03]

--min-seq-id FLOAT List matches above this sequence identity (for clustering) (range 0.0-1.0) [0.700]

--min-aln-len INT Minimum alignment length (range 0-INT_MAX) [0]

--seq-id-mode INT 0: alignment length 1: shorter, 2: longer sequence [0]

--alt-ali INT Show up to this many alternative alignments [0]

--max-rejected INT Maximum rejected alignments before alignment calculation for a query is stopped [2147483647]

--max-accept INT Maximum accepted alignments before alignment calculation for a query is stopped [2147483647]

--score-bias FLOAT Score bias when computing SW alignment (in bits) [0.000]

--realign BOOL Compute more conservative, shorter alignments (scores and E-values not changed) [0]

--realign-score-bias FLOAT Additional bias when computing realignment [-0.200]

--realign-max-seqs INT Maximum number of results to return in realignment [2147483647]

--corr-score-weight FLOAT Weight of backtrace correlation score that is added to the alignment score [0.000]

--gap-open TWIN Gap open cost [aa:11,nucl:5]

--gap-extend TWIN Gap extension cost [aa:1,nucl:2]

--zdrop INT Maximal allowed difference between score values before alignment is truncated (nucleotide alignment only) [40]

clust:

--cluster-mode INT 0: Set-Cover (greedy)

1: Connected component (BLASTclust)

2,3: Greedy clustering by sequence length (CDHIT) [0]

--max-iterations INT Maximum depth of breadth first search in connected component clustering [1000]

--similarity-type INT Type of score used for clustering. 1: alignment score 2: sequence identity [2]

--single-step-clustering BOOL Switch from cascaded to simple clustering workflow [0]

--cluster-steps INT Cascaded clustering steps from 1 to -s [3]

--cluster-reassign BOOL Cascaded clustering can cluster sequence that do not fulfill the clustering criteria.

Cluster reassignment corrects these errors [0]

kmermatcher:

--weights STR Weights used for cluster priorization

--cluster-weight-threshold FLOAT Weight threshold used for cluster priorization [0.900]

--kmer-per-seq INT k-mers per sequence [21]

--kmer-per-seq-scale TWIN Scale k-mer per sequence based on sequence length as kmer-per-seq val + scale x seqlen [aa:0.000,nucl:0.200]

--adjust-kmer-len BOOL Adjust k-mer length based on specificity (only for nucleotides) [0]

--hash-shift INT Shift k-mer hash initialization [67]

--include-only-extendable BOOL Include only extendable [0]

--ignore-multi-kmer BOOL Skip k-mers occurring multiple times (>=2) [0]

profile:

--pca Pseudo count admixture strength

--pcb Pseudo counts: Neff at half of maximum admixture (range 0.0-inf)

misc:

--taxon-list STR Taxonomy ID, possibly multiple values separated by ','

--rescore-mode INT Rescore diagonals with:

0: Hamming distance

1: local alignment (score only)

2: local alignment

3: global alignment

4: longest alignment fulfilling window quality criterion [0]

--search-mode INT 0: sequence search with MMseqs2, 1: structure comparison with Foldseek [0]

common:

--sub-mat TWIN Substitution matrix file [aa:blosum62.out,nucl:nucleotide.out]

--max-seq-len INT Maximum sequence length [65535]

--db-load-mode INT Database preload mode 0: auto, 1: fread, 2: mmap, 3: mmap+touch [0]

--threads INT Number of CPU-cores used (all by default) [128]

--compressed INT Write compressed output [0]

-v INT Verbosity level: 0: quiet, 1: +errors, 2: +warnings, 3: +info [3]

--remove-tmp-files BOOL Delete temporary files [0]

--force-reuse BOOL Reuse tmp filse in tmp/latest folder ignoring parameters and version changes [0]

--mpi-runner STR Use MPI on compute cluster with this MPI command (e.g. "mpirun -np 42")

expert:

--filter-hits BOOL Filter hits by seq.id. and coverage [0]

--sort-results INT Sort results: 0: no sorting, 1: sort by E-value (Alignment) or seq.id. (Hamming) [0]

references:

- Steinegger M, Soding J: MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology, 35(11), 1026-1028 (2017)

> spacedust clustersearch -h

$ spacedust clustersearch -h

usage: spacedust clustersearch <i:querySetDB> <i:targetSetDB> <o:output[.tsv]> [options]

By Ruoshi Zhang ruoshi.zhang@mpinat.mpg.de & Milot Mirdita milot@mirdita.de

options: prefilter:

--comp-bias-corr INT Correct for locally biased amino acid composition (range 0-1) [1]

--comp-bias-corr-scale FLOAT Correct for locally biased amino acid composition (range 0-1) [1.000]

--add-self-matches BOOL Artificially add entries of queries with themselves (for clustering) [0]

--seed-sub-mat TWIN Substitution matrix file for k-mer generation [aa:VTML80.out,nucl:nucleotide.out]

-s FLOAT Sensitivity: 1.0 faster; 4.0 fast; 7.5 sensitive [5.700]

-k INT k-mer length (0: automatically set to optimum) [0]

--k-score TWIN k-mer threshold for generating similar k-mer lists [seq:2147483647,prof:2147483647]

--alph-size TWIN Alphabet size (range 2-21) [aa:21,nucl:5]

--max-seqs INT Maximum results per query sequence allowed to pass the prefilter (affects sensitivity) [300]

--split INT Split input into N equally distributed chunks. 0: set the best split automatically [0]

--split-mode INT 0: split target db; 1: split query db; 2: auto, depending on main memory [2]

--split-memory-limit BYTE Set max memory per split. E.g. 800B, 5K, 10M, 1G. Default (0) to all available system memory [0]

--diag-score BOOL Use ungapped diagonal scoring during prefilter [1]

--exact-kmer-matching INT Extract only exact k-mers for matching (range 0-1) [0]

--mask INT Mask sequences in k-mer stage: 0: w/o low complexity masking, 1: with low complexity masking [1]

--mask-prob FLOAT Mask sequences is probablity is above threshold [0.900]

--mask-lower-case INT Lowercase letters will be excluded from k-mer search 0: include region, 1: exclude region [0]

--min-ungapped-score INT Accept only matches with ungapped alignment score above threshold [15]

--spaced-kmer-mode INT 0: use consecutive positions in k-mers; 1: use spaced k-mers [1]

--spaced-kmer-pattern STR User-specified spaced k-mer pattern

--local-tmp STR Path where some of the temporary files will be created

--disk-space-limit BYTE Set max disk space to use for reverse profile searches. E.g. 800B, 5K, 10M, 1G. Default (0) to all available disk space in the temp folder [0]

align:

-a BOOL Add backtrace string (convert to alignments with mmseqs convertalis module) [1]

--alignment-mode INT How to compute the alignment:

0: automatic

1: only score and end_pos

2: also start_pos and cov

3: also seq.id

4: only ungapped alignment [2]

--alignment-output-mode INT How to compute the alignment:

0: automatic

1: only score and end_pos

2: also start_pos and cov

3: also seq.id

4: only ungapped alignment

5: score only (output) cluster format [0]

--wrapped-scoring BOOL Double the (nucleotide) query sequence during the scoring process to allow wrapped diagonal scoring around end and start [0]

-e DOUBLE List matches below this E-value (range 0.0-inf) [1.000E+01]

--min-seq-id FLOAT List matches above this sequence identity (for clustering) (range 0.0-1.0) [0.000]

--min-aln-len INT Minimum alignment length (range 0-INT_MAX) [30]

--seq-id-mode INT 0: alignment length 1: shorter, 2: longer sequence [0]

--alt-ali INT Show up to this many alternative alignments [0]

-c FLOAT List matches above this fraction of aligned (covered) residues (see --cov-mode) [0.800]

--cov-mode INT 0: coverage of query and target

1: coverage of target

2: coverage of query

3: target seq. length has to be at least x% of query length

4: query seq. length has to be at least x% of target length

5: short seq. needs to be at least x% of the other seq. length [2]

--max-rejected INT Maximum rejected alignments before alignment calculation for a query is stopped [2147483647]

--max-accept INT Maximum accepted alignments before alignment calculation for a query is stopped [2147483647]

--score-bias FLOAT Score bias when computing SW alignment (in bits) [0.000]

--realign BOOL Compute more conservative, shorter alignments (scores and E-values not changed) [0]

--realign-score-bias FLOAT Additional bias when computing realignment [-0.200]

--realign-max-seqs INT Maximum number of results to return in realignment [2147483647]

--corr-score-weight FLOAT Weight of backtrace correlation score that is added to the alignment score [0.000]

--gap-open TWIN Gap open cost [aa:11,nucl:5]

--gap-extend TWIN Gap extension cost [aa:1,nucl:2]

--zdrop INT Maximal allowed difference between score values before alignment is truncated (nucleotide alignment only) [40]

--gap-pc INT Pseudo count for calculating position-specific gap opening penalties [10]

--exhaustive-search-filter INT Filter result during search: 0: do not filter, 1: filter [0]

profile:

--pca Pseudo count admixture strength

--pcb Pseudo counts: Neff at half of maximum admixture (range 0.0-inf)

--mask-profile INT Mask query sequence of profile using tantan [0,1] [1]

--e-profile DOUBLE Include sequences matches with < E-value thr. into the profile (>=0.0) [1.000E-03]

--wg BOOL Use global sequence weighting for profile calculation [0]

--filter-msa INT Filter msa: 0: do not filter, 1: filter [1]

--filter-min-enable INT Only filter MSAs with more than N sequences, 0 always filters [0]

--max-seq-id FLOAT Reduce redundancy of output MSA using max. pairwise sequence identity [0.0,1.0] [0.900]

--qid STR Reduce diversity of output MSAs using min.seq. identity with query sequences [0.0,1.0]

Alternatively, can be a list of multiple thresholds:

E.g.: 0.15,0.30,0.50 to defines filter buckets of ]0.15-0.30] and ]0.30-0.50] [0.0]

--qsc FLOAT Reduce diversity of output MSAs using min. score per aligned residue with query sequences [-50.0,100.0] [-20.000]

--cov FLOAT Filter output MSAs using min. fraction of query residues covered by matched sequences [0.0,1.0] [0.000]

--diff INT Filter MSAs by selecting most diverse set of sequences, keeping at least this many seqs in each MSA block of length 50 [1000]

--pseudo-cnt-mode INT use 0: substitution-matrix or 1: context-specific pseudocounts [0]

--num-iterations INT Number of iterative profile search iterations [1]

--exhaustive-search BOOL For bigger profile DB, run iteratively the search by greedily swapping the search results [0]

--lca-search BOOL Efficient search for LCA candidates [0]

misc:

--taxon-list STR Taxonomy ID, possibly multiple values separated by ','

--rescore-mode INT Rescore diagonals with:

0: Hamming distance

1: local alignment (score only)

2: local alignment

3: global alignment

4: longest alignment fulfilling window quality criterion [0]

--allow-deletion BOOL Allow deletions in a MSA [0]

--min-length INT Minimum codon number in open reading frames [30]

--max-length INT Maximum codon number in open reading frames [32734]

--max-gaps INT Maximum number of codons with gaps or unknown residues before an open reading frame is rejected [2147483647]

--contig-start-mode INT Contig start can be 0: incomplete, 1: complete, 2: both [2]

--contig-end-mode INT Contig end can be 0: incomplete, 1: complete, 2: both [2]

--orf-start-mode INT Orf fragment can be 0: from start to stop, 1: from any to stop, 2: from last encountered start to stop (no start in the middle) [1]

--forward-frames STR Comma-separated list of frames on the forward strand to be extracted [1,2,3]

--reverse-frames STR Comma-separated list of frames on the reverse strand to be extracted [1,2,3]

--translation-table INT 1) CANONICAL, 2) VERT_MITOCHONDRIAL, 3) YEAST_MITOCHONDRIAL, 4) MOLD_MITOCHONDRIAL, 5) INVERT_MITOCHONDRIAL, 6) CILIATE

  1. FLATWORM_MITOCHONDRIAL, 10) EUPLOTID, 11) PROKARYOTE, 12) ALT_YEAST, 13) ASCIDIAN_MITOCHONDRIAL, 14) ALT_FLATWORM_MITOCHONDRIAL

  2. BLEPHARISMA, 16) CHLOROPHYCEAN_MITOCHONDRIAL, 21) TREMATODE_MITOCHONDRIAL, 22) SCENEDESMUS_MITOCHONDRIAL

  3. THRAUSTOCHYTRIUM_MITOCHONDRIAL, 24) PTEROBRANCHIA_MITOCHONDRIAL, 25) GRACILIBACTERIA, 26) PACHYSOLEN, 27) KARYORELICT, 28) CONDYLOSTOMA

  4. MESODINIUM, 30) PERTRICH, 31) BLASTOCRITHIDIA [1]

--translate INT Translate ORF to amino acid [0]

--use-all-table-starts BOOL Use all alternatives for a start codon in the genetic table, if false - only ATG (AUG) [0]

--id-offset INT Numeric ids in index file are offset by this value [0]

--add-orf-stop BOOL Add stop codon '*' at complete start and end [0]

--sequence-overlap INT Overlap between sequences [0]

--sequence-split-mode INT Sequence split mode 0: copy data, 1: soft link data and write new index, [1]

--headers-split-mode INT Header split mode: 0: split position, 1: original header [0]

--search-type INT Search type 0: auto 1: amino acid, 2: translated, 3: nucleotide, 4: translated nucleotide alignment [0]

--start-sens FLOAT Start sensitivity [4.000]

--sens-steps INT Number of search steps performed from --start-sens to -s [1]

--simple-best-hit BOOL Update the p-value by a single best hit, or by best and second best hits [1]

--suboptimal-hits INT Include sub-optimal hits of query sequence up to a factor of its E-value. 0: only include one best hit [0]

--alpha FLOAT Set alpha for combining p-values during aggregation [1.000]

--aggregation-mode INT Combined P-values computed from 0: multi-hit, 1: minimum of all P-values, 2: product-of-P-values, 3: truncated product [0]

--filter-self-match BOOL Remove hits between the same set. 0: do not filter, 1: filter [0]

--multihit-pval FLOAT Multihit P-value threshold for cluster match output [0.010]

--cluster-pval FLOAT Clustering and Ordering P-value threshold for cluster match output [0.010]

--max-gene-gap INT Maximum number of genes allowed between two clusters to merge [3]

--cluster-size INT Minimum number of genes to define cluster [2]

--cluster-use-weight BOOL Use weighting factor to calculate cluster match score [0]

--profile-cluster-search BOOL Perform profile(target)-sequence searches in clustersearch [0]

--search-mode INT 0: sequence search with MMseqs2, 1: structure comparison with Foldseek [0]

common:

--sub-mat TWIN Substitution matrix file [aa:blosum62.out,nucl:nucleotide.out]

--max-seq-len INT Maximum sequence length [65535]

--db-load-mode INT Database preload mode 0: auto, 1: fread, 2: mmap, 3: mmap+touch [0]

--threads INT Number of CPU-cores used (all by default) [128]

--compressed INT Write compressed output [0]

-v INT Verbosity level: 0: quiet, 1: +errors, 2: +warnings, 3: +info [3]

--mpi-runner STR Use MPI on compute cluster with this MPI command (e.g. "mpirun -np 42")

--force-reuse BOOL Reuse tmp filse in tmp/latest folder ignoring parameters and version changes [0]

--remove-tmp-files BOOL Delete temporary files [0]

expert:

--filter-hits BOOL Filter hits by seq.id. and coverage [0]

--sort-results INT Sort results: 0: no sorting, 1: sort by E-value (Alignment) or seq.id. (Hamming) [0]

--create-lookup INT Create database lookup file (can be very large) [0]

--chain-alignments INT Chain overlapping alignments [0]

--merge-query INT Combine ORFs/split sequences to a single entry [1]

--strand INT Strand selection only works for DNA/DNA search 0: reverse, 1: forward, 2: both [1]

--db-output BOOL Return a result DB instead of a text file [1]

references:

- Steinegger M, Soding J: MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology, 35(11), 1026-1028 (2017)

Foldseekデータベース

いくつか利用できる。画像はレポジトリより。

#1 PDBをpdbとしてカレントディレクトリにダウンロードする。
foldseek databases PDB pdb tmpFolder

#2 spacedustの構造DBに変換
spacedust aa2foldseek setDB refFoldseekDB tmpFolder

実行方法

ゲノムのfastaファイルまたはGFF3、あるいはProdigalのタンパク質FASTAファイル(.faaまたは.faa.gz)を指定する。

spacedust clustersearch querySetDB targetSetDB result.tsv tmpFolder

タブ区切りのテキストが出力される。

作成中

引用

De novo discovery of conserved gene clusters in microbial genomes with Spacedust

Ruoshi Zhang, Milot Mirdita, Johannes Söding

bioRxiv, Posted October 03, 2024.

関連

比較ゲノム解析では、ゲノムの遺伝子座のアラインメントを可視化することがよくある。PythonやRのライブラリからスタンドアローンのGUIまで、このタスクのためにいくつかのソフトウェアツールが利用可能であるが、高速で自動化された使用法と出版可能なベクター画像の作成を提供するツールが不足している。

ここでは、LoVis4uを紹介する。LoVis4uは、複数のゲノム遺伝子座を高度にカスタマイズ可能かつ高速に可視化するために設計されたコマンドラインツールとPython APIである。LoVis4uは、GenBankまたはGFFファイルからのアノテーションデータに基づいて、PDF形式のベクター画像を生成する。LoVis4uは、原核生物ゲノムのプラスミドやユーザー定義領域だけでなく、バクテリオファージの全ゲノムを可視化することができる。さらに、LoVis4uは、入力配列中のアクセサリー遺伝子やコア遺伝子を同定し、ハイライトするためのオプションのデータ処理ステップを提供する。

LoVis4uはPython3で実装されており、LinuxとMacOS上で動作する。コマンドラインインターフェースは、ほとんどの実用的なユースケースをカバーし、提供されるPython APIは、Pythonプログラム内での使用、外部ツールへの統合、追加のカスタマイズを可能にする。ソースコードはGitHubページで入手できる：github.com/art-egorov/lovis4u。例によるガイドを含む詳細なドキュメントは、ソフトウェアのホームページから入手できる： art-egorov.github.io/lovis4u。

HP

https://art-egorov.github.io/lovis4u/

Gallery

https://art-egorov.github.io/lovis4u/Gallery/gallery/

インストール

ubuntu22.04にcondaで環境を作ってインストールした。また、WSLのubuntu22でもテストした。

#PyPI ( link )
mamba create -n lovis4u python=3.11 -y
conda activate lovis4u
python3 -m pip install lovis4u

#Linuxマシンを使っている場合、インストール後に'lovis4u --linux` コマンドを実行して、Linuxのmmseqsバイナリに切り替える必要がある
lovis4u --linux

#mmseqsが入ってない場合は導入する
mamba install -c conda-forge -c bioconda mmseqs2 -y

> lovis4u -h

LoVis4u (version 0.0.9):

Home page and documentation: https://github.com/art-egorov/lovis4u

The Atkinson Lab 4U | AE

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

COMMAND-LINE PARAMETERS

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

[POST-INSTALL STEPS]

--data

Creates the 'lovis4u_data' folder in the current working directory.

The folder contains adjustable configuration files used by lovis4u

(e.g. config, palettes...)

--linux

Replaces the mmseqs path in the pre-made config file from the MacOS

version [default] to the Linux version.

--mac

Replaces the mmseqs path in the pre-made config file from the Linux

version [default] to the MacOS version.

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

[MANDATORY ARGUMENTS]

-gff

Path to a folder containing extended gff files.

Each gff file should contain corresponding nucleotide sequence.

(designed to handle pharokka produced annotation files).

OR

-gb

Path to a folder containing genbank files.

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

[OPTIONAL ARGUMENTS | DATA PROCESSING]

-ufid, --use-filename-as-id

Use filename (wo extension) as track (contig) id instead

of the contig id written in the gff/gb file.

-laf, --locus-annotation-file

Path to the locus annotation table.

(See documentation for details)

-faf, --feature-annotation-file

Path to the feature annotation table.

(See documentation for details)

-mmseqs-off, --mmseqs-off

Deactivate mmseqs clustering of proteomes of loci.

-cl-owp, --cluster-only-window-proteins

Cluster only proteins that are overlapped with

the visualisation windows, not all.

-fv-off, --find-variable-off

Deactivate annotation of variable or conserved protein clusters.

-cl-off, --clust_loci-off

Deactivate defining locus order and using similarity based hierarchical

clustering of proteomes.

-oc, --one-cluster

Consider all sequences to be members of one cluster but use clustering

dendrogram to define the optimal order.

-reorient_loci, --reorient_loci

Auto re-orient loci (set new strands) if they are not matched.

(Function tries to maximise co-orientation of homologous features.)

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

[OPTIONAL ARGUMENTS | LOCUS VISUALISATION]

-sgc-off, --set-group-colour-off

Deactivate auto-setting of feature fill and stroke colours.

(Pre-set colours specified in feature annotation table will be kept.)

-sgcf, --set-group-colour-for <feature_group1 [feature group2 ...]>

Space-separated list of feature groups for which colours should be set.

[default: variable, labeled]

-scc, --set-category-colour

Set category colour for features and plot category colour legend.

-cct, --category-colour-table

Path to the table with colour code for categories.

Default table can be found in lovis4u_data folder.

-lls, --locus-label-style <id|description|full>

Locus label style based on input sequence annotation.

-llp, --locus-label-position <left|bottom>

Locus label position on figure.

-safl, --show-all-feature-labels

Display all feature labels.

-sflf, --show-feature-label-for <feature_group1 [feature group2 ...]>

Space-separated list of feature groups for which label should be shown.

[default: variable, labeled]

-sfflf, --show-first-feature-label-for <feature_group1 [feature group2 ...]>

Space-separated list of feature group types for which label will be displayed

only for the first occurrence of feature homologues group.

[default: shell/core]

-ifl, --ignored-feature-labels <feature_label1 [feature_label2 ...]>

Space-separated list of feature names for which label won't be shown.

[default: hypothetical protein, unknown protein]

-sxa, --show-x-axis

Plot individual x-axis for each locus track.

-hix, --hide-x-axis

Do not plot individual x-axis for each locus track.

-dml, --draw-middle-line

Draw middle line for each locus.

-mm-per-nt, --mm-per-nt <float value>

Scale which defines given space for each nt cell on canvas.

[default: 0.05]

-fw, --figure-width <float value>

Output figure width in mm.

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

[OPTIONAL ARGUMENTS | ADDITIONAL TRACKS]

-hl, --homology-links

Draw homology link track.

-slt, --scale-line-track

Draw scale line track.

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

[OPTIONAL ARGUMENTS | OTHERS]

-o

Output dir name. It will be created if it does not exist.

[default: lovis4u_{current_date}; e.g. uorf4u_2022_07_25-20_41]

--pdf-name

Name of the output pdf file (will be saved in the output folder).

[default: lovis4u.pdf]

-c <standard|<file.cfg>

Path to a configuration file or name of a pre-made config file

[default: standard]

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

[MISCELLANEOUS ARGUMENTS]

-h, --help

Show this help message and exit.

-v, --version

Show program version.

--debug

Provide detailed stack trace for debugging purposes.

--parsing-debug

Provide detailed stack trace for debugging purposes

for failed reading of gff/gb files.

-q, --quiet

Don't show progress messages.

テストラン

テスト用のデータをカレントにコピーするオプションが用意されている。

lovis4u --data
cd lovis4u_data/guide/

実行するにはGFFファイルのディレクトリを指定する。

cd lovis4u_data/guide/
lovis4u -gff gff_files/ -hl --set-category-colour -c A4p2

結果はディレクトリに保存される。

lovis4u.pdf

レイアウトは便利なプリセットパラメータが準備されている。-cで指定する。

（gallaryより転載）

単一配列の視覚化。相同タンパク質群を異なる色でハイライトする。

cd lovis4u/lovis4u/lovis4u_data/guide/
lovis4u -gff single_gff_file/ -hl --set-category-colour -c A4p2 --set-group-colour-for conserved

78の大腸菌ファージのBASEL phage collection（Maffei et.al. PLOS Biology）

 lovis4u -gff BaselCollection/  -hl --set-category-colour -c A4p2 -fw 500

（途中まで）

その他

引用

LoVis4u: Locus Visualisation tool for comparative genomics

Artyom A. Egorov, Gemma C. Atkinson

bioRxiv, Posted September 14, 2024.

関連

レポジトリより

これはPCRプライマーを作成するための小さくて簡単なツールである。用途はDNAアセンブリである。Primer3の代わりにプライマーを選択する理由は以下の通り：

インストール

Github

#PyPI
pip install primers

> primers -h

usage: primers [-h] [--version] {score,create} ...

Create or score PCR primers

positional arguments:

{score,create}

score score existing primers

create create primers to amplify a sequence

optional arguments:

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

--version show program's version number and exit

> primers score -h

$ primers score -h

usage: primers score [-h] [-s SEQ] [-t SEQ] [-j | --json | --no-json] primer [primer ...]

positional arguments:

primer primer sequences

optional arguments:

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

-s SEQ the sequence that was amplified

-t SEQ sequence to check for off-target binding sites

-j, --json, --no-json

write the primers to a JSON array

> primers create -h

$ primers create -h

usage: primers create [-h] [-f SEQ] [-fl INT INT] [-r SEQ] [-rl INT INT] [-t SEQ] [-j | --json | --no-json] SEQ

positional arguments:

SEQ create primers to amplify this sequence

optional arguments:

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

-f SEQ additional sequence to add to FWD primer (5' to 3')

-fl INT INT space separated min-max range for the length to add from '-f' (5' to 3')

-r SEQ additional sequence to add to REV primer (5' to 3')

-rl INT INT space separated min-max range for the length to add from '-r' (5' to 3')

-t SEQ sequence to check for off-target binding sites

-j, --json, --no-json

write the primers to a JSON array

実行方法

primersは、長さ、tm、GC比、二次構造、オフターゲット結合を最適化しながらペアを選ぶ。最も単純なケースでは、増幅したい配列を渡すだけでデザインできる。

primers create CTACTAATAGCACACACGGGGACTAGCATCTATCTCAGCTACGATCAGCATC

出力例

設計するプライマーの5'末端を任意の配列にする（ターゲットにはない配列；例えば制限酵素でカットされる配列など）。

primers create CTACTAATAGCACACACGGGGACTAGCATCTATCTCAGCTACGATCAGCATC -f GGG -r CCC

5'末端に付加する最適な部分配列を選択させたい場合、-fl or -rlでmin-maxサイズを指定する。

primers create CTACTAATAGCACACACGGGGACTAGCATCTATCTCAGCTACGATCAGCATC -fl 2 5 -f CCCCC

CLIでは固定のデザインパラメータが適用されますが、pthonではより細かい設定を指定してプライマーを設計する事ができます。レポジトリで簡単な例が説明されています。

ほか、レポジトリより

引用

https://github.com/Lattice-Automation/primers

関連

2024/10/09追記

生のシーケンスデータの処理に特化したバイオインフォマティクスツールの大部分は、k-mersの概念を多用している。これにより、データの冗長性（ひいてはメモリの圧迫）を減らし、シーケンスエラーを破棄し、操作可能で容易に比較できる固定サイズのオブジェクトを処分することができる。欠点は、各k-merとそれが属する元の配列セットとの間のリンクが一般に失われることである。この文脈で考慮されるデータ量を考えると、この関連性を探し出すにはコストがかかる。この研究では、「back_to_sequences 」を紹介する。「back_to_sequences 」は、興味のあるk-merの集合にインデックスを付け、配列の集合をストリーミングし、インデックス付けされたk-merを少なくとも1つ含むものを抽出するように設計されたシンプルなツールである。さらに、配列中のk-merの出現数も提供する。back_to_sequencesは1ミリ秒あたり約200のショートリードをストリームし、数億のリード中のk-merを数分で検索できる。

Documentation

https://b2s-doc.readthedocs.io/en/latest/index.html

My first paper on @JOSS_TheOJ ! https://t.co/kbhzvWh0Xp

— Pierre Marijon 🏳️‍🌈 (@pierre_marijon) 2024年9月24日

Hey kmers' fans!
The "back_to_sequences" tool https://t.co/IgCdQbOOzd is now published in JOSS https://t.co/12CulElV0v.
Find sequences (reads, unitigs, genes) related to a set of kmers in large datasets, in a matter of seconds.
More to read in this thread: https://t.co/Eh7QI33y3w pic.twitter.com/OhK19xJEcF

— Pierre Peterlongo (@p_peterlongo) 2024年9月23日

インストール

ubuntu22.04LTSでテストした（rustc --version => rustc 1.78.0）。

To compile from source

Github

git clone https://github.com/pierrepeterlongo/back_to_sequences.git
cd back_to_sequences
RUSTFLAGS="-C target-cpu=native" cargo install --path .

> back_to_sequences -h

Back to sequences: find the origin of kmers

Usage: back_to_sequences [OPTIONS] --in-kmers <IN_KMERS>

Options:

--in-kmers <IN_KMERS>

Input fasta file containing the original kmers

Note: back_to_sequences considers the content as a set of kmers

This means that a kmer is considered only once,

even if it occurs multiple times in the file.

If the stranded option is not used (default), a kmer

and its reverse complement are considered as the same kmer.

--in-sequences <IN_SEQUENCES>

Input fasta or fastq [.gz] file containing the original sequences (eg. reads).

The stdin is used if not provided

(and if `--in_filelist` is not provided neither) [default: ]

--in-filelist <IN_FILELIST>

Input txt file containing in each line a path to a fasta or fastq [.gz] file

containing the original sequences (eg. reads).

Note1: if this option is used, the `--out_filelist` option must be used.

The number of lines in out_filelist must be the same as in_filelist

Note2: Incompatible with `--in_sequences` [default: ]

--out-sequences <OUT_SEQUENCES>

Output file containing the filtered original sequences (eg. reads).

It will be automatically in fasta or fastq format depending on the input file.

If not provided, only the in_kmers with their count is output [default: ]

--out-filelist <OUT_FILELIST>

Output txt file containing in each line a path to a fasta or fastq [.gz] file

that will contain the related output file from the input files list [default: ]

--out-kmers <OUT_KMERS>

If provided, output a text file containing the kmers that occur in the reads

with their

* number of occurrences

or

* their occurrence positions if the --output_kmer_positions option is used

Note: if `--in_filelist` is used the output counted kmers are

those occurring the last input file of that list [default: ]

--counted-kmer-threshold <COUNTED_KMER_THRESHOLD>

If out_kmers is provided, output only reference kmers whose number of occurrences

is at least equal to this value.

If out_kmers is not provided, this option is ignored [default: 0]

--output-kmer-positions

If out_kmers is provided, either only count their number of occurrences (default)

or output their occurrence positions (read_id, position, strand) if this option is used

--output-mapping-positions

If provided, output matching positions on sequences in the

out_sequence file(s)

-k, --kmer-size <KMER_SIZE>

Size of the kmers to index and search [default: 31]

-m, --min-threshold <MIN_THRESHOLD>

Output sequences are those whose ratio of indexed kmers is in ]min_threshold; max_threshold]

Minimal threshold of the ratio (%) of kmers that must be found in a sequence to keep it (default 0%).

Thus by default, if no kmer is found in a sequence, it is not output. [default: 0]

--max-threshold <MAX_THRESHOLD>

Output sequences are those whose ratio of indexed kmers is in ]min_threshold; max_threshold]

Maximal threshold of the ratio (%) of kmers that must be found in a sequence to keep it (default 100%).

Thus by default, there is no limitation on the maximal number of kmers found in a sequence. [default: 100]

--stranded

Used original kmer strand (else canonical kmers are considered)

--query-reverse

Query the reverse complement of reads. Useless without the --stranded option

--no-low-complexity

Do not index low complexity kmers (ie. with a Shannon entropy < 1.0)

-t, --threads

Number of threads

Note: if not provided, the number of threads is set to the number of logical cores [default: 0]

-h, --help

Print help

-V, --version

Print version

実行方法

探索するk-mer配列のfastaと探索対象のリード、出力として返すリード名、さらに任意でk-merカウント数を示したテキストを指定する。

back_to_sequences --in-kmers compacted_kmers.fasta --in-sequences reads.fasta --out-sequences filtered_reads.fasta  --out-kmers counted_kmers.txt

出力例

filtered_reads.fastaは、入力にfastqフォーマットのリードを指定した時はfastqフォーマットとして出力される。

> head counted_kmers.txt

$ head counted_kmers.txt

CGCTCTATCTGAAACGCGGCGATACGATTTA 1

CGGAATTTGGAGAGAGTATTTTGGTCACGGC 1

ACCGACACGAAGCGACCCTCGCCCTGATCAT 2

CGCAACGTCGCGCGGATGACCGTGCTGCTCC 1

CTGCTCGATCTGTTCTTTGGTGACGGGCTTA 1

CCGCGCCAGCTTGATCGCCATAGTGCCAATC 1

ACATCGCCCACTTCGCGCAGACGAAAGCCCG 1

AACCTTCCGCAAAGGCAAGCGCCATACCGAC 1

ATTTTTCGGTCGCCATTGGCGCAAGCGCTGA 1

CGGCGGTACGTTCACATTCTCAGCGAGCTCC 1

compacted_kmers.fastaに50%以上、多くても70%以下のkmerが含まれる配列のみ出力するフィルターを指定。

back_to_sequences --in-kmers compacted_kmers.fasta --in-sequences reads.fasta --out-sequences filtered_reads.fasta  --out-kmers counted_kmers.txt --min-threshold 50 --max-threshold 70

ほかにもstrand指定や逆相補の配列をターゲットにするオプションなどがある（Document参照）。

論文より

コメント

このプログラムを適切なk値で使うと、関心ある配列をシークエンスしたリードをfastq全部から濃縮することもできると思います（Mirabaitのように）。配列からk-merを発生させるにはjellyfishが使えます。jellyfishは過去に紹介しています（link）。

引用

Back to sequences: Find the origin of 𝑘-mers

Anthony Baire, Pierre Marijon, Francesco Andreace, Pierre Peterlongo

JOSS, Published: 23 September 2024

関連

MetabuliのGUIアプリがリリースされているので簡単に紹介します。

https://github.com/steineggerlab/Metabuli-App/releases/tag/v1.0.0

"これはMetabuli Appの最初のリリースで、これまでコマンドライン経由でのみ利用可能だったMetabuliメタゲノム分類ツールをデスクトップユーザーに提供する。このアプリにより、ユーザーは数回クリックするだけで分類学的プロファイリングを実行したり、過去のレポートをアップロードして素早く視覚化することができる。インタラクティブでカスタマイズ可能なサンキープロットとkronaチャートも含まれており、データの可視化と結果の探索が強化されている。また、特定の分類群に分類されたリードを抽出し、サンキー画像をダウンロードすることもできる。さらに、このアプリには過去10件のジョブ履歴が保存され、すばやくアクセスできるようになっている。"

Our metagenomic classifier Metabuli, which combines AA/DNA information, is now available as app on Windows, MacOS, and Linux. The app developed by @haysunny_hi allows classifying and analyzing metagenomic samples on a notebook! Read more about the app here https://t.co/pfim5rZabc pic.twitter.com/mAldCubH9R

— Martin Steinegger 🇺🇦 (@thesteinegger) October 1, 2024

インストール

M1 macstudioでテストした。

Platforms Supported

Github

https://github.com/steineggerlab/Metabuli-App/releasesよりダウンロードする。

実行方法

上に貼られている動画で使い方は説明されています。ここでは流れだけ簡単に確認します。

ダウンロード後、Applicationsにコピーして実行する。

セキュリティの関係で開けないという警告が出る場合、システム環境設定＝＞セキュリティで”ここまま開く”を選ぶ。

起動した。シングルエンドとペアエンドのショートリード、ロングリードを分析できる。

初回はデータベースをダウンロードする必要がある（MetabuliのDBをダウンロードしているならそれも使える）。

ダウンロードできるDB。目的に合わせて選ぶ。

ここでは以前MetabuliのサイトからダウンロードしたGTDB R214を指定した。出力と最大RAM使用量も指定する。

計算にはしばらく時間がかかる（M1 macstudioだと200MBx2のペアエンドfastqで15秒ほど）。

ここでは”LOAD SAMPLE DATA”の結果を見てみる。

rootから種までアサインされたリードの比率が示されている。

Sankey plot

Configure

図はHTMLや画像として保存できる。

Krona plot

右端のextract readsボタンからはその分類群のfastqを取り出して書き出すことが出来る。

引用

Metabuli: sensitive and specific metagenomic classification via joint analysis of amino acid and DNA

Jaebeom Kim & Martin Steinegger

Nature Methods, Published: 20 May 2024

関連

リファレンスゲノム配列に対するリードのアライメントは、次世代シーケンサー（NGS）技術によって得られたヒト全ゲノムシーケンスデータの解析における重要なステップの1つである。遺伝的変異の臨床的解釈の結果やゲノムワイド関連研究GWASの結果など、その後の解析ステップの質は、アライメントの結果としてリードの位置が正しく特定されるかどうかに依存している。ヒトNGS全ゲノムシーケンスデータの量は常に増加している。世界中で多くのヒトゲノム配列決定プロジェクトが行われており、その結果、配列決定されたヒトゲノムの遺伝子変異に関する大規模なデータベースが作成されている。このような既知の遺伝子変異に関する情報は、例えばゲノムグラフを作成することで、新しい個体について得られたシーケンスデータを解析する際のリードアライメント段階におけるアライメントの質を向上させるために利用することができる。リードを線形リファレンスゲノムにアライメントする既存の方法はアライメント速度が速いが、リードをゲノムグラフにアライメントする方法はゲノムの可変領域において精度が高い。線形リファレンス配列のインデックスに既知の遺伝子変異を考慮したリードアライメント法を開発することで、両手法の利点を組み合わせることができる。

本論文では、minimap2_index_modifierツールを紹介する。このツールは、既知の一塩基バリアントとヒト集団特有の挿入/欠失（indel）を用いて、リファレンスゲノムの修正インデックスを構築することができる。modified minimap2インデックスの使用により、バイオインフォマティクスパイプラインを変更することなく、また計算オーバーヘッドを大幅に追加することなく、バリアントコールの品質が向上する。PrecisionFDA Truth Challenge V2のベンチマークデータ（GRCh38リニアリファレンス（GCA_000001405.15）にアライメントされたHG002ショートリードデータ、パラメータk = 27、w = 14）を用いて、Human Pangenome Reference Consortiumの遺伝的バリアントを用いてインデックスを修正すると、偽陰性の遺伝的バリアントの数が9500以上減少し、偽陽性の数が7000以上減少することが実証された。

~minimap2ツールの主な特徴の1つは、minimizersを使用することである。minimizersとは、与えられた配列内に位置する短い部分文字列（または、そうでなければk-measure）のことである。minimizers（最小化子）という用語は、アライメントの過程で比較する必要のある部分文字列の数を最小化するという考え方に由来する。minimizerは、与えられた配列ウィンドウ内の各k-merに与えられたハッシュ関数を適用し、ユニークな整数値を作成することで計算される。次に整数値をソートし、最小値を持つk-measureを最小化子として選択する（論文図2）。ミニマイザーを使用することで、アライメント処理が高速になる。

線形リファレンス配列にインデックスを付ける場合、ハッシュテーブルが作成され、キーは最小化子配列にハッシュ関数を適用した結果であり、値はゲノム参照配列の位置のリストである。

minimap2ツールのインデックスを変更する可能性は、ハッシュテーブルが線形リファレンス配列の所与の位置における最小化因子の数に制約を課さないという事実によって提供され、遺伝的バリアントに関する情報を追加しても、その後のアラインメントアルゴリズムには影響しない。言い換えれば、線形リファレンス配列インデックスは、ゲノムグラフと同様に、遺伝的バリアントの追加によって誘導される分岐の追加を可能にする（論文図3）。（以下省略）

インストール

依存

To compile from source, use this version of tools:

Github

#HTSlibが必要（samtools解説）先にminimap2_index_modifierをcloneして、minimap2_index_modifierのrootで実行する
git clone --recurse-submodules https://github.com/samtools/htslib.git
cd htslib/
autoheader
autoreconf -i
./configure --prefix=/usr/local
make -j10
make install
cd ../

#本体(minimap2のfolk)*1（注;特別なoptionを使わない限り本家minimap2と同じものだが、パスを通す時は本家と混同しないように注意）
git clone https://github.com/ispras/minimap2_index_modifier.git
cd minimap2_index_modifier/
make -j10

#docker image(hub)
docker pull egorguga/minimap2_index_modifier:2.24

> ./minimap2

Usage: minimap2 [options] <target.fa>|<target.idx> [query.fa] [...]

Options:

Indexing:

-H use homopolymer-compressed k-mer (preferrable for PacBio)

-k INT k-mer size (no larger than 28) [15]

-w INT minimizer window size [10]

-I NUM split index for every ~NUM input bases [4G]

-d FILE dump index to FILE

--vcf-file-with-variants FILE pass VCF FILE to modify index

--parse-haplotype parse haplotypes from VCF to generate real SNP combinations, otherwise use all.

Mapping:

-f FLOAT filter out top FLOAT fraction of repetitive minimizers [0.0002]

-g NUM stop chain enlongation if there are no minimizers in INT-bp [5000]

-G NUM max intron length (effective with -xsplice; changing -r) [200k]

-F NUM max fragment length (effective with -xsr or in the fragment mode) [800]

-r NUM[,NUM] chaining/alignment bandwidth and long-join bandwidth [500,20000]

-n INT minimal number of minimizers on a chain [3]

-m INT minimal chaining score (matching bases minus log gap penalty) [40]

-X skip self and dual mappings (for the all-vs-all mode)

-p FLOAT min secondary-to-primary score ratio [0.8]

-N INT retain at most INT secondary alignments [5]

Alignment:

-A INT matching score [2]

-B INT mismatch penalty (larger value for lower divergence) [4]

-O INT[,INT] gap open penalty [4,24]

-E INT[,INT] gap extension penalty; a k-long gap costs min{O1+k*E1,O2+k*E2} [2,1]

-z INT[,INT] Z-drop score and inversion Z-drop score [400,200]

-s INT minimal peak DP alignment score [80]

-u CHAR how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n]

Input/Output:

-a output in the SAM format (PAF by default)

-o FILE output alignments to FILE [stdout]

-L write CIGAR with >65535 ops at the CG tag

-R STR SAM read group line in a format like '@RG\tID:foo\tSM:bar'

-c output CIGAR in PAF

--cs[=STR] output the cs tag; STR is 'short' (if absent) or 'long' [none]

--MD output the MD tag

--eqx write =/X CIGAR operators

-Y use soft clipping for supplementary alignments

-t INT number of threads [3]

-K NUM minibatch size for mapping [500M]

--version show version number

Preset:

-x STR preset (always applied before other options; see minimap2.1 for details)

- map-pb/map-ont - PacBio CLR/Nanopore vs reference mapping

- map-hifi - PacBio HiFi reads vs reference mapping

- ava-pb/ava-ont - PacBio/Nanopore read overlap

- asm5/asm10/asm20 - asm-to-ref mapping, for ~0.1/1/5% sequence divergence

- splice/splice:hq - long-read/Pacbio-CCS spliced alignment

- sr - genomic short-read mapping

See `man ./minimap2.1' for detailed description of these and other advanced command-line options.

ビルド済みindex

２つのヒトゲノムアセンブリバージョンでのmodified indexが利用できる。

Nextcloud

テストラン

test/README.mdの通り進めます。

1、indexing

cd minimap2_index_modifier/tests/
#比較のため通常のindexとバリアントを反映させたindexを作る
#A normal
../minimap2 -ax sr -d test.mni test.fasta

#B modified index (このminimap2 folkでないと実行できない点に注意)
../minimap2 -ax sr -d test.modified.mni --vcf-file-with-variants test_long_chr1_not_the_same.vcf.gz test.fasta

=>diffコマンドで異なるか確認
diff test.mni test.modified.mni

２、mapping

../minimap2 -ax sr -t 4 test.modified.mni test.1.fq test.2.fq > mapped.sam

アラインメントフォーマットには変化はないが、バリアントを含む部位でのアラインメントが改善されている可能性がある。

レポジトリと論文より

コメント

パイプラインに影響を与えないため、最小限の手間でバリアントコールパイプラインを修正することが可能になります。気になるのは、図７のグラフで、わずかな差ではありますが、デフォルトindexのminimap2はrecallの指標でbwa mem2に負けていて、ラージコホートのバリアントを使ったmodified indexを使ったminimap2はbwa mem2を上回る点でしょうか。当時minimap2がショートリードマッピングでわずかに劣る可能性が指摘されてましたが、それが見えているのかと思いました

引用

Enhancing SNV identification in whole-genome sequencing data through the incorporation of known genetic variants into the minimap2 index

Egor Guguchkin, Artem Kasianov, Maksim Belenikin, Gaukhar Zobkova, Ekaterina Kosova, Vsevolod Makeev & Evgeny Karpulevich

BMC Bioinformatics volume 25, Article number: 238 (2024)

*1

デフォルトでは"../htslib"にある HTSlibに対してビルドされるので、取ってきたminimap2_index_modifierの直下でHTSlibライブラリをビルドした。

関連

微生物群集に含まれる細菌種は、ゲノムの小さな変異によって区別される菌株の混合物であることが多い。ショートリード法は、菌株間の小規模な変異を検出するために使用できるが、これらの変異を連続したハプロタイプにphasing（位相を揃える）することはできない。ロングリードメタゲノムアセンブラは、連続した細菌chromosomeを作成できるが、しばしば、種レベルのコンセンサスを優先して菌株レベルの変異を抑制する。ここでは、NanoporeおよびPacBioリードから株レベルのメタゲノム・アセンブリとphasingを行うアルゴリズム、Strainyを紹介する。Strainyでは、de novoメタゲノム・アセンブリを入力として、株のバリアントを同定し、それをフェーズ処理し、連続したハプロタイプにアセンブルする。NanoporeおよびPacBioメタゲノム・データのシミュレーションおよび模擬データを用いて、Strainyが正確で完全な菌株ハプロタイプをアセンブルすることを示し、その完全性と正確性において、現行のNanoporeベースの手法を上回り、PacBioベースのアルゴリズムと同等であることを示す。次に、Strainyを使用して、複雑な環境メタゲノムの菌株ハプロタイプをアセンブルし、細菌種における明確な菌株分布と変異パターンを明らかにした。

Tutorial

https://github.com/katerinakazantseva/strainy?tab=readme-ov-file#strainy-tutorial

インストール

Github

git clone https://github.com/katerinakazantseva/strainy
cd strainy
git submodule update --init
make -C submodules/Flye
mamba env create -f environment.yml -n strainy
conda activate strainy
./strainy.py -h

#docker imageも公開されている
docker pull mkolmogo/strainy:1.0

>./strainy.py -h

usage: strainy.py [-h] -o OUTPUT -g GFA -m {hifi,nano} -q FASTQ [-s {phase,transform,e2e}] [--snp SNP] [-t THREADS] [-b BAM] [--link-simplify] [--debug] [--unitig-split-length UNITIG_SPLIT_LENGTH] [--only-split ONLY_SPLIT] [-d CLUSTER_DIVERGENCE] [-a ALLELE_FREQUENCY]

[--min-unitig-length MIN_UNITIG_LENGTH] [--min-unitig-coverage MIN_UNITIG_COVERAGE] [--max-unitig-coverage MAX_UNITIG_COVERAGE] [-v]

options:

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

-s {phase,transform,e2e}, --stage {phase,transform,e2e}

stage to run: either phase, transform or e2e (phase + transform) (default: e2e)

--snp SNP path to vcf file with SNP calls to use (default: None)

-t THREADS, --threads THREADS

number of threads to use (default: 4)

-b BAM, --bam BAM path to indexed alignment in bam format (default: None)

--link-simplify Enable agressive graph simplification (default: False)

--debug Generate extra output for debugging (default: False)

--unitig-split-length UNITIG_SPLIT_LENGTH

The length (in kb) which the unitigs that are longer will be split, set 0 to disable (default: 50)

--only-split ONLY_SPLIT

Do not run stRainy, only split long gfa unitigs (default: False)

-d CLUSTER_DIVERGENCE, --cluster-divergence CLUSTER_DIVERGENCE

cluster divergence (default: 0)

-a ALLELE_FREQUENCY, --allele-frequency ALLELE_FREQUENCY

Set allele frequency for internal caller only (pileup) (default: 0.2)

--min-unitig-length MIN_UNITIG_LENGTH

The length (in kb) which the unitigs that are shorter will not be phased (default: 1)

--min-unitig-coverage MIN_UNITIG_COVERAGE

The minimum coverage threshold for phasing unitigs, unitigs with lower coverage will not be phased (default: 20)

--max-unitig-coverage MAX_UNITIG_COVERAGE

The maximum coverage threshold for phasing unitigs, unitigs with higher coverage will not be phased (default: 500)

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

Required named arguments:

-o OUTPUT, --output OUTPUT

output directory (default: None)

-g GFA, --gfa GFA input gfa to uncollapse (default: None)

-m {hifi,nano}, --mode {hifi,nano}

type of reads (default: None)

-q FASTQ, --fastq FASTQ

fastq file with reads to phase / assemble (default: None)

テストラン

StrainyはPacBio HiFi、Nanopore R9（Guppy5+）、R10シーケンサーをサポートしている。

（metaFlyeの）GFAファイルとアセンブル/phasingのリードを指定する。

cd strainy/
./strainy.py -g test_set/toy.gfa -q test_set/toy.fastq.gz -o out_strainy -m nano -t 8

出力例

出力についてはレポジトリのOutput filesで説明されています。

引用

Strainy: phasing and assembly of strain haplotypes from long-read metagenome sequencing

Ekaterina Kazantseva, Ataberk Donmez, Maria Frolova, Mihai Pop & Mikhail Kolmogorov

Nature Methods (2024), Published: 26 September 2024

関連