macでインフォマティクス

macでインフォマティクス

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

IQ-TREE 2

 

2020年の論文

IQ-TREE(http://www.iqtree.org)は、最尤法を用いた系統推論を行うための、ユーザーフレンドリーで広く利用されているソフトウェアパッケージである。2014年にバージョン1をリリースして以来、本著者らはIQ-TREEを継続的に拡張し、ゲノムデータを扱うための配列進化の新しいモデルや系統推定の効率的な計算アプローチを多数統合してきた。ここでは、IQ-TREEバージョン2の特筆すべき機能について説明し、他のソフトウェアと比較した際の主な利点を強調する。

 

2015年の論文

大規模な系統学データセットでは、特に最尤(ML)系統樹のための高速な樹推定法が必要とされる。高速なプログラムも存在するが、最適な樹形を見つけるinherent heuristicsのため、最適なツリーが見つかっているかどうかは明らかではない。従って、MLツリーを見つけるために異なる探索戦略を採用し、同時に現在利用可能なMLプログラムと同程度に高速なアプローチを追加する必要がある。本著者らは、ヒルクライムアプローチとstochastic perturbation methodの組み合わせが時間効率よく実装できることを示す。RAxMLやPhyMLと同じCPU時間であれば、IQ-TREEは62.2%から87.1%より高い尤度を発見し、効率的にツリー空間を探索した。IQ-TREEの停止ルールを用いた場合、DNAアラインメントでは75.7%、47.1%、タンパク質アラインメントでは42.2%、100%でRAxMLとPhyMLの方が高速となる。しかし、IQ-TREEにより高い尤度が得られる範囲は73.3-97.1%に改善された。

 

Manual

http://www.iqtree.org/doc/

Beginner's Tutorial

http://www.iqtree.org/doc/Tutorial

Advanced Tutorial

http://www.iqtree.org/doc/Advanced-Tutorial

Substitution models

http://www.iqtree.org/doc/Substitution-Models

コマンドリファレンス

http://www.iqtree.org/doc/Command-Reference

 

v1系とv2以降ではコマンドが微妙に異なっています。ここではv2を簡単に紹介します。

インストール

Github

https://github.com/iqtree/iqtree2

#link(iqtreeと指定すると通常v2が導入される)
conda install bioconda::iqtree

HPからもそれぞれのプラットフォームのバージョンのIQ-TREEをダウンロードできる。2024年現在v2.3.4となっている。

アクセスして、一番下にスクロールする。

> iqtree2 -h
IQ-TREE multicore version 2.3.4 COVID-edition for Linux x86 64-bit built Jun  5 2024
Developed by Bui Quang Minh, James Barbetti, Nguyen Lam Tung, Olga Chernomor,
Heiko Schmidt, Dominik Schrempf, Michael Woodhams, Ly Trong Nhan, Thomas Wong

Usage: iqtree [-s ALIGNMENT] [-p PARTITION] [-m MODEL] [-t TREE] ...

GENERAL OPTIONS:
  -h, --help           Print (more) help usages
  -s FILE[,...,FILE]   PHYLIP/FASTA/NEXUS/CLUSTAL/MSF alignment file(s)
  -s DIR               Directory of alignment files
  --seqtype STRING     BIN, DNA, AA, NT2AA, CODON, MORPH (default: auto-detect)
  -t FILE|PARS|RAND    Starting tree (default: 99 parsimony and BIONJ)
  -o TAX[,...,TAX]     Outgroup taxon (list) for writing .treefile
  --prefix STRING      Prefix for all output files (default: aln/partition)
  --seed NUM           Random seed number, normally used for debugging purpose
  --safe               Safe likelihood kernel to avoid numerical underflow
  --mem NUM[G|M|%]     Maximal RAM usage in GB | MB | %
  --runs NUM           Number of indepedent runs (default: 1)
  -v, --verbose        Verbose mode, printing more messages to screen
  -V, --version        Display version number
  --quiet              Quiet mode, suppress printing to screen (stdout)
  -fconst f1,...,fN    Add constant patterns into alignment (N=no. states)
  --epsilon NUM        Likelihood epsilon for parameter estimate (default 0.01)
  -T NUM|AUTO          No. cores/threads or AUTO-detect (default: 1)
  --threads-max NUM    Max number of threads for -T AUTO (default: all cores)

CHECKPOINT:
  --redo               Redo both ModelFinder and tree search
  --redo-tree          Restore ModelFinder and only redo tree search
  --undo               Revoke finished run, used when changing some options
  --cptime NUM         Minimum checkpoint interval (default: 60 sec and adapt)

PARTITION MODEL:
  -p FILE|DIR          NEXUS/RAxML partition file or directory with alignments
                       Edge-linked proportional partition model
  -q FILE|DIR          Like -p but edge-linked equal partition model
  -Q FILE|DIR          Like -p but edge-unlinked partition model
  -S FILE|DIR          Like -p but separate tree inference
  --subsample NUM      Randomly sub-sample partitions (negative for complement)
  --subsample-seed NUM Random number seed for --subsample

LIKELIHOOD/QUARTET MAPPING:
  --lmap NUM           Number of quartets for likelihood mapping analysis
  --lmclust FILE       NEXUS file containing clusters for likelihood mapping
  --quartetlh          Print quartet log-likelihoods to .quartetlh file

TREE SEARCH ALGORITHM:
  --ninit NUM          Number of initial parsimony trees (default: 100)
  --ntop NUM           Number of top initial trees (default: 20)
  --nbest NUM          Number of best trees retained during search (defaut: 5)
  -n NUM               Fix number of iterations to stop (default: OFF)
  --nstop NUM          Number of unsuccessful iterations to stop (default: 100)
  --perturb NUM        Perturbation strength for randomized NNI (default: 0.5)
  --radius NUM         Radius for parsimony SPR search (default: 6)
  --allnni             Perform more thorough NNI search (default: OFF)
  -g FILE              (Multifurcating) topological constraint tree file
  --fast               Fast search to resemble FastTree
  --polytomy           Collapse near-zero branches into polytomy
  --tree-fix           Fix -t tree (no tree search performed)
  --treels             Write locally optimal trees into .treels file
  --show-lh            Compute tree likelihood without optimisation
  --terrace            Check if the tree lies on a phylogenetic terrace

ULTRAFAST BOOTSTRAP/JACKKNIFE:
  -B, --ufboot NUM     Replicates for ultrafast bootstrap (>=1000)
  -J, --ufjack NUM     Replicates for ultrafast jackknife (>=1000)
  --jack-prop NUM      Subsampling proportion for jackknife (default: 0.5)
  --sampling STRING    GENE|GENESITE resampling for partitions (default: SITE)
  --boot-trees         Write bootstrap trees to .ufboot file (default: none)
  --wbtl               Like --boot-trees but also writing branch lengths
  --nmax NUM           Maximum number of iterations (default: 1000)
  --nstep NUM          Iterations for UFBoot stopping rule (default: 100)
  --bcor NUM           Minimum correlation coefficient (default: 0.99)
  --beps NUM           RELL epsilon to break tie (default: 0.5)
  --bnni               Optimize UFBoot trees by NNI on bootstrap alignment

NON-PARAMETRIC BOOTSTRAP/JACKKNIFE:
  -b, --boot NUM       Replicates for bootstrap + ML tree + consensus tree
  -j, --jack NUM       Replicates for jackknife + ML tree + consensus tree
  --jack-prop NUM      Subsampling proportion for jackknife (default: 0.5)
  --bcon NUM           Replicates for bootstrap + consensus tree
  --bonly NUM          Replicates for bootstrap only
  --tbe                Transfer bootstrap expectation

SINGLE BRANCH TEST:
  --alrt NUM           Replicates for SH approximate likelihood ratio test
  --alrt 0             Parametric aLRT test (Anisimova and Gascuel 2006)
  --abayes             approximate Bayes test (Anisimova et al. 2011)
  --lbp NUM            Replicates for fast local bootstrap probabilities

MODEL-FINDER:
  -m TESTONLY          Standard model selection (like jModelTest, ProtTest)
  -m TEST              Standard model selection followed by tree inference
  -m MF                Extended model selection with FreeRate heterogeneity
  -m MFP               Extended model selection followed by tree inference
  -m ...+LM            Additionally test Lie Markov models
  -m ...+LMRY          Additionally test Lie Markov models with RY symmetry
  -m ...+LMWS          Additionally test Lie Markov models with WS symmetry
  -m ...+LMMK          Additionally test Lie Markov models with MK symmetry
  -m ...+LMSS          Additionally test strand-symmetric models
  --mset STRING        Restrict search to models supported by other programs
                       (raxml, phyml, mrbayes, beast1 or beast2)
  --mset STR,...       Comma-separated model list (e.g. -mset WAG,LG,JTT)
  --msub STRING        Amino-acid model source
                       (nuclear, mitochondrial, chloroplast or viral)
  --mfreq STR,...      List of state frequencies
  --mrate STR,...      List of rate heterogeneity among sites
                       (e.g. -mrate E,I,G,I+G,R is used for -m MF)
  --cmin NUM           Min categories for FreeRate model [+R] (default: 2)
  --cmax NUM           Max categories for FreeRate model [+R] (default: 10)
  --merit AIC|AICc|BIC  Akaike|Bayesian information criterion (default: BIC)
  --mtree              Perform full tree search for every model
  --madd STR,...       List of mixture models to consider
  --mdef FILE          Model definition NEXUS file (see Manual)
  --modelomatic        Find best codon/protein/DNA models (Whelan et al. 2015)

PARTITION-FINDER:
  --merge              Merge partitions to increase model fit
  --merge greedy|rcluster|rclusterf
                       Set merging algorithm (default: rclusterf)
  --merge-model 1|all  Use only 1 or all models for merging (default: 1)
  --merge-model STR,...
                       Comma-separated model list for merging
  --merge-rate 1|all   Use only 1 or all rate heterogeneity (default: 1)
  --merge-rate STR,...
                       Comma-separated rate list for merging
  --rcluster NUM       Percentage of partition pairs for rcluster algorithm
  --rclusterf NUM      Percentage of partition pairs for rclusterf algorithm
  --rcluster-max NUM   Max number of partition pairs (default: 10*partitions)

SUBSTITUTION MODEL:
  -m STRING            Model name string (e.g. GTR+F+I+G)
                 DNA:  HKY (default), JC, F81, K2P, K3P, K81uf, TN/TrN, TNef,
                       TIM, TIMef, TVM, TVMef, SYM, GTR, or 6-digit model
                       specification (e.g., 010010 = HKY)
             Protein:  LG (default), Poisson, cpREV, mtREV, Dayhoff, mtMAM,
                       JTT, WAG, mtART, mtZOA, VT, rtREV, DCMut, PMB, HIVb,
                       HIVw, JTTDCMut, FLU, Blosum62, GTR20, mtMet, mtVer, mtInv, FLAVI,
                        Q.LG, Q.pfam, Q.pfam_gb, Q.bird, Q.mammal, Q.insect, Q.plant, Q.yeast
     Protein mixture:  C10,...,C60, EX2, EX3, EHO, UL2, UL3, EX_EHO, LG4M, LG4X
              Binary:  JC2 (default), GTR2
     Empirical codon:  KOSI07, SCHN05
   Mechanistic codon:  GY (default), MG, MGK, GY0K, GY1KTS, GY1KTV, GY2K,
                       MG1KTS, MG1KTV, MG2K
Semi-empirical codon:  XX_YY where XX is empirical and YY is mechanistic model
      Morphology/SNP:  MK (default), ORDERED, GTR
      Lie Markov DNA:  1.1, 2.2b, 3.3a, 3.3b, 3.3c, 3.4, 4.4a, 4.4b, 4.5a,
                       4.5b, 5.6a, 5.6b, 5.7a, 5.7b, 5.7c, 5.11a, 5.11b, 5.11c,
                       5.16, 6.6, 6.7a, 6.7b, 6.8a, 6.8b, 6.17a, 6.17b, 8.8,
                       8.10a, 8.10b, 8.16, 8.17, 8.18, 9.20a, 9.20b, 10.12,
                       10.34, 12.12 (optionally prefixed by RY, WS or MK)
      Non-reversible:  STRSYM (strand symmetric model, equiv. WS6.6),
                       NONREV, UNREST (unrestricted model, equiv. 12.12)
                       NQ.pfam, NQ.bird, NQ.mammal, NQ.insect, NQ.plant, NQ.yeast
           Otherwise:  Name of file containing user-model parameters

STATE FREQUENCY:
  -m ...+F             Empirically counted frequencies from alignment
  -m ...+FO            Optimized frequencies by maximum-likelihood
  -m ...+FQ            Equal frequencies
  -m ...+FRY           For DNA, freq(A+G)=1/2=freq(C+T)
  -m ...+FWS           For DNA, freq(A+T)=1/2=freq(C+G)
  -m ...+FMK           For DNA, freq(A+C)=1/2=freq(G+T)
  -m ...+Fabcd         4-digit constraint on ACGT frequency
                       (e.g. +F1221 means f_A=f_T, f_C=f_G)
  -m ...+FU            Amino-acid frequencies given protein matrix
  -m ...+F1x4          Equal NT frequencies over three codon positions
  -m ...+F3x4          Unequal NT frequencies over three codon positions

RATE HETEROGENEITY AMONG SITES:
  -m ...+I             A proportion of invariable sites
  -m ...+G[n]          Discrete Gamma model with n categories (default n=4)
  -m ...*G[n]          Discrete Gamma model with unlinked model parameters
  -m ...+I+G[n]        Invariable sites plus Gamma model with n categories
  -m ...+R[n]          FreeRate model with n categories (default n=4)
  -m ...*R[n]          FreeRate model with unlinked model parameters
  -m ...+I+R[n]        Invariable sites plus FreeRate model with n categories
  -m ...+Hn            Heterotachy model with n classes
  -m ...*Hn            Heterotachy model with n classes and unlinked parameters
  --alpha-min NUM      Min Gamma shape parameter for site rates (default: 0.02)
  --gamma-median       Median approximation for +G site rates (default: mean)
  --rate               Write empirical Bayesian site rates to .rate file
  --mlrate             Write maximum likelihood site rates to .mlrate file

POLYMORPHISM AWARE MODELS (PoMo):
  -s FILE              Input counts file (see manual)
  -m ...+P             DNA substitution model (see above) used with PoMo
  -m ...+N<POPSIZE>    Virtual population size (default: 9)
  -m ...+WB|WH|S]      Weighted binomial sampling
  -m ...+WH            Weighted hypergeometric sampling
  -m ...+S             Sampled sampling
  -m ...+G[n]          Discrete Gamma rate with n categories (default n=4)

