macでインフォマティクス

macでインフォマティクス

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

リモートヌクレオチド配列を高速に検索して大規模MSAを構築する Riboseek

 

 構造予測の対象はRNAへと拡大したが、現在、深層アラインメントの生成がボトルネックとなっている。本著者らは、配列を重複するダイマー(2塩基単位)として表現する、RNAおよびDNA向けの検索・アラインメントツールRiboseekを開発した。RiboseekはBLASTNやnhmmerよりも高い感度で相同性を検出でき、アラインメント生成速度はnhmmerやrMSAと比較してそれぞれ250倍、376倍高速である。構造を考慮した再アラインメントを行うことで、rMSAと同等の塩基対復元精度を実現しつつ、処理速​​度は29倍高速に保たれている。また、構造を持つと推定されるRNAについて、ウェブサーバーおよび170万件の事前計算済みアラインメントを提供している。Riboseekは https://github.com/steineggerlab/riboseekにて無償で利用できる。

 

インストール

Win11 WSLでテストした。

Github

# Linux AVX2 build (check using: cat /proc/cpuinfo | grep avx2)
wget https://mmseqs.com/riboseek/riboseek-linux-avx2.tar.gz; tar xvzf riboseek-linux-avx2.tar.gz; export PATH=$(pwd)/riboseek/bin/:$PATH

# Linux AVX2 & GPU build (req. glibc >= 2.17 and nvidia driver >= 525.60.13)
wget https://mmseqs.com/riboseek/riboseek-linux-gpu.tar.gz; tar xvzf riboseek-linux-gpu.tar.gz; export PATH=$(pwd)/riboseek/bin/:$PATH

#conda
mamba create -n riboseek -c conda-forge -c bioconda riboseek -y
conda activate riboseek

> riboseek

Riboseek enables fast RNA homology search using covariance models and structure-aware alignment.

 

riboseek Version: 1.0.0

© 

 

usage: riboseek <command> [<args>]

 

Easy workflows for plain text input/output

  easy-search       Sensitive homology search

  easy-cluster      Slower, sensitive clustering

  easy-linclust     Fast linear time cluster, less sensitive clustering

  easy-proteomecluster Cluster proteomes and identify reference proteomes

  easy-proteomesearch Calculate pairwise proteome similarity

  easy-taxonomy     Taxonomic classification

  easy-rbh          Find reciprocal best hit

 

Main workflows for database input/output

  search            RNA homology search with dinucleotide profiles

  map               Map nearly identical sequences

  rbh               Reciprocal best hit search

  linclust          Fast, less sensitive clustering

  cluster           Slower, sensitive clustering

  clusterupdate     Update previous clustering with new sequences

  taxonomy          Taxonomic classification

 

Input database creation

  databases         List and download databases

  createdb          Convert FASTA/Q file(s) to a sequence DB

  createindex       Store precomputed index on disk to reduce search overhead

  convertmsa        Convert Stockholm/PFAM MSA file to a MSA DB

  msa2profile       Convert a MSA DB to a profile DB

 

Format conversion for downstream processing

  convertalis       Convert alignment DB to BLAST-tab, SAM or custom format

  createtsv         Convert result DB to tab-separated flat file

  convert2fasta     Convert sequence DB to FASTA format

  taxonomyreport    Create a taxonomy report in Kraken or Krona format

 

An extended list of all modules can be obtained by calling 'riboseek -h'.

 

> riboseek easy-search -h

usage: riboseek easy-search <i:queryFastaFile1[.gz|.bz2]> ... <i:queryFastaFileN[.gz|.bz2]>|<i:stdin> <i:targetFastaFile[.gz]>|<i:targetDB> <o:alignmentFile> <tmpDir> [options]