COMPLEX MODELS:
  -m "MIX{m1,...,mK}"  Mixture model with K components
  -m "FMIX{f1,...fK}"  Frequency mixture model with K components
  --mix-opt            Optimize mixture weights (default: detect)
  -m ...+ASC           Ascertainment bias correction
  --tree-freq FILE     Input tree to infer site frequency model
  --site-freq FILE     Input site frequency model file
  --freq-max           Posterior maximum instead of mean approximation

TREE TOPOLOGY TEST:
  --trees FILE         Set of trees to evaluate log-likelihoods
  --test NUM           Replicates for topology test
  --test-weight        Perform weighted KH and SH tests
  --test-au            Approximately unbiased (AU) test (Shimodaira 2002)
  --sitelh             Write site log-likelihoods to .sitelh file

ANCESTRAL STATE RECONSTRUCTION:
  --ancestral          Ancestral state reconstruction by empirical Bayes
  --asr-min NUM        Min probability of ancestral state (default: equil freq)

TEST OF SYMMETRY:
  --symtest               Perform three tests of symmetry
  --symtest-only          Do --symtest then exist
  --symtest-remove-bad    Do --symtest and remove bad partitions
  --symtest-remove-good   Do --symtest and remove good partitions
  --symtest-type MAR|INT  Use MARginal/INTernal test when removing partitions
  --symtest-pval NUMER    P-value cutoff (default: 0.05)
  --symtest-keep-zero     Keep NAs in the tests

CONCORDANCE FACTOR ANALYSIS:
  -t FILE              Reference tree to assign concordance factor
  --gcf FILE           Set of source trees for gene concordance factor (gCF)
  --df-tree            Write discordant trees associated with gDF1
  --scf NUM            Number of quartets for site concordance factor (sCF)
  --scfl NUM           Like --scf but using likelihood (recommended)
  -s FILE              Sequence alignment for --scf
  -p FILE|DIR          Partition file or directory for --scf
  --cf-verbose         Write CF per tree/locus to cf.stat_tree/_loci
  --cf-quartet         Write sCF for all resampled quartets to .cf.quartet

ALISIM: ALIGNMENT SIMULATOR