options:                                

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

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

                               

 --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 [3]

 --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]

 --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) [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]

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

 --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]

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

                               

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

                               

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

 --prefilter-mode INT            prefilter mode: 0: kmer/ungapped 1: ungapped, 2: nofilter, 3: ungapped&gapped [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]

 --translation-mode INT          Translation AA seq from nucleotide by 0: ORFs, 1: full reading frames [0]

 --format-mode INT               Output format:

                                 0: BLAST-TAB

                                 1: SAM

                                 2: BLAST-TAB + query/db length

                                 3: Pretty HTML

                                 4: BLAST-TAB + column headers

                                 BLAST-TAB (0) and BLAST-TAB + column headers (4) support custom output formats (--format-output) [0]

 --format-output STR             Choose comma separated list of output columns from: query,target,evalue,gapopen,pident,fident,nident,qstart,qend,qlen

                                 tstart,tend,tlen,alnlen,raw,bits,cigar,qseq,tseq,qheader,theader,qaln,taln,qframe,tframe,mismatch,qcov,tcov

                                 qset,qsetid,tset,tsetid,taxid,taxname,taxlineage,qorfstart,qorfend,torfstart,torfend,ppos [query,target,fident,alnlen,mismatch,gapopen,qstart,qend,tstart,tend,evalue,bits]

 --greedy-best-hits BOOL         Choose the best hits greedily to cover the query [0]

                               

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

 --compressed INT                Write compressed output [0]

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

 --gpu-server INT                Use GPU server [0]

 --gpu-server-wait-timeout INT   Wait for GPU server for 0: don't wait -1: no wait limit: >0 this many seconds [600]

 

examples:

 # Search multiple FASTA against FASTA (like BLASTP, TBLASTN, BLASTX, BLASTN --search-type 3, TBLASTX --search-type 2)

 mmseqs easy-search examples/QUERY.fasta examples/QUERY.fasta examples/DB.fasta result.m8 tmp

 

 # Iterative profile search from stdin (like PSI-BLAST)

 cat examples/QUERY.fasta | mmseqs easy-search stdin examples/DB.fasta result.m8 tmp --num-iterations 2

 

 # Profile search against small databases (e.g. PFAM, eggNOG)

 mmseqs databases PFAM pfam_db tmp

 mmseqs easy-search examples/QUERY.fasta pfam_db res.m8 tmp

 

 # Exhaustive search against sequences or profiles (works for large DBs)

 mmseqs easy-search examples/QUERY.fasta targetProfiles res.m8 tmp --exhaustive-search

 

 # Increasing sensitivity search (from 2 to 7 in 3 steps)

 mmseqs easy-search examples/QUERY.fasta examples/DB.fasta result.m8 tmp --start-sens 2 -s 7 --sens-steps 3

 

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)

 - Mirdita M, Steinegger M, Soding J: MMseqs2 desktop and local web server app for fast, interactive sequence searches. Bioinformatics, 35(16), 2856-2858 (2019)

 - Kallenborn F, Chacon A, Hundt C, Sirelkhatim H, Didi K, Cha S, Dallago C, Mirdita M, Schmidt B, Steinegger M: GPU-accelerated homology search with MMseqs2. Nature Methods (2025)

 

Show an extended list of options by calling 'riboseek easy-search -h'.

> riboseek search -h

usage: riboseek search <i:queryDB> <i:targetDB> <o:resultDB> <tmpDir> [options]

options:                                

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

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

                               

 -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 [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]

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

 --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]

                               

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

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

                               

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

 --prefilter-mode INT            prefilter mode: 0: kmer/ungapped 1: ungapped, 2: nofilter, 3: ungapped&gapped [0]

 --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]

 --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]

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

 --start-sens FLOAT              Start sensitivity [4.000]

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

 --translation-mode INT          Translation AA seq from nucleotide by 0: ORFs, 1: full reading frames [0]

                               

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

 --compressed INT                Write compressed output [0]

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

 --gpu INT                       Use GPU (CUDA) if possible [0]

 --gpu-server INT                Use GPU server [0]

 --gpu-server-wait-timeout INT   Wait for GPU server for 0: don't wait -1: no wait limit: >0 this many seconds [600]

 

examples:

 # Search RNA query against RNA target database

 riboseek search queryDB targetDB resultDB tmp

 

 # Iterative profile search

 riboseek search queryDB targetDB resultDB tmp --num-iterations 3

 

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)

 

Show an extended list of options by calling 'riboseek search -h'.

 

> riboseek map -h

usage: riboseek map <i:queryDB> <i:targetDB> <o:alignmentDB> <tmpDir> [options]

 By Milot Mirdita <milot@mirdita.de> & Martin Steinegger <martin.steinegger@snu.ac.kr>