Usage: iqtree --alisim <OUTPUT_PREFIX> [-m MODEL] [-t TREE] ...

  --alisim OUTPUT_ALIGNMENT Activate AliSim and specify the output alignment filename
  -t TREE_FILE              Set the input tree file name
  --length LENGTH           Set the length of the root sequence
  --num-alignments NUMBER   Set the number of output datasets
  --seqtype STRING          BIN, DNA, AA, CODON, MORPH{NUM_STATES} (default: auto-detect)
                            For morphological data, 0<NUM_STATES<=32
  --m MODEL_STRING          Specify the evolutionary model. See Manual for more detail
  --mdef FILE               Name of a NEXUS model file to define new models (see Manual)
  --fundi TAXA_LIST,RHO     Specify a list of taxa, and Rho (Fundi weight) for FunDi model
  --indel <INS>,<DEL>       Set the insertion and deletion rate of the indel model,
                            relative to the substitution rate
  --indel-size <INS_DIS>,<DEL_DIS> Set the insertion and deletion size distributions
  --sub-level-mixture       Enable the feature to simulate substitution-level mixture model
  --no-unaligned            Disable outputing a file of unaligned sequences
                            when using indel models
  --root-seq FILE,SEQ_NAME  Specify the root sequence from an alignment
  -s FILE                   Specify the input sequence alignment
  --no-copy-gaps            Disable copying gaps from input alignment (default: false)
  --site-freq <OPTION>      Specify the option (MEAN (default), or SAMPLING, or MODEL)
                            to mimic the site-frequencies for mixture models from
                            the input alignment (see Manual)
  --site-rate <OPTION>      Specify the option (MEAN (default), or SAMPLING, or MODEL)
                            to mimic the discrete rate heterogeneity from
                            the input alignment (see Manual)
  -t RANDOM{MODEL,NUM_TAXA} Specify the model and the number of taxa to generate a random tree
  -rlen MIN MEAN MAX        Specify three numbers: minimum, mean and maximum branch lengths
                            when generating a random tree
  -p FILE                   NEXUS/RAxML partition file
                            Edge-linked proportional partition model
  -q FILE                   Like -p but edge-linked equal partition model
  -Q FILE                   Like -p but edge-unlinked partition model
  --distribution FILE       Supply a definition file of distributions,
                            which could be used to generate random model parameters
  --branch-distribution DIS Specify a distribution, from which branch lengths of the input trees
                            are randomly generated and overridden.
  --branch-scale SCALE      Specify a value to scale all branch lengths
  --single-output           Output all alignments into a single file
  --write-all               Enable outputting internal sequences
  --seed NUM                Random seed number (default: CPU clock)
                            Be careful to make the AliSim reproducible,
                            users should specify the seed number
  -gz                       Enable output compression but taking longer running time
  -af phy|fasta             Set the output format (default: phylip)
  User Manual is available at http://www.iqtree.org/doc/alisim

ANALYSIS WITH GENTRIUS ALGORITHM:
  --gentrius FILE      File must contain either a single species-tree or a set of subtrees.
  -pr_ab_matrix FILE   Presence-absence matrix of loci coverage.
  -s FILE              PHYLIP/FASTA/NEXUS/CLUSTAL/MSF alignment file(s)
  -p FILE              NEXUS/RAxML partition file
  -g_stop_t NUM        Stop after NUM species-trees were generated, or use 0 to turn off this stopping rule. Default: 1MLN trees.
  -g_stop_i NUM        Stop after NUM intermediate trees were visited, or use 0 to turn off this stopping rule. Default: 10MLN trees.
  -g_stop_h NUM        Stop after NUM hours (CPU time), or use 0 to turn off this stopping rule. Default: 7 days.
  -g_non_stop          Turn off all stopping rules.
  -g_query FILE        Species-trees to test for identical set of subtrees.
  -g_print             Write all generated species-trees. WARNING: there might be millions of trees!
  -g_print_lim NUM     Limit on the number of species-trees to be written.
  -g_print_induced     Write induced partition subtrees.
  -g_print_m           Write presence-absence matrix.
  -g_rm_leaves NUM     Invoke reverse analysis for complex datasets.

TIME TREE RECONSTRUCTION:
  --date FILE          File containing dates of tips or ancestral nodes
  --date TAXNAME       Extract dates from taxon names after last '|'
  --date-tip STRING    Tip dates as a real number or YYYY-MM-DD
  --date-root STRING   Root date as a real number or YYYY-MM-DD
  --date-ci NUM        Number of replicates to compute confidence interval
  --clock-sd NUM       Std-dev for lognormal relaxed clock (default: 0.2)
  --date-no-outgroup   Exclude outgroup from time tree
  --date-outlier NUM   Z-score cutoff to remove outlier tips/nodes (e.g. 3)
  --date-options ".."  Extra options passing directly to LSD2
  --dating STRING      Dating method: LSD for least square dating (default)

 

 

実行方法

実行するには、多重整列されたPHYLIPフォーマットかFASTAフォーマット、あるいはNEXUS/CLUSTAL/MSFの配列を準備する必要がある。

IQ-treeチュートリアルより

 

IQ-TREEを実行するには、多重整列された配列セットを指定する。デフォルト設定の最尤法の系統推定が行われる。

#PHYLIP
iqtree -s example.phy

#FASTA
iqtree -s example.fa
  • -s     PHYLIP/FASTA/NEXUS/CLUSTAL/MSF alignment file(s)

 

 

#PHYLIP
iqtree -s example.phy

#FASTA
iqtree -s example.fa

 

全プロセスを再実行する。

iqtree -s example.aln --redo
  • --redo           Redo both ModelFinder and tree search
  • --redo-tree   Restore ModelFinder and only redo tree search
  • --undo          Revoke finished run, used when changing some options 

 

 

IQ-TREEは自動で樹形推定プロセスを進めるため、ログを読むことが重要(.logとしてファイル保存されている)

勉強会用の資料を貼っておきます。

 

引用

IQ-TREE 2: New Models and Efficient Methods for Phylogenetic Inference in the Genomic Era 
Bui Quang Minh, Heiko A Schmidt, Olga Chernomor, Dominik Schrempf, Michael D Woodhams, Arndt von Haeseler, Robert Lanfear
Molecular Biology and Evolution, Volume 37, Issue 5, May 2020, Pages 1530–1534

 

IQ-TREE: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies
Lam-Tung Nguyen, Heiko A. Schmidt, Arndt von Haeseler, and Bui Quang Minh

Molecular Biology and Evolution. 2015 Jan; 32(1): 268–274.

 

ModelFinder: fast model selection for accurate phylogenetic estimates

Subha Kalyaanamoorthy, Bui Quang Minh, Thomas K F Wong, Arndt von Haeseler & Lars S Jermiin 
Nature Methods volume 14, pages587–589 (2017)

 

Ultrafast Approximation for Phylogenetic Bootstrap 
Bui Quang Minh, Minh Anh Thi Nguyen, Arndt von Haeseler
Molecular Biology and Evolution, Volume 30, Issue 5, May 2013, Pages 1188–1195

 

参考

http://www.tezuru-mozuru.com/?cat=200

 

NCBI BLASTのClusteredNR database

 

近年、配列決定技術の高度化によってNCBIのタンパク質NRデータベースは急速に成長しており、特定の種の生物のタンパク質は過剰に公開されている。このため、特に過剰に読まれた生物の配列かそれに進化的に近縁な生物の配列を使ってBLASTサーチを実行すると、非常に近縁な生物種から得られた同じ種類のタンパク質が結果の大部分を占めることがある。しかし、BLASTサーチでは、必ずしもクエリの配列に近い情報の重要性が高いわけではない。クエリの配列とは進化的に遠かったり、少しだけ似ている別の機能を持った配列の情報を取得することが重要な事も多い。BLAST実行前にヒット数を増やしてジョブを投げればより遠縁なヒットも得られるが、NCBI BLASTはリソース制限を行っており、制限を超えたジョブは途中で止まってしまう。また、ヒット数が増える事は抜本的な解決にはなっておらず、却って視認性を悪くしてしまう。

NCBI ClusteredNRデータベースは、2022年に登場した、MMseqs2を用いて類似した配列をクラスタリングして得られた、冗長性を減らしたNRデータベースとなっている。各クラスタには互いに90%以上同一で、最長の配列長の90%以内の長さのタンパク質が含まれている。 2024年現在でもExperimental だが、データベースをクラスタリングすることで、元のデータベースに含まれる生物やタンパク質の多様性をよりよく表現したデータベースとなっていて、検索にかかる時間も短くなってるなど、状況によってはデフォルトのNRデータベースより扱いやすい。通常のBLASTとは結果の見方が異なるので、使い方を簡単に確認しておきます。

 

NCBI insights

https://ncbiinsights.ncbi.nlm.nih.gov/2022/05/02/clusterednr_1/

 

webサービス

NCBI BLASTPかBLASTXにアクセスする。

 

Choose search set => Experimental databaseでClusteredNRを選択する。

Organisumの欄で指定することで、BLASTサーチの対象範囲を特定の分類のみに制限できる。制限するには、生物の一般名、属名+種小名の二命名法、またはNCBI taxIDで指定する。

 

ここではE.coliの60S ribosomal protein L21(link)のアミノ酸配列を使用した。パラメータはデフォルトとした。

 

出力例

通常のBLAST結果と異なり、結果はクラスターごとにまとめられる。

 

クラスタの情報として、アノテーションの充実した代表タンパク質の情報が表示されている。

クラスタは90%以上同一で、最長の配列長の90%以内の長さでクラスタリングしているので、複数の生物(種)の配列が含まれる場合がある。

 

右端のPer identityの列を見ると、トップヒットが100%、次が91%となっていて冗長なヒットが抑制されていることが分かる。

 

例えば上の画像の3つ目のヒットには 8 members, 58 organismとある。クリックすると、そのクラスタに含まれる配列の情報が表示される。

 

このパネルの右上のボタンから、全配列をダウンロードできる。

また、クラスターの全メンバーのMSAを実行したり、含まれる配列に対してBLASTのジョブを投げたりもできる。

 

8配列含まれるクラスタへのBLAST結果(BLAST alignmentボタン)

 

多重整列結果

 

クラスタは事前定義されており、探索は高速に実行できます。クエリと類似した配列を集めてきたい時にも便利だと思います。

引用

Database resources of the national center for biotechnology information
Eric W Sayers, Evan E Bolton, J Rodney Brister, Kathi Canese, Jessica Chan, Donald C Comeau, Ryan Connor, Kathryn Funk, Chris Kelly, Sunghwan Kim, Tom Madej, Aron Marchler-Bauer, Christopher Lanczycki, Stacy Lathrop, Zhiyong Lu, Francoise Thibaud-Nissen, Terence Murphy, Lon Phan, Yuri Skripchenko, Tony Tse, Jiyao Wang, Rebecca Williams, Barton W Trawick, Kim D Pruitt, Stephen T Sherry

Nucleic Acids Res. 2022 Jan 7;50(D1):D20-D26. doi:0.1093/nar/gkab1112.

 

 

ゲノム領域を柔軟に調整して視覚化と比較ができるユーザーフレンドリーなアプリケーション GenoFig

 

 生物の分子進化の歴史を理解するには、通常、近縁種や系統のゲノム領域を視覚的に比較する必要がある。このタスクを達成するためのアプリケーションはすでにいくつか存在するが、それらは古すぎたり、限定的すぎたり、あるいは複雑すぎたりして、ほとんどのユーザーのニーズには合わない。GenoFigは、原核生物のゲノム領域を視覚化するためのグラフィカルなアプリケーションであり、可能な限り使いやすく、様々なニーズに適応できる柔軟性を意図している。GenoFigは、正規表現を用いて、GenBankファイルから抽出されたアノテーションを、配列間で一貫した方法でパーソナライズされた形で表現することができる。また、配列間の相同領域の表示を最適化するユニークなオプションや、配列のGCパーセンテージやGC-skew表現のような、より古典的な機能も提供する。要約すると、GenoFigは、原核生物における特定のゲノム領域の進化を探索し、出版に耐えうる図を作成するための、シンプルで、無料で、高度に設定可能なツールである。Genofigは、GPL 3.0ライセンスのもと、https://forgemia.inra.fr/public-pgba/genofigで利用できる。

 

インストール

ubuntu22でcondaで環境を作ってテストした。レポジトリではリリース1.1のWindows版とMacOS版のコンパイル済みバージョンのダウンロードリンクも用意されている。

Github

git clone https://forgemia.inra.fr/public-pgba/genofig.git
cd genofig/
#linux
mamba env create -f extras/requirements.yml
conda activate genofig
extras/SETUP.sh #link
#windows (WSLではない)
mamba env create -f extras/requirements_windows.yml


#パスを通す
export PATH="'$(pwd)':$PATH"
#or,
echo 'export PATH="'$(pwd)':$PATH"' >> ~/.bashrc && source ~/.bashrc

> Genofig


GenoFigはGenBankフォーマットのアノテーションを使う。GenoFigは小さなゲノム領域(せいぜい数百Kbp)を比較するように設計されているので、ゲノム全体を読み込むには向いていない。数個の遺伝子や遺伝子クラスターなどを可視化する。

 

チュートリアルではNCBI nucleotideから遺伝子フィーチャーを可視化して、比較したい領域をGenBankでダウンロードしている。

GenBankは左端のプラスのボタンから読み込む。

 

CREATE FIGURE ボタンを押すと視覚化されてSVG形式で保存される。

アノテーションがついているCDSは灰色のボックスでプロットされる。例外として、アノテーションが hypothetical proteinsは白色のボックスでプロットされる。

 

特定のCDSだけ色を変更できる。Featureタブに移動し、プラスマークで新しいfeatureを追加、type列で色を変更したいCDSアノテーション名(同じ名前だと一括して適用される)を入力し、productを選択、さらに右の方のcolorで色を指定する。

視覚化すると色が変更された。フィルターは正規表現をサポートしているので、柔軟に特定のfeatureだけ強調したりできる。

 

レポジトリで紹介されている代表的な機能を見ていく。

 

配列タブでreverseにチェックを入れると全ての矢印の向きが反転する(全体が反転するのではなく、個々のfeatureが180度回転する)。

SequencesタブのBとIをチェックして左端のラベルを太字と斜体にする。

 

Sequencesタブの左端のラベルの色を緑にして、さらにtypeをorganismにして生物名がラベルされるようにする。

 

Sequencesタブの左端のラベルのサイズを20に下げる。

Sequencesタブには、他にもラベルの位置やサイズ、幅などを細かく調整できるようになっている。

 

ラベルがはみ出さないように、上のパネルで左側のマージンを400に増やす。

 

2つ目の配列を追加した。

 

ラベルが被っている。ほとんどのフィーチャーは同じなので、Sequencesタブの右端近くのfeat.label列のGeneralのチェックをはずし、最初の配列にだけチェックを入れる。

一番上のgenbank配列にのみアノテーションが表示された。

 

Featuresタブの右端の方のlabel列のチェックを特定のフィーチャー以外は外した。

 

Homologies パネルでRun blastnを実行し、それから描画した。

配列間の相同なCDSシンテニーで示される(すべての配列間でBLASTnが実行されるため配列数が多い時は注意)。画像では、上下は同じ配列のため、全領域にシンテニーが見つかっている。

 

レジェンドタブで凡例が表示されるようにする。サイズはスケール係数で、場所はpositonで指定できる。

カスタムした設定はsaveから保存できるが、 そのApplicationメニューが正しく表示されないという問題が発生しているらしい。試した時は表示されなかった(M1 macでrosetta2使用)。

 

レポジトリでは情報が整理されてずっと綺麗な作図がされています。確認して下さい。

引用

GenoFig: a user-friendly application for the visualisation and comparison of genomic regions 

Maxime Branger,   Sébastien O Leclercq

Bioinformatics, Published: 13 June 2024

 

関連

遺伝子クラスターを比較してインタラクティブな図で視覚化する clinker

細菌・古細菌の環状ゲノムプロットを出力する GenoVi

 

タンパク質の機能的アノテーションを行う AnnoPRO

 

 タンパク質の機能アノテーションは生物科学における長年の課題の一つであり、様々な計算手法が開発されてきた。しかし、既存の方法では、GOファミリーの数が多く、アノテーションされたタンパク質が少ないという深刻なロングテール問題に悩まされている。そこで、配列に基づくマルチスケールタンパク質表現、事前学習によるデュアルパスタンパク質エンコーディング、長期短期記憶に基づくデコーディングによる機能アノテーションを可能にする、AnnoPROと名付けられた革新的な戦略を構築した。様々なベンチマークに基づくケーススタディを実施し、AnnoPROが利用可能な手法の中で優れた性能を持つことを確認した。ソースコードとモデルは https://github.com/idrblab/AnnoPROhttps://zenodo.org/records/10012272にあり、自由に利用できる。

 

手順

step 1: 入力タンパク質のシーケンス
step 2: Profeatによる特徴抽出
step 3: 特徴の対距離計算 --> cosine, correlation, jaccard
Step4: 特徴の2次元埋め込み --> umap, tsne, mds
step5: 特徴の格子配置 --> grid, scatter
Step6: 変換 --> minmax, standard