options: prefilter:                  

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

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

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

 --target-search-mode INT     target search mode (0: regular k-mer, 1: similar k-mer) [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:25,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) [0]

 --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 prefilter stage with tantan: 0: w/o low complexity masking, 1: with low complexity masking [0]

 --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]

 --mask-n-repeat INT          Repeat letters that occur > threshold in a rwo [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.950]

 --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]

 --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]

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

 --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]

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 [2]

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

 --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

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

                              15) BLEPHARISMA, 16) CHLOROPHYCEAN_MITOCHONDRIAL, 21) TREMATODE_MITOCHONDRIAL, 22) SCENEDESMUS_MITOCHONDRIAL

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

                               29) 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]

 --start-sens FLOAT           Start sensitivity [4.000]

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

common:                     

 --sub-mat TWIN               Substitution matrix file [aa:dinuc.out,nucl:dinuc.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) [10]

 --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) [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)

 - Steinegger M, Mirdita M, Soding J: Protein-level assembly increases protein sequence recovery from metagenomic samples manyfold. Nature Methods, 16(7), 603-606 (2019)

> riboseek rbh -h

sage: riboseek rbh <i:queryDB> <i:targetDB> <o:alignmentDB> <tmpDir> [options]

 By Eli Levy Karin

options: prefilter:                      

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

 --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:dinuc.out,nucl:dinuc.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]

 --target-search-mode INT         target search mode (0: regular k-mer, 1: similar k-mer) [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:25,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 prefilter stage with tantan: 0: w/o low complexity masking, 1: with low complexity masking [0]

 --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]

 --mask-n-repeat INT              Repeat letters that occur > threshold in a rwo [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) [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 [3]

 --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.000]

 --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]

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

 --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]

 --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]

 --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]

 --profile-output-mode INT        Profile output mode: 0: binary log-odds 1: human-readable frequencies [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 ','

 --prefilter-mode INT             prefilter mode: 0: kmer/ungapped 1: ungapped, 2: nofilter, 3: ungapped&gapped [0]

 --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 [10]

 --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

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

                                  15) BLEPHARISMA, 16) CHLOROPHYCEAN_MITOCHONDRIAL, 21) TREMATODE_MITOCHONDRIAL, 22) SCENEDESMUS_MITOCHONDRIAL

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

                                   29) 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]

 --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]

 --translation-mode INT           Translation AA seq from nucleotide by 0: ORFs, 1: full reading frames [0]

common:                         

 --sub-mat TWIN                   Substitution matrix file [aa:dinuc.out,nucl:dinuc.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) [10]

 --compressed INT                 Write compressed output [0]

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

 --gpu INT                        Use GPU (CUDA) if possible [0]

 --gpu-server INT                 Use GPU server [0]

 --gpu-server-wait-timeout INT    Wait for GPU server for 0: don't wait -1: no wait limit: >0 this many seconds [600]

 --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]

 

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)

> riboseek cluster -h

usage: riboseek cluster <i:sequenceDB> <o:clusterDB> <tmpDir> [options]

 By Martin Steinegger <martin.steinegger@snu.ac.kr> & Lars von den Driesch

options: prefilter:                      

 --seed-sub-mat TWIN              Substitution matrix file for k-mer generation [aa:dinuc.out,nucl:dinuc.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]

 --target-search-mode INT         target search mode (0: regular k-mer, 1: similar k-mer) [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:25,nucl:5]

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

 --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 prefilter stage with tantan: 0: w/o low complexity masking, 1: with low complexity masking [0]

 --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]

 --mask-n-repeat INT              Repeat letters that occur > threshold in a rwo [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 [3]

 --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.000]

 --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]

 --set-mode BOOL                  0: Cluster by each entry

                                  1: Cluster by set [0]

 --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]

 --cluster-version INT            Cluster version: 1: Cluster1, 2: Cluster2 [1]

 --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]

 --include-count-table BOOL       Include count table [1]

 --num-count-table INT            Number of count table based center swapping [2]

 --include-adjacency BOOL         Include adjacency based center swapping [1]

 --num-adjacency INT              Number of adjacency based center swapping [3]

 --use-parallelism BOOL           Enable or disable parallel execution for group assignment and related k-mer processing steps [0]

 --linclust-version INT           Linclust version: 1: Linclust1, 2: Linclust2 [2]

 --clust-hash BOOL                Use clusthash before kmermatcher in linclust [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]

common:                         

 --sub-mat TWIN                   Substitution matrix file [aa:dinuc.out,nucl:dinuc.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) [10]

 --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]

 --filter-cludb-file STR          Specify a cluDB file to filter

 --filter-seqdb-file STR          Specify a seqDB of cluDB file to filter

 

examples:

 # Cascaded clustering of FASTA file

 mmseqs cluster sequenceDB clusterDB tmp

 

 #                  --cov-mode 

 # Sequence         0    1    2

 # Q: MAVGTACRPA  60%  IGN  60%

 # T: -AVGTAC---  60% 100%  IGN

 # Cutoff -c 0.7    -    +    -

 #        -c 0.6    +    +    +

 

 # Cascaded clustering with reassignment

 # - Corrects criteria-violations of cascaded merging

 # - Produces more clusters and is a bit slower

 mmseqs cluster sequenceDB clusterDB tmp --cluster-reassign

 

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)

 - Hauser M, Steinegger M, Soding J: MMseqs software suite for fast and deep clustering and searching of large protein sequence sets. Bioinformatics, 32(9), 1323-1330 (2016)

 - Steinegger M, Soding J: Clustering huge protein sequence sets in linear time. Nature Communications, 9(1), 2542 (2018)

> riboseek linclust -h

usage: riboseek linclust <i:sequenceDB> <o:clusterDB> <tmpDir> [options]

 By Martin Steinegger <martin.steinegger@snu.ac.kr>

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]

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

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

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

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

 --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]

 --mask-n-repeat INT              Repeat letters that occur > threshold in a rwo [0]

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

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

align:                          

 -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 [3]

 --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]

 -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]

 --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]

 --set-mode BOOL                  0: Cluster by each entry

                                  1: Cluster by set [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]

 --include-count-table BOOL       Include count table [1]

 --num-count-table INT            Number of count table based center swapping [2]

 --include-adjacency BOOL         Include adjacency based center swapping [1]

 --num-adjacency INT              Number of adjacency based center swapping [3]

 --use-parallelism BOOL           Enable or disable parallel execution for group assignment and related k-mer processing steps [0]

 --linclust-version INT           Linclust version: 1: Linclust1, 2: Linclust2 [2]

 --clust-hash BOOL                Use clusthash before kmermatcher in linclust [0]

profile:                        

 --pca                            Pseudo count admixture strength

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

misc:                           

 --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]

common:                         

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

 --compressed INT                 Write compressed output [0]

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

 --sub-mat TWIN                   Substitution matrix file [aa:dinuc.out,nucl:dinuc.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]

 --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]

 --filter-cludb-file STR          Specify a cluDB file to filter

 --filter-seqdb-file STR          Specify a seqDB of cluDB file to filter

 

examples:

 # Linear-time clustering of FASTA file

 mmseqs linclust sequenceDB clusterDB tmp

 

                    --cov-mode 

 # Sequence         0    1    2

 # Q: MAVGTACRPA  60%  IGN  60%

 # T: -AVGTAC---  60% 100%  IGN

 # Cutoff -c 0.7    -    +    -

 #        -c 0.6    +    +    +

 

 # Cluster nucleotide sequences 

他のコマンドは割愛。1記事あたりの文字数制限(5万文字以内)があるため。

 

実行方法

QUERY.fastaに含まれる配列を、DB.fastaの配列データベースに対して検索する。

riboseek easy-search query.fasta DB.fasta output_aln.m8 tmp

easy-searchは、入力されたFASTAから内部データベースを自動作成して、そのまま検索、読みやすい表形式出力まで行う便利なコマンドとなっている。

 

NCBI ntのように、同じDB.fastaを何度も検索する場合は、riboseek createdbコマンドを使って事前にデータベースを作成しておく方が効率的となる。

#1.DB作成
riboseek createdb DB.fasta target_db
riboseek createindex target_db tmp_index

#2.ここではクエリもDB作成
riboseek createdb query.fasta query_db

#3.Search
riboseek search query_db target_db search_result tmp_search -a

#4.検索結果を表形式に変換、convertalisコマンドを使う
riboseek convertalis query_db target_db search_result output.m8

#5.MSAに変換。 result2msaコマンドを使う
riboseek result2msa query_db target_db search_result msa_db --msa-format-mode 6

#6.MSAファイルを取り出す、unpackdbコマンドを使う
riboseek unpackdb msa_db msa_a3m --unpack-suffix a3m --unpack-name-mode 0

 