インストール

ubuntu22でcondaで環境を作ってテストした(RTX3090使用)。

Github

git clone https://github.com/idrblab/AnnoPRO.git
cd AnnoPRO
conda create -n annopro python=3.8
conda activate annopro
pip install .

> annopro -h

usage: annopro [-h] [--fasta_file FASTA_FILE] [--output OUTPUT] [--used_gpu USED_GPU] [--disable_diamond] [--overwrite] [--version]

 

Arguments for AnnoPRO

 

optional arguments:

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

  --fasta_file FASTA_FILE, -i FASTA_FILE

                        The protein sequences file

  --output OUTPUT, -o OUTPUT

                        Output directory

  --used_gpu USED_GPU   GPU device selected, default is CPU

  --disable_diamond     Disable blast with diamond

  --overwrite           Overwrite existed output

  --version             Show version

 

 

実行方法

タンパク質のfastaファイルを指定する。配列名に特殊文字が含まれているとエラーを起こす可能性があるので注意する。

annopro -i test_proteins.fasta -o output
  • --used_gpu     GPU device selected, default is CPU

初回はモデルダウンロードに時間がかかる。試した時は、1度目は途中で止まり、2度目は1時間ほどでダウンロードできた。

 

出力例

result.csvが最終出力となる。Gene OntologyのBP、MF、CCのカテゴリに分けて3つのファイルとして出力されている。

 

> csvlook cc_result.csv |head -n 20

1列目:タンパク質名、2列目:GO term、3列目:スコア

 

レポジトリのexample下にも出力ファイル例がある。

https://github.com/idrblab/AnnoPRO/tree/main/example

引用

AnnoPRO: a strategy for protein function annotation based on multi-scale protein representation and a hybrid deep learning of dual-path encoding

Lingyan Zheng, Shuiyang Shi, Mingkun Lu, Pan Fang, Ziqi Pan, Hongning Zhang, Zhimeng Zhou, Hanyu Zhang, Minjie Mou, Shijie Huang, Lin Tao, Weiqi Xia, Honglin Li, Zhenyu Zeng, Shun Zhang, Yuzong Chen, Zhaorong Li & Feng Zhu 

Genome Biology volume 25, Article number: 41 (2024) 

 

関連

https://kazumaxneo.hatenablog.com/entry/2024/04/15/004855

 

生物間の遺伝子機能の類似点と相違点をインタラクティブに探索するウェブツール Comparative Genome Dashboard

 

 Comparative Genome Dashboardは、生物間の遺伝子機能の類似点と相違点をインタラクティブに探索するためのウェブベースのソフトウェアツールである。このツールは細胞機能のハイレベルなグラフィカルな調査を提供し、興味のあるサブシステムをより詳細に調べるためにドリルダウン(注;データのレベルを掘り下げてさらに詳細を調べる操作)することを可能にする。比較ダッシュボードの最も高いレベルには、生合成、エネルギー代謝、輸送、刺激に対する応答などの細胞システムのパネルが含まれている。各パネルには、そのパネルのサブシステムのセットについての各生物の化合物または遺伝子産物の数をプロットした棒グラフのセットが含まれている。ユーザーは、興味のあるサブシステムにフォーカスするためにインタラクティブにドリルダウンし、各生物によって生成または消費される化合物のグリッド、特定のGO termの割り当て、パスウェイダイアグラム、およびより詳細な比較ページへのリンクを見ることができる。例えば、ダッシュボードでは、一組の生物が合成できる補酵素、輸送できる金属イオン、DNA損傷修復能力、バイオフィルム形成遺伝子、ウイルス応答タンパク質を比較することができる。このダッシュボードにより、ユーザーは様々な詳細レベルでの包括的な比較を素早く行うことができる。

 

BioCyc Guided Tour (2023)

https://youtu.be/t3xuKH7_Txo?si=Uw4jkzwOzit1GLNt

 

HPより

Comparative Genome Dashboardは、それぞれのゲノムとパスウェイのアノテーションによってコードされる、生物または生物群の全体的な生物学的能力を視覚化するためのツールである。全ての細胞システムの迅速な調査を容易にし、生物間の類似点と相違点を素早く特定することを可能にする。

webサービス

BioCyc web portal: https://biocyc.orgにアクセスする。

右上からアカウントを作ってログインする。すぐに作れる。

 

アクセスするには、上のToolsメニューからComparative Genome Dashboard(右端の上から2つ目)を選択する。 

 

Comparative Genome Dashboardに入ったら、比較する生物を選ぶウィンドウが出てくる。

このウィンドウで生物を選択する。数が多いと表示が遅くなるため、10以下に留めることが推奨されている。

画像では細菌を選んでいるが、真核生物も選べる(ヒトなど)。

 

ダッシュボード・データの生成には1~2分かかる。

 

出力例

 