このあとも、通常検索の結果から、RNAの二次構造情報を考慮した再検索を行うこともできる。この場合、MSAを展開したmsa_a3mではなく、3で作成した検索結果データベースsearch_resultを使う。

#7.riboseek cmbuildで検索結果から共分散モデル(CM)を構築
riboseek cmbuild query_db target_db search_result cm_model --cmlite-msa-eval 1e-3 --threads 4

#8.riboseek cmsearchを使い、作成したCMから候補配列を再検索
riboseek cmsearch cm_model target_db search_result cm_search_result --cm-region 3.0 --threads 4

#9.CM検索の結果からMSAを作成。result2msaコマンドを使う
riboseek result2msa query_db target_db cm_search_result cm_msa_db --msa-format-mode 6

#10.MSAファイルを取り出す、unpackdbコマンドを使う
riboseek unpackdb cm_msa_db cm_msa_a3m --unpack-suffix a3m --unpack-name-mode 0

#11.CM検索の結果を表形式に変換
riboseek convertalis query_db target_db cm_search_result output_cm.m8

 

テストラン

クエリの配列QUERY.fastaを、配列データベースDB.fastaに対してサーチする。ここではmmseq系でお馴染みのeasy-searchコマンドを使う。easy-searchは、FASTA配列を入力して、相同性検索と検索結果の表形式への変換までを一括で行う。

git clone https://github.com/steineggerlab/riboseek.git
cd riboseek/
riboseek easy-search query.fasta example/DB.fasta output_aln.m8 tmp

 

出力例

出力ファイルには、querとtargetの配列名、配列一致率、アラインメント長、ミスマッチ数、ギャップ数が記録される。続いて、queryとtarget上の一致位置、E-value、bit scoreが出力される。

(codexより) カラムについてはMMseqのマニュアルで説明されている(リンク)。

 

helpをみると、easyコマンドは以下のような種類が用意されている。

  • easy-search:配列データの読み込みから、データベース作成、相同性検索、検索結果の表形式出力まで一括して行う。
  • easy-cluster:配列データの読み込みから、感度重視の配列クラスタリング、クラスタ結果の出力まで行う。
  • easy-linclust:配列データの読み込みから、高速な配列クラスタリング、クラスタ結果の出力まで行う。
  • easy-proteomecluster:複数のプロテオーム配列データの読み込みから、プロテオーム間のクラスタリングと代表プロテオームの選定まで行う。
  • easy-proteomesearch:プロテオーム配列データの読み込みから、プロテオーム間の類似性検索と比較結果の出力まで行う。
  • easy-taxonomy:配列データの読み込みから、相同性検索、分類群の割り当て、分類結果のレポート作成まで行う。
  • easy-rbh:2つの配列データセットの検索から、相互ベストヒット(RBH)の抽出まで行う。

RiboseekはMMseqs2を基盤として開発されており、MMseqs2由来の汎用コマンドも同じ実行ファイルに含まれている。easy-proteomesearchやeasy-proteomeclusterが表示さているが、Riboseekの中心機能がタンパク質向けという意味ではない。RNAを中心としたヌクレオチド配列向けのツールとなる。

 

その他

  • ヌクレオチド相同性検索に使われるBLASTNのような高速なツールは、ヒューリスティックなワードヒットのシード化と伸長に依存しているため、その感度はシード長、伸長パラメータ、および十分に保存された局所的な一致領域が存在するかどうかに左右される(4、5)。その結果、置換、挿入・欠失、あるいは配列マスキングによってシードとして利用できる領域が少なくなると、遠縁のヌクレオチドホモログが見逃される可能性がある。
  • ヌクレオチド相同性検索の別の手法として、単一のクエリ配列を配列ファミリーの統計モデルに置き換えることによって感度を高めるプロファイルベースの手法がある。現在の標準的なツールであるnhmmerは、この枠組みをDNA/RNAの相同性検索へ拡張し、遠縁のヌクレオチドホモログの検出能力を向上させた(6)。しかし、プロファイル隠れマルコフモデル(profile HMM)を用いるヌクレオチド検索は、最速のシードベース手法よりも2桁以上遅く、大規模なデータベースや多数のクエリへの拡張性が制限されている。
  • 高感度側の別のアプローチは、多くのRNAファミリーでは相補的置換を通じて二次構造が保存されるため、塩基対形成相互作用を維持したまま一次配列が多様化し得るという事実に基づいている(7、8)。共分散モデルは、配列保存性と構造的共変動を単一の確率的枠組みのなかで同時にモデル化することにより、このシグナルを捉える(7)。Infernal(9)に実装されたことで、共分散モデルはRfam(10)などのリソースにおける高感度なRNAファミリー検出の標準手法となった。Infernalのフィルターパイプラインと、HMMによってバンド制限された共分散モデルアラインメントは、共分散モデル検索の速度を大幅に向上させたが、それでもなお計算負荷は非常に高い。
  • RNA構造予測では、多数の相同配列を使った深いMSAが精度に関わる。現在では、RNA構造予測における深いアラインメントの作成が大きなボトルネックになっている (アラインメントの深さがモデルの精度に影響する)。
  • AlphaFold 3(11)とProtenix-v2(12)は、いずれもRNA MSAを入力として受け付ける。AlphaFold 3のパイプラインでは、nhmmerを用いてRNAcentral(13)およびNT(14)を検索することによりRNA MSAを構築する。さらに低速な方法としてrMSA(15)があり、nhmmerおよびBLASTNによる検索に、Rfamに基づくフィルタリングと構造を考慮したフィルタリングを組み合わせている。大規模データベースを用いたMSA生成はクエリごとに実行され、数時間を要し主要なボトルネックとなっている。
  • MMseqs2(16)のようなツールは、タンパク質相同性検索における速度と感度のトレードオフを大きく変え、AlphaFold2(17)やColabFold(18)によるタンパク質構造予測を含むさまざまな応用を可能にした。しかし、ヌクレオチド配列については、これに相当する進歩は達成されていない。本研究では、プロファイルベースのダイマー採点、高速なプレフィルタリング、およびギャップを含むアラインメントを組み合わせ、大規模に相同配列を同定する配列ベースの検索エンジンRiboseekを提示する(論文図1)。
  • 本研究のベンチマーク全体において、RiboseekはBLASTN、nhmmer、およびペアワイズSmith–Waterman検索よりも高い感度を示しながら、nhmmerおよびInfernalよりも2桁高速だった。
  • Riboseekで作成したMSAをProtenix-v2に入力すると、nhmmerのMSAを使った場合より、構造予測の指標が改善した。また、約173万本のRNAについて事前計算済みMSAを作成し、Web APIやColabFoldから利用できるようにした。
  • (レポジトリより) Riboseekは、大規模なRNA配列セットに対する高速かつ高感度な検索を可能にする。本ツールは、事前フィルタリングにジヌクレオチド(2塩基)単位のアルファベットを用い、NTデータベース規模の塩基配列データベース検索においてGPUによる高速化をサポートしている。検索機能に加え、Riboseekは後続のRNA解析に用いるための多重配列アラインメントや共分散モデルの生成も行う(つまり、簡単に言うと、Riboseekは、RNAやDNAの大規模な配列データベースから検索対象の配列と相同な配列を高速に見つけ、多重配列アラインメント(MSA)の作成や、RNAの構造情報を考慮した再検索まで行えるツール)
  • テストデータに含まれる8,791配列を対象にCM検索を実行したところ、メモリ不足で強制終了した。 (テスト環境:Win11、64 GBメモリ)。CM検索は通常の配列検索よりも多くのメモリと計算時間を必要とするため、大規模なデータベースを検索する場合は、十分な計算資源が必要。
  • 普通の検索(easy-searchやsearch)では、DNA配列もFASTA形式で使用できる。ただしRiboseekはRNA向けに設計されたヌクレオチド配列検索ツールのため、CM検索では、内部的にTをUとして扱い、RNAの二次構造を考慮した検索となる。
  • easy-searchは、FASTA入力から通常の相同性検索を行い、結果をm8形式で出力するところまで行う。現時点では、CM構築や構造を考慮した再検索までを一括実行するeasy-cmsearchのようなラッパーコマンドは用意されていない。

引用

Fast remote nucleotide sequence alignment with Riboseek

Sukhwan Park,  Kieran Didi,  Andrew Favor,  Anton Bushuiev,  Soohyun Kim,  View ORCID ProfileMilot Mirdita,  Martin Steinegger

bioRxiv, Posted July 31, 2026.

 

関連

 

参考