(マニュアルより)ゲノムダッシュボードは、分解/利用/同化といった細胞システムを表すパネルで構成されている。その構成は、オミックスダッシュボード(オミックスデータを可視化するツール:https://biocyc.org/dashboard/dashboard-intro.shtml?orgid=ECOLI)と同じである。各代謝パネルには独自のY軸があり、各生物によって分解または利用された代謝物の総数を示している。異なる色のプロットは生物に対応している(画像では5つの生物)。

代謝物のカウントは、その生物に存在する代謝経路の有無に基づいている(ゲノムアノテーションの質に依存)。

 

代謝物ではないパネルでは、棒グラフはそのシステムにアノテーションされた遺伝子数を表している。

 

プロットをクリックすると、そのサブシステムを構成するサブシステムについての新しいパネルが表示される。Cofactor Synをクリックした。

出てきたサブパネルをさらにクリックすると、各代謝物の代謝の有無を示したタイル状のプロットが表示される。

 

行は代謝物(文脈によって合成、分解、輸送)またはその他の生物学的能力を表し、列は生物を表す。

この例のように、検索後の初期パネルでは代謝物の数だけ表示されていたが、クリックしてドリルダウンする事で、どの代謝物の合成や分解が可能であるのかを詳しく調べる事ができるように設計されている。

(注;上の図のボックスは、その生物のデータベースに該当する合成パスウェイが1つ以上含まれている場合のみ色が付けられている(複数のバリアント経路が定義されている場合があることに注意)。色の薄いボックスは、その生物に経路が存在すると予測されたとしても1つ以上の経路の穴があることを示している。色が薄いほど欠けている経路の割合が大きい。)

 

棒グラフやカラーボックスにマウスカーソルを合わせると、化合物名とその化合物および生物に関連するすべてのパスウェイがリストアップされたウィンドウが表示される。

 

カラーボックスをクリックすると、その生物での該当する化合物の代謝に関連するすべてのバリアント経路のパスウェイ図表示される。

パスウェイ図の化合物をクリックすると詳細を確認できる。

 

また、サブパネルのカラーパネルの上でカーソルをホバーするとshow pathway comparisonが出てくる。これをクリックすると、

新しいウィンドウで、生物名(行)と遺伝子の象限として(それぞれの象限には構造式か遺伝子名)、その代謝反応の有無を比較できる。

 

ここからBIOCYCにジャンプしてどのような反応であるかを詳しく調べることが出来る。原核生物でオペロンになっているならオペロン構造なども調べられる。

 

検索後のトップページ(ゲノムダッシュボード)上部にある2つの検索ボタン:Search Compound/Pathway/GO-Term ボタンをクリックすると、化合物やパスウェイ、GO termの名前を入力して検索し、該当するパスウェイが比較している生物のどのカテゴリーに存在するのか調べられる。ヒットした場合、ここから該当するパスウェイに素早くアクセスできる。



 

その他・マニュアルより

  • 生物の生物学的性質を反映するだけでなく、ダッシュボードに表示される機能はゲノムアノテーションの質に依存する。
  •  化合物名の上、または色の付いていないボックスの上にマウスを置くと、ツールチップに詳細なパスウェイ比較ページへのリンクが表示され、違いが実際の生物学的根拠を持つ可能性が高いのか、またはアノテーションやパスウェイ予測の質の違いによる人工物なのかを評価するのに役立つ。
  • オプション・メニューには、その他のコマンドも用意されている。これには、パネルの内容をSVGまたはPNG画像ファイルにエクスポートする機能、比較表ページまたはダウンロード可能なデータテーブルを生成する機能などが含まれている。トップレベルのパネルについては、指定されたパネルを非表示にするコマンドもある(ページ下部のボタンで非表示のパネルを復元できる)。メイン表示ページから、パネルのタイトルをマウスでつかみ、希望の位置までドラッグすることで、トップレベルパネルの並び替えもできる。

  • 通化合物およびユニーク化合物の数を表示する、オプションを有効にすると、サブシステムのプロットを横切る黒いバーは、すべての生物(隠れた生物を含む)に共通する化合物の数を示す。分析が2つ以上の生物を含む場合、白いバーは、少なくとも1つの他の生物と共有する化合物の数(バーの下)と、その生物に固有の化合物の数(バーの上)の両方を示す。

  • Organism Preferencesセクションでは、ダッシュボードでの各生物の表示方法をカスタマイズできる。生物を選択的に非表示または表示したり、ラベルや色を編集したり、並べ替えたりできる。1つまたは複数の生物を選択的に非表示にするには、対応するボックスのチェックを外します。これにより、表示から生物のみが非表示になる。

 

コメント

アクセスできない時はメンテナンス中である可能性があります。

引用

The Comparative Genome Dashboard

Suzanne Paley, Ron Caspi, Paul O’Maille, Peter D. Karp

bioRxiv, posted June 12, 2024

 

高忠実度なin silicoモデリングによるRNA-Seqシミュレータ BEERS2

 

 RNA-seqリードのシミュレーションは、バイオインフォマティクスツールの評価、比較、ベンチマーク、開発において極めて重要である。しかし、RNA-seqシミュレータの分野は過去10年間ほとんど進歩していない。このニーズに応えるため、本著者らは柔軟で高度に設定可能な設計と、ライブラリー調製とシーケンスパイプライン全体の詳細なシミュレーションを組み合わせたBEERS2を開発した。BEERS2は、カスタマイズ可能な入力またはCAMPAREEでシミュレーションされたRNAサンプルから、入力転写産物(通常、polyAテールを持つ全長のメッセンジャーRNA転写産物)を取り込む。BEERS2は、これらの転写産物のリアルなリードをFASTQ、SAMまたはBAMフォーマットで生成する。また、真の転写産物レベルの定量値も生成する。BEERS2は、柔軟で高度なコンフィギュレーションが可能な設計と、ライブラリー調製とシーケンスパイプライン全体の詳細なシミュレーションを兼ね備えており、ribosomal depletionのためのpolyA選択とRiboZeroの効果、ヘキサマープライミング配列のバイアス、ポリメラーゼ連鎖反応(PCR)増幅におけるGC-contentのバイアス、バーコードの読み取りエラー、PCR増幅中のエラーを含むように設計されている。これらの特徴を併せ持つBEERS2は、これまでで最も完全なRNA-seqのシミュレーションである。最後に、一般的なSalmon擬似アライメントアルゴリズムに対するいくつかの設定の効果を測定することで、BEERS2の利用を実証する。

 

インストール

Github

https://github.com/itmat/BEERS2

mamba create -n beers2 python=3.11 -y
conda activate beers2
pip install git+https://github.com/itmat/BEERS2

beers -h

$ beers -h

usage: beers --configfile CONFIGFILE

beers: error: the following arguments are required: --configfile

 

テストラン

curl https://s3.amazonaws.com/itmat.data/BEERS2/examples/baby_mouse_example.tar.gz | tar -xz
=> baby_mouse_example/ができる

解凍して出来たbaby_mouse_example/baby.config.yamlを使う。実際の使用時にもこのyamlファイルをテンプレートとして使うように書かれている。

> cat baby.config.yaml

###### BEERS CONFIGURATION ######

# The following file configures BEERS 2.

 

 

# This file is setup to run the example 'baby' data, which is

# meant as an example and to verify that BEERS2 work properly.

# The data comes from mouse, with data stripped to several small

# genomic regions.

 

# For the data files necesary to use this file, please run:

# wget -c https://s3.amazonaws.com/itmat.data/BEERS2/examples/baby_mouse_example.tar.gz -O - | tar -xz

# cd baby_mouse_example/

 

# The seed value initializes the random number generator

# so that runs are reproducible. Any integer works.

seed: 100

 

output:

    # Output types, set to True to generate these data.

    # The results appear in the 'results' directory

    # FastQ Files are named as S1_L8_R2.fastq for example

    # for sample 1, lane 8, read 2.

    # Moreover, an 'unidentified' file is produced for each sample

    # that contains reads that originated from that sample but whose

    # demultiplexing step failed due to errors in the barcode reads.

    output_fastq: true

 

    # Set to true for sam/bam files containing the true alignments to

    # the reference genome. Note that such alignments may no longer

    # be replicable even in theory from just the reads, for exampe if

    # most or all of the read is a PolyA tail and so cannot be aligned.

    # This file still indicates where it came from

    # Output is in the 'results' directory and named like S1_L8.sam

    # for sample 1, lane 8. Like fastq, also has the potential for

    # unidentified reads, stored in separate files.

    output_sam: true

    output_bam: false

 

    # Set to true to receive full logs of intermediate steps

    # These can be very large for a full-sized BEERS run (100s of GB)

    # since they document every intermediate molecule throughout library prep

    full_logs: false

 

global_config:

    # global configuration contains configuration values that are available to all

    # steps of the Library Prep and Sequencing Pipeline.

 

    samples:

        # Specify sample-specific information

        # Samples are given identfiers, which are strings like '1' or '2'

        '1':

            # Sample one configuration

            barcodes:

                # Barcodes are used by AdapterLigationStep

                # Each sample should have unique (i5, i7) pair of barcodes

                # These get sequenced to determine which sample the read came from

                # Typical examples of these can be founda at:

                # https://support-docs.illumina.com/SHARE/AdapterSeq/Content/SHARE/AdapterSeq/TruSeq/UDIndexes.htm

                i5: AGCGCTAG

                i7: AACCGCGG

        '2':

            # Sample 2 configuration, same as above

            barcodes:

                i5: GATATCGA

                i7: TTATAACC

 

    molecule_maker_parameters:

        # BEERS has two different methods for generating molecules

        # Either they can be provided molecule-by-molecule in the molecule file

        # Or they can be generated on demand using CAMPAREE output sample directories

        # that contain all the information to synthesize random molecules

        # These options affect this second on-the-fly molecules.

 

        # The range of polyA tails to generate. Selected uniformly within this range.

        min_polyA_tail_length: 50

        max_polyA_tail_length: 250

 

    resources:

        # Resoruces contain general-use information

 

        # Adapters are added to the molecules in the AdapterLigationStep

        # They are also read in the sequencing by synthesis step where they

        # are used to initiate transcription.

        # Each adapter flanks the corresponding sample barcode, which all flanks the original molecule

        # So the molecule ends up looking like (5' to 3'):

        #  (pre_i5_adapter) (i5) (post_i5_adapter) (molecule sequence) (pre_i7_adapter) (i7) (post_i7_adapter)

        # These ones have been obtained from:

        # https://support-docs.illumina.com/SHARE/AdapterSeq/Content/SHARE/AdapterSeq/TruSeq/UDIndexes.htm

        pre_i5_adapter: AATGATACGGCGACCACCGAGATCTACAC

        post_i5_adapter: ACACTCTTTCCCTACACGACGCTCTTCCGATCT

        pre_i7_adapter: GATCGGAAGAGCACACGTCTGAACTCCAGTCAC

        post_i7_adapter: ATCTCGTATGCCGTCTTCTGCTTG

 

        # Reference genome FASTA is used to generate output BAM files, which require chromosome lengths.

        # Use the same FASTA as the BEERS input (such as from CAMPAREE) was

        # generated from.

        reference_genome_fasta: input_data/baby_genome.mm10/baby_genome.mm10.oneline_seqs.fa

 

library_prep_pipeline:

    # The library prep pipeline takes input molecules and walks them through a library prep

    # to get them ready for sequencing.

    input:

        # Input the Library Prep step, which is also the input to BEERS in general

        # The results of this will be fed to the sequencing pipeline to generate the output

 

        # There are TWO ways to provide input to BEERS.

        # First, we can provide molecule files, such as output by CAMPAREE

        # These specify exactly the molecules to prepare and sequence

        # Specify these by a path to a directory containing the molecule packets

        # as if generated by CAMPAREE:

        #

        # /my/directory/path:

        #   - sample1:

        #       molecule_file1.txt

        #       molecule_file2.txt

        #   - sample2:

        #       molecule_file1.txt

        #       molecule_file2.txt

        directory_path: input_data/

 

        # Second, we can provide CAMPAREE sample directories

        # Which instead contain the information to generate new molecules

        # based off a sample's variants, gene, isoform, and allele distributions

        from_distribution_data:

            # Specify, for each sample, what molecule data to generate (if any)

            # Samples here should correspond to the 'samples' entry above

            '1':

                # Number of packets to generate.

                # If 0, don't generate any and just use the provided molecule files

                num_packets: 1

 

                # Number of molecules to generate per packet.

                # Increasing this number increases the amount of work done per job.

                # Generally want to keep this under 200000 to not run out of memory.

                # If memory usage errors occur, reducing this can help.

                num_molecules_per_packet: 200

 

                # Directory containing CAMPAREE output for the sample

                # This data is used to generate the resulting molecules

                sample_data_directory: input_data/sample1/

            '2':

                num_molecules_per_packet: 200

                num_packets: 1

                sample_data_directory: input_data/sample2/

 

    steps:

        # The 'steps' list specifies all the steps that the Library Prep Pipeline

        # will run. To specify a custom step, include it here.

        # Python files with the corresponding code are looked up according to

        # the pattern {module_name}.{class_name} pattern, so that the

        # {module_name}.py file will be checked for {class_name} class

 

    -   step_name: polya_step.PolyAStep

        # PolyA selection step removes RNA that does not have an poly A tail

        # of sufficient length. Enriches for mRNA.

        parameters:

            # Chance per base of fragmentation prior to selection

            # Increasing this above 0 induces a 3' bias.

            # A value of 0.001 induces a reasonably high 3' bias

            breakpoint_prob_per_base: 0.0

 

            # Probability of retention is computed as a value between

            # min_retention_prob and max_retention_prob

            # For every base of the polyA tail beyond min_polya_tail_length

            # the probability increases linearly from min_retention_prob

            # by length_retention_prob, up to a max of max_polya_tail_length.

            max_retention_prob: 1.0

            min_retention_prob: 0.0

            min_polya_tail_length: 40

            length_retention_prob: 0.05

 

    -   step_name: fragment_step.FragmentStep

        # Fragmentation step breaks each molecule into pieces

        parameters:

            # Fragmentation methods are available.

            # 'uniform' fragments at each base equally and is the default

            method: uniform

            # This 'rate' is the rate parameter for an exponential distribution

            # which determines the time it takes until a molecule to fragments.

            # Molecules which would take longer then 'runtime' to fragment, do not fragment.

            # Molecules may also fragment multiple times if rate is high enough.

            rate: 0.005

            runtime: 1

            # Since fragmentation generates many very small fragments that will

            # be quickly lost in the following steps, we have the option to

            # drop those fragments immediately. Setting this value can significantly

            # decrease runtime and memory requirements.

            min_frag_size: 20

 

            # If method == 'beta', then the fragmentation sites can be biased within

            # the molecule, according to a beta distribution.

            # NOTE: using 'beta' can generate unusual coverage plots since there are

            # significant edge effects around the transcript. However, it can be used

            # to give a more realistic fragment distribution.

            #

            # If using 'beta', you must also specify the following:

            #

            # The parameters for the beta distribution. Set these so that

            # A = B > 0 to bias towards fragmentation in the middle of the fragment,

            # with larger values biasing further towards the middle.

            # If A > B, then would bias towards the 5' end instead.

            # beta_A: 3.0

            # beta_B: 3.0

            #

            # And the N factor allows a non-linear fragmentation rate depending

            # upon the length of the molecule. Values >1 indicate that longer molecules

            # are more likely to fragment than smaller ones, biasing towards larger

            # fragment sizes.

            # beta_N: 2.0

 

    -   step_name: first_strand_synthesis_step.FirstStrandSynthesisStep

        # First Strand Synthesis performs (hexamer) priming on the RNA

        # and then generates the cDNA from this. Priming sites can be

        # biased according the position probability matrix (PPM) that

        # specifies the weighting of each base according to the position

        # relative to the first base of the primer.

        # The 5' most primer is then used to extend the DNA and create

        # the cDNA, thereby potentially losing some of the 5' end of

        # the molecule.

        parameters:

            # If 'perfect_priming' is true, then we always prime

            # exactly on the 5' most end of the molecule and the

            # cDNA matches perfectly the original molecule

            # If true, all other parameters are ignored.

            perfect_priming: false

 

            # PPM matrix describing biases of the priming sites

            # Values must sum to 1 in each base. Length of the PPM

            # determines the primer lengths.

            # See Kasper et all, 2010 PMID: 20395217 for example

            # observation of these biases in fragment start positions.

            # If no bias is desired, set all values to 0.25.

            position_probability_matrix:

                A: [0.50, 0.1, 0.40, 0.30, 0.25, 0.15]

                C: [0.20, 0.5, 0.3 , 0.25, 0.25, 0.15]

                G: [0.15, 0.1, 0.15, 0.25, 0.25, 0.20]

                T: [0.15, 0.3, 0.15, 0.20, 0.25, 0.50]

 

            # Rate of priming: approximate number of primers that

            # will bind (if bias-free) per kilobase of the fragment

            # Higher numbers will loose less of the 5' side.

            primes_per_kb: 50

 

    -   step_name: second_strand_synthesis_step.SecondStrandSynthesisStep

        # Second Strand synthesis copies the single-stranded cDNA from

        # the First Strand Synthesis step into double-stranded cDNA.

        # All parameters are same as the First Strand Synthesis step.

        # If not perfect_priming, will loose some of the 5' end of the

        # first strand of cDNA, which is the 3' end of the original fragment

        parameters:

            perfect_priming: false

            position_probability_matrix:

                A: [0.50, 0.1, 0.40, 0.30, 0.25, 0.15]

                C: [0.20, 0.5, 0.3 , 0.25, 0.25, 0.15]

                G: [0.15, 0.1, 0.15, 0.25, 0.25, 0.20]

                T: [0.15, 0.3, 0.15, 0.20, 0.25, 0.50]

            primes_per_kb: 50

 

    -   step_name: sizing_step.SizingStep

        # Sizing step throws out molecules that do not fit the

        # desired size distribution

        parameters:

            #    Probability of retention by length

            #                   __________________

            #    1|         ___/                  \___

            # p   |    ____/                          \____

            #    0|___/                                    \____

            #     -----------------------------------------------

            #         |          |               |          |

            #     min_length     |               |      max_length

            #                    |   select_all  |

            #                  start            end

            #

            # All molecules with lengths outside of this range

            # will be discarded:

            min_length: 100

            max_length: 400

            # Molecules whose length is inside this range will

            # always be retained. Retention probability raises

            # linearly from 0 to 1 between the min/max values

            # above and this range here.

            select_all_start_length: 200

            select_all_end_length: 300

 

    -   step_name: adapter_ligation_step.AdapterLigationStep

        # Adpater Ligation attaches adapters to each end of the molecules

        # These are used for the PCR step to initiate PCR and include the

        # sample identifying barcodes.

        #

        # This uses the adapters specified in the resources section

        # as well as the sample i5/i7 barcodes in the samples section

        parameters: {}

 

    -   step_name: pcr_amplification_step.PCRAmplificationStep

        # PCR amplification step creates many copies of each molecule

        # Since the number of molecules generated is quite large and

        # since in real sequencing, only a small fraction of the molecules

        # prepped end up forming clusters on the flowcell and being

        # sequenced, we give the option to remove most of the PCR

        # generated molecules prior to the end of the step.

        # This provides a very large increase in speed and decrease

        # in memory usage.

        parameters:

            # Number of cycles to use

            # Generates 2^n fragments for each input fragment

            number_cycles: 10

 

            # Retain only this percent of the data

            # NOTE: scale is 0-100, not 0-1, so 0.08 means 0.0008 of the generated molecules are kep

            # This value can be approximated from real data with UMI tags by the following:

            # Compute the PCR duplication rate using the UMI tags.

            # Let N be the number of PCR steps. Assuming that each molecule generates 2^N

            # PCR descendants, if all were sequenced, we would expect all molecules to be duped

            # 2^N times. Instead, choose retention_percentage with the following:

            #

            # '''

            # import scipy.stats

            # pcr_rate = scipy.stats.binom(p=retention_percentage / 100, n=2^N).sf(1)

            # # Choose retention_percentage so that pcr_rate is approximately the observed PCR rate

            # '''

            #

            # For example, retention_percentage = 0.08 and number_cycles=10 gives dupe rate

            # of about 20%, which is pretty typical. If number_cycles changes, this should

            # be modified too to maintain dupe rate.

            retention_percentage: 0.08

 

            # Induce a GC bias by discarding some PCR duplicates

            # according to their GC bias. All fragments overall

            # GC content is computed and then the following three

            # parameters are used to compute a probability of retention

            # gc_bias_constant + gc_bias_linear*(gc - 0.5) + gc_bias_quadratic*(gc - 0.5)^2

            # clipped to always be within 0 and 1

            # For no bias, set to 1, 0, and 0 for const, linear, and quadratic, respectively.

            # To bias towards GC=0.5 content, set gc_bias_quadratic to be negative,

            # such as -100 for a large GC bias.

            gc_bias_constant: 1.0

            gc_bias_linear: 0.0

            gc_bias_quadratic: -100

 

            # During PCR, we have some chance of mis-copying the molecules

            # These are specified here as probability per-base, so 0.001 is one error per 1kb

            deletion_rate: 0.0001

            insertion_rate: 0.0001

            substitution_rate: 0.001

 

sequence_pipeline:

    # The Sequence Pipeline takes clusters bound to the flowcell and turns them into

    # final, sequenced results, including fastq files or SAM files with the 'true'

    # alignments of the simulated molecules to the reference genome.

    steps:

    -   step_name: bridge_amplification_step.BridgeAmplificationStep

        # Bridge amplification takes a seed of a cluster and generates

        # many copies of the molecule to form a cluster.

        parameters:

            # Number of cycles to perform. Generates 2^N molecules in the cluster

            # via a PCR-like bridge amplification process.

            cycles: 10

            # Substitution rate per base in the cluster formation

            # No indels are generated here, for computation simplicity.

            substitution_rate: 0.01

 

    -   step_name: sequence_by_synthesis_step.SequenceBySynthesisStep

        # Sequencing by Synthesis then synthesizes complementary strands to each

        # molecule in the clusters with flourescence readings at each base.

        # If some of the molecules have errors in them (from bridge amplification),

        # the flourescence will be imperfect. Moreover, errors also happen in this step

        # including phasing, where some molecules get ahead or behind in the sequencing

        # and so are displaying the wrong base.

        # From these readings, the final sequenced values are generated.

        parameters:

            # Determines which read is the forward and which is the reverse read

            # NOTE: currently only 'true' is implemented

            forward_is_5_prime: true

            # Whether to sequence both ends or just one

            # NOTE: currently only 'true' is implemented

            paired_ends: true

            # Length of the reads to generate, in bases

            read_length: 100

            # The rate of phasing, either forward (skip) or backwards (drop)

            # per base per molecule

            skip_rate: 0.002

            drop_rate: 0.002

 

    flowcell:

        # BEERS simulates one flowcell

        # Layout of the flowcell is defined here

 

        # Coordinate strategy defines how coordinates are assigned

        # Currently only 'random' is supported

        coordinate_strategy: random

 

        # Flowcell geometry defines the dimensions of the flowcell

        # Example is given for a typical HiSeq 2500 flowcell

        flowcell_geometry:

            # Range of the lane numbers

            min_lane: 1

            max_lane: 8

            # Range of the tile values

            max_tile: 2228

            min_tile: 1101

            # Range of the x-coordinates

            min_x: 1012

            max_x: 32816

            # Range of the y-coordinates

            min_y: 998

            max_y: 49247

        # List of lanes that samples will be written into

        # Molecules are distributed evenly across lanes

        # These must all be within [min_lane, max_lane] range

        lanes_to_use: [1, 2]

 

 

パラメータが決まったら実行する。

cd baby_mouse_example/
#BEERS2のラン
beers --configfile baby.config.yaml --jobs 1

テスト時は1分以内に終了した。

出力例

デフォルトでは、fastqの他にゴールデンスタンダードなsamも出力される("output_sam: true"の部分)。

 

引用

BEERS2: RNA-Seq simulation through high fidelity in silico modeling 

Thomas G Brooks, Nicholas F Lahens, Antonijo Mrčela, Dimitra Sarantopoulou, Soumyashant Nayak, Amruta Naik, Shaon Sengupta, Peter S Choi, Gregory R Grant

Briefings in Bioinformatics, Volume 25, Issue 3, May 2024

タンパク質言語モデルの配列表現の直接比較に基づいて遠隔相同性検出を行う pLM-BLAST

 

 配列比較による相同性の検出は、タンパク質の機能と進化の研究における典型的な最初のステップである。この研究では、タンパク質言語モデルのこのタスクへの適用可能性を探る。pLM-BLASTはBLASTにインスパイアされたツールであり、タンパク質言語モデルProtT5から得られた単一配列表現(embeddings)を比較することにより、遠方相同性を検出する。本著者らのベンチマークにより、pLM-BLASTはHHsearchと同程度の精度を維持し、非常に類似した配列(50%以上の同一性)および著しく分岐した配列(30%未満の同一性)の両方において、有意に高速であることが明らかになった。さらに、pLM-BLASTはローカルアラインメントを計算する能力により、他の埋め込みベースのツールの中で際立っている。pLM-BLASTによって生成されたローカルアラインメントが、多くの場合、高度に乖離したタンパク質を連結することを示し、それによってpLM-BLASTが、これまで発見されていなかった相同関係を発見し、タンパク質アノテーションを改善する可能性を強調する。

pLM-BLASTは、MPI Bioinformatics Toolkitを介して、事前計算されたデータベースを検索するためのウェブサーバーとしてアクセスできる (https://toolkit.tuebingen.mpg.de/tools/plmblast)。また、カスタムデータベースを構築し、バッチ検索を実行するためのスタンドアロンツールとしても利用できる (https://github.com/labstructbioinf/pLM-BLAST)。

 

webサービス

https://toolkit.tuebingen.mpg.de/tools/plmblastにアクセスする。

 

配列を選択してサブミットする。

 

出力例

 

ローカルマシンへのインストール

Github

mamba create --name plmblast python=3.10 -y
conda activate plmblast
git clone https://github.com/labstructbioinf/pLM-BLAST.git
cd pLM-BLAST/
pip install -r requirements.txt

python embeddings.py -h

$ python embeddings.py -h

usage: embeddings.py [-h] {start,resume} ...

 

Embedding script create embeddings from sequences via desired embedder

by default `seq` column in used as embedder input. Records are stored

as list maintaining dataframe order. 

In python load via: 

>>> import torch

>>> torch.load(..)

or 

>>> import pickle

>>> with open(.., 'rb') as f:

>>>     embs = pickle.load(f)

example use:

python embeddings.py start data.csv data.pt -cname seqfull

# for fasta input

python embeddings.py start data.fasta data.pt

# for file per embedding output

python embeddings.py start data.fasta data --asdir

# resume interrupted calculations

python embeddings.py resume data

 

options:

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

 

options:

  {start,resume}

    start         starting new calculations

    resume        continue calculations from checkpoint, checkpoint is automatically created and stored as emb_checkpoint.json in output directory this is only available when using asdir flag

 

 

データベース

wget http://ftp.tuebingen.mpg.de/pub/protevo/toolkit/databases/plmblast_dbs

 

実行方法

デフォルトではembeddings.pyはProtT5(-embedder pt alias for Rostlab/prot_t5_xl_half_uniref50-enc)を使用する。(pytouchが)GPUを認識するなら--gpuをつける。

python embeddings.py start database.fasta database -embedder pt --gpu -bs 0 --asdir
  • 大きなデータベースを扱う場合、embeddings.pyを実行する前に入力配列を長さでソートする。

 

論文より

  • TM-vecやDEDALとは対照的に、pLM-BLASTは教師なしアプローチに基づいており、特別なディープラーニングモデルのトレーニングや、構造的に類似または相同なタンパク質ペアに基づくポジティブラベルの定義を必要としない。さらに、グローバルアラインメントのみを提供するTM-vecやEBAのような手法とは異なり、pLM-BLASTはローカルアラインメントとグローバルアラインメントの両方を計算する能力を持ち、これは短いサブドメイン断片の保存に依存する可能性のある遠い相同性関係を決定するために不可欠である(Alva et al. 2015, Kolodny et al)。

  • ECODタンパク質分類データベース(Cheng et al. 2014)のドメインペアを用いたベンチマークにおいて、pLM-BLASTは高速なペアワイズ比較とデータベース検索に適しており、遠方の相同性を検出し、高精度のアラインメントを生成することが実証された。

引用

pLM-BLAST: distant homology detection based on direct comparison of sequence representations from protein language models 

Kamil Kaminski, Jan Ludwiczak, Kamil Pawlicki, Vikram Alva, Stanislaw Dunin-Horkawicz Author Notes

Bioinformatics, Volume 39, Issue 10, October 2023

 

関連

MPI Bioinformatics ToolkitのPSI-BLASTサービス