macでインフォマティクス

macでインフォマティクス

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

バイオインフォのツールが動かない問題を試行錯誤しながら直す

2021 9/23 誤字修正

 

 複数の方からsamtoolsやbowtie2が動作しないという連絡をいただいたので、今日はその対処方法を検討します。通常、このような問題をいただいた場合、調べる側としては、その方が該当するツールをどのような環境にどのような方法で導入したのかにまず注目するかと思います。今回のケースでは、連絡をいただいた方の1名は、インテル版のmacにパッケージマネージャのcondaを使ってsamtoolsを導入したところ、エラーが発生したという事でした。condaコマンドにはM1 maclinux版もありますが、インテルmacのcondaコマンドでは、ローカルのライブラリのチェックをした上で、macos向けにコンパイルされたプログラムのバイナリと足りないライブラリを導入(もしくはアップデート)しようとします(*1)。回りくどい書き方をしましたが、condaを使って導入したなら、バイナリの実行コード自体はOS互換と思われます(*2)。

 導入方法自体は問題なさそうなので、あとはエラー情報を見て考えることにします。samtoolsのコマンドの実行時にどんなエラーが出たのかエラーメッセージを送ってもらいました。以下のメッセージが出たそうです。

$ samtools
samtools: error while loading shared libraries: libcrypto.so.1.0.0: cannot open shared object file: No such file or directory

 

これは共有ライブラリのlibcrypto.so.1.0.0が参照できないというエラーです(.soはshared.libraryの略。macでは.dylib)。自分も昔condaを使ってsamtoolsやbcftoolsをインストールした時に発生した記憶があります。condaで導入すると通常ライブラリを共有するdynamic linked(説明)のバイナリとしてインストールされるので、インストール時、samtoolsが依存している共有ライブラリのlibcrypto.so.1.0.0が正しく導入されなかったのだと推測されます。

ではここからは自分の環境でエラーを再現できるか試していきます。まずは自分のローカル環境にcondaコマンドでインストールしたsamtoolsを検索します。私のように色々なツールを導入している場合、(samtoolsが様々なツールで使われているため)複数のsamtoolsバイナリがヒットするはずです。

find ~/miniconda/ -name samtools

出力は省略しますが、確かに、自分のマシンの~/miniconda/以下から複数のsamtoolsバイナリがヒットしました。そこで、lddコマンドを使ってそれぞれのsamtoolsで共有ライブラリは全て見えているか順番に調べていったところ、version1.7.0のsamtoolsで同様の現象が発生していました。

#lddコマンドで共有ライブラリを調べる(例;ldd /bin/cat)
ldd <path>/<to>/samtools1.7

f:id:kazumaxneo:20210923171648p:plain

 

このバージョンのsamtoolsを、質問された方と同じcondaコマンドを使って導入します。いつものように、ubuntu18で検証用の仮想環境を立ち上げて行います。conda互換のmambaを使います。pythonのバージョンは3.9です。

mamba create -n samtools_env -y
conda activate samtools_env
mamba install -c bioconda samtools==1.7.0 -y

 

動作するか確認します。

samtools

samtools: error while loading shared libraries: libcrypto.so.1.0.0: cannot open shared object file: No such file or directory

version1.7のsamtoolsをcondaで入れると同じエラーが再現されました。共有ライブラリのlibcrypto.so.1.0.0が見つからないというエラーです。

lddコマンドを使ってこのsamtoolsの共有ライブラリが見えているかどうか確認します。

ldd $(which samtools)

f:id:kazumaxneo:20210923171648p:plain

samtoolsの動作に必要な共有ライブラリの中で、libcrypto.so.1.0.0だけライブラリのパスが見えていません(写真は上と同一のもの)。

 

ldconfigコマンドにオプションの"-p"をつけると、システム側から見えている全ての共有ライブラリを確認できます。

> sudo ldconfig -p

しかしこれだと多すぎて見にくいので、grepを使って"libcrypto.so"があるか調べます。

sudo ldconfig -p |grep libcrypto.so

f:id:kazumaxneo:20210923180218p:plain

libcrypto.so.1.1は見つかりましたが、やはりlibcrypto.so.1.0.0は見つかりません。

さて、ここからどうするかですが、マイナーバージョンアップの例えばlibcrypto.so.1.0.1があれば、シンボリックリンクでlibcrypto.so.1.0.0として、これを/usr/libなどに追加すれば認識するかもしれません。しかし、ここではもう少し粘ります。まず、/etc/ld.so.cacheに含まれない場所に見つからないか、念のため確認します(参考ページ)。検索対象は、各ユーザーが共通して利用するライブラリやプログラムのバイナリを置く/usr(ユーザーディレクトリ、参考ページ)とします。

find /usr -type f -name "libcrypto.so.1.0*"

(出力は省略)

 

見つかりませんでした。どうしましょうか?このライブラリ名でググったところ、libcrypto.so.1.0はオープンソースの暗号通信プロトコルであるOpenSSL(リンク)でも使われていました。バイオインフォの情報共有サイトでは、このことを利用して、condaでOpenSSL(version1.0)をインストールしてsamtoolsのエラーを回避する方法が提案されていました(リンク)。OpenSSLの依存だけcondaで導入するなら

> mamba install -c bioconda samtools openssl=1.0 -y --only-deps

 

しかし、この方法を試したところ、私の場合は、他のツールとの依存関係の問題で導入できませんでした。そこで、OpenSSLの配布ファイルをbiocondaから直接ダウンロードします。

リンク

https://anaconda.org/conda-forge/openssl/files?page=5

linux-64/openssl-1.0.2h-0.tar.bz2をダウンロードして解凍し、lib/のlibcrypto.so.1.0.0を確認します。

f:id:kazumaxneo:20210923191056p:plain

このlibcrypto.so.1.0.0をcondaの共有ライブラリの配置場所である~/miniconda/lib/(ここでは~/miniconda/env/samtools_env/lib/)、もしくは/usr/lib/x86_64-linux-gnu/に配置します。

> cp lib/libcrypto.so.1.0.0 ~/miniconda/env/samtools_env/lib/

 

コピーしたらsamtoolsがランできるか確認します。

>samtools --version

f:id:kazumaxneo:20210923204801p:plain

バージョンが表示されます。ランできるようになりました。

 

今回はcondaでツールを管理していたため、LD_LIBRARY_PATHに追加することは避けました。自分でソースコードからビルドする場合、LD_LIBRARY_PATHに共有ライブラリのパスを共有ライブラリの優先度を考慮して追加し、"ldconfig"して反映させれば良いのではないでしょうか(参考サイト)。

 

いつになるか未定ですが、次回はperlのライブラリエラーについて考えます。

参考


 

*1

2021年9月現在、condaコマンドでは、numpyなどを除いてほとんどのバイオインフォマティクスのツールはarmのM1 macosにはネイティブにはインストールできません(導入時に弾かれるはず)。M1 macバイオインフォマティクスのツールをネイティブに動かすには、有志の方が公開されているバイナリを探してダウンロードするか (例: bwaのstatic binary)、もしくは自分で頑張ってビルドするということになります(参考ページ)。その一方で、インテルCPUが搭載されたmacは、人気のバイオインフォマティクスのツールの多くがインストールできます。macにこだわる以上は、インテルCPUが搭載されたmacの方が現状有利です (armのM1 macのdockerはrosetta2に依存しており(docker)、現在は速度面でかなり不利)。

 

*2

致命的なバグがあったり、対応するライブラリのバージョンがなかったりするバージョンは、コミュニティへの情報提供を通じて削除されます。例えばbowtie2の2.2.3はbioconda channelには見当たりません(biostarsスレッド)。

f:id:kazumaxneo:20210923152226p:plain

 

 

16S rRNA塩基配列データから超可変領域を抽出する HyperEx

 

 16SリボソームRNA遺伝子は、生物学において最も研究されている遺伝子の一つである。この16SリボソームRNAの重要性は、細菌や古細菌の系統学や分類学上の解明に広く応用されていることによる。実際、16SリボソームRNAは、ほとんどすべてのバクテリア古細菌に存在し、多くの有用な特性に加え、低い変異率が特徴にある。16SリボソームRNAは、9つの超可変領域から構成されており、ハイスループットシーケンシング技術を用いた同定や、メタバーコーディング研究などのコミュニティ研究では、これらの領域が一般的にターゲットとなる。残念ながら、超可変領域は、すべての細菌の分類学上の分解能が同じではない。このため、特定の研究で対象とする最適な超可変領域を決定するために、事前にインシリコ解析を行う必要がある。しかし、著者らの知る限り、16S rRNAシーケンスデータから超可変領域を抽出するプライマーベースの自動化されたオープンソースツールは存在しない。ここでは、埋め込まれたプライマーまたはユーザーに与えられたプライマーに基づいて、目的の超可変領域を効率的に抽出するHyperExを紹介する。HyperExは、正確なペアワイズシーケンスアラインメントのためのMyersアルゴリズムを実装している。HyperExは、オペレーティングシステムに依存しないRustコマンドラインツールとして、MITライセンスの下、https://github.com/Ebedthan/hyperex および https://crates.io から自由に入手することができる。

 

インストール

cargoでインストールできなかったので、リリースよりプリビルドされたバイナリをダウンロードしてテストした。

  •  

Github

cargo install hyperex

#From source
git clone https://github.com/Ebedthan/hyperex.git
cd hyperex
cargo build --release
cargo test

> ./hyperex -h

f:id:kazumaxneo:20210922115924p:plain

 

 

実行方法

16SrRNAの配列を指定する。gzip圧縮ファイルにも対応している。

hyperex file.fa

出力

hyperex_out.gff

f:id:kazumaxneo:20210922121506p:plain

hyperex_out.fa

f:id:kazumaxneo:20210922121542p:plain

logも出力される。

 

可変領域を指定する。

hyperex --region v3v4 file.fa.xz

 

引用

HyperEx: A Tool to Extract Hypervariable Regions from 16S rRNA Sequencing Data
Anicet Ebou,  Dominique Koua,  Adolphe Zeze

bioRxiv, Posted September 05, 2021

 

機械学習の手法でエミュレートされたBWA-MEM: BWA-MEME

 

  次世代シーケンサーの普及やシーケンサースループットの向上に伴い、効率的なショートリードのアライメントが求められているが、その中でもシーディングは主要な性能ボトルネックの一つとなっている。Seeding phaseのキーとなるチャレンジは、リファレンスDNA配列中からショートリードの部分配列の正確な一致を検索することである。しかし、既存のアルゴリズムでは、頻繁にメモリにアクセスするため、性能に限界があった。

 本論文では、seedingを効率的に行うための完全一致検索問題を解決するために、学習したインデックスを活用した初の本格的なショートリードアライメントソフトウェアであるBWA-MEMEを紹介する。BWA-MEMEは、seeding phaseで多用されるSMEM探索に学習済みインデックスを活用するという課題を解決する、接尾辞配列探索アルゴリズムをベースとした実用的かつ効率的なseedingアルゴリズムである。評価の結果、BWA-MEMEは、BWA-MEM2と同一のSAM出力を確保しつつ、命令数を4.60倍、メモリアクセスを8.77倍、LLCミスを2.21倍に削減することで、BWA-MEM2と比較してシーディングのスループットを最大で3.45倍に高速化した。

(中略)

BWA-MEMEは、新しい学習済みインデックス構造とアルゴリズムを導入することで、これらの課題に対処する。まず、部分的に3層の再帰的モデルインデックス(P-RMI)を提示し、接尾辞の不均衡な分布にうまく適応し、正確な予測を行う。この過程で、入力された部分文字列や接尾辞を数値キーにエンコードするアルゴリズムを設計する。この数値キーはP-RMIへの入力として与えられ、P-RMIは接尾辞配列内の予測位置と予測の誤差範囲を提供する。P-RMIは、接尾辞配列内の予測位置とその誤差範囲を提供し、誤差範囲内でバイナリ検索を行い、部分文字列が一致する参照位置を見つける。(以下略)

 

 

インストール

  • 3990x CPUを積んだUbuntuマシンでビルドした(make -j arch=avx2)。

Github

#from source
git clone https://github.com/kaist-ina/BWA-MEME.git BWA-MEME
cd BWA-MEME

# If AVX512BW (512-bit SIMD) is supported
make clean
make -j<num_threads> arch=avx512

# If AVX2 (256-bit SIMD) is supported
make clean
make -j arch=avx2

# If SSE4.1 (128-bit SIMD) is supported (default)
make -j<num_threads>

> ./bwa-mem2.avx2 

Usage: bwa-mem2 <command> <arguments>
Commands:
  index         create index
  mem           alignment
  version       print version number

 

 

実行方法

レポジトリではヒトゲノムのアセンブリ配列をダウンロードして、indexを付け、訓練してマッピングするまでの流れが書かれています。ここでもその通り勧めます。

 

1、ダウンロード

wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/human_g1k_v37.fasta.gz
gunzip human_g1k_v37.fasta.gz

 

2、indexing

bwa-mem2.avx2 index -a meme -t32 human_g1k_v37.fasta

 

3、P-RMIのトレーニン

./build_rmis_dna.sh human_g1k_v37.fasta

 

4、マッピング

bwa-mem2 mem -Y -K 100000000 -t32 -7 human_g1k_v37.fasta input_1.fastq -o out.sam
  • -Y    use soft clipping for supplementary alignments
  • -K    process INT input bases in each batch regardless of nThreads (for reproducibility)
  • -7    Use Learned index for seeding
  • -t    number of threads [1]
  • -o   Output SAM file name

 

この後、レポジトリでは、オリジナルのbwa-mem2も走らせ、diffコマンドを使ってsam形式ファイルに違いがないか確認する流れについて書かれています。

引用

BWA-MEME: BWA-MEM emulated with a machine learning approach

Youngmok Jung,  Dongsu Han

bioRxiv, Posted September 01, 2021

 

関連


elPrep 5を使ったバリアントコール

 

 GATK Best Practices for variant callingに完全対応したelPrep5 (紹介) には、大きく分けて2つのモードが用意されています。1つ目は完全にRAM内で動作する(フィルタ)モードで、これは中間ファイルを全く書き出さず完全にRAM内で計算を進めるため、非常に高速です(GATK4で順番に進める場合の11−15倍高速)。2つ目は、ユーザーが提供したsam/bamファイルから実行するsfmモードで、このsfmモードでは、前もって中間ファイルであるsamやbamを作っておく必要がありますが、染色体領域ごとにデータを分割して進めるため、完全にRAM内で動作するフィルタモードよりもメモリ使用量を抑えることができます(GATK4で順番に進める場合の6-7.5倍高速)。利用可能な計算資源の量に対応して使い分けできるようになっているという事になりますGIthubの説明によると、sfmモードは、ヒトのWGS30xで128GBのメモリ使用量でランでき、ヒトWES30xでは24GBメモリ使用できるようです( WES30xのメモリモードだと96GGB使用)。elPrep5の論文ではm5.24xlargeインスタンスリンク)というかなりリッチな環境を使ってベンチマークなどを行なっていますが、このsfmモードなら、より少ない計算資源でも安定してランできそうです。

 

インストール

論文が出たからだと思いますが、ソースコードが公開され、condaで導入できるバージョンもv5になっています。

#anaconda (link)
mamba install -c bioconda elprep -y

 

 

実行方法

1、BQSR(参考サイト)を行うには、補正部位から除外する既知バリアントサイトのリスト (vcf or bed)も提供する必要がある。このvcf またはbed形式のファイルは、予めelprep vcf-to-elsites コマンドを走らせてelprepの形式に変換しておかないと使用できない。

ここでは時間の関係でdbSNP common SNV(生殖細胞由来、少なくとも1つの主要な集団においてマイナーアレル頻度(MAF)が>>0.01、少なくとも2人の血縁関係のない人がマイナーアレルを持っている)を使う。elprep vcf-to-elsitesコマンドは、入力のVCF、出力ファイルの順番に指定する。

elprep vcf-to-elsites dbSNP_common_all.vcf.gz  dbSNP_common_all.elsites

このコマンドは1回だけ実行すれば良いのですが、VCFが巨大だとかなりメモリを使います(最近のバージョンのdbSNP allだと300GB以上)。注意して下さい。

 

2、BQSRする際に必要なリファレンスfastaファイルも、elPrepの形式に変換しておく必要がある。elprep fasta-to-elfastaコマンドを使う。

elprep fasta-to-elfasta ucsc.GRCh38.fasta GRCh38.elfasta

 

3、elPrep5をランする。ここではメモリ使用量を抑えられるsfmモードを使う。sfmモードでは、ユーザーの方でアラインされたリードのsam かbam形式のファイルを提供する必要がある(論文中ではbwa memを使用)。

elprep sfm input.bam output.bam \
--mark-duplicates --mark-optical-duplicates NA12878.output.metrics\
--sorting-order coordinate \
--bqsr output.recal --known-sites dbSNP_common_all.elsites \
--reference GRCh38.elfasta --haplotypecaller output.vcf.gz \
--target-regions targetedregions.bed \
--filter-mapping-quality 0 \
--nr-of-threads 12 \

必要なら"--filter-mapping-quality 0"オプションも使います。 exome解析では、"--target-regions"オプションでbed形式ファイルを指定します。"filter-non-overlapping-reads”オプションを立てると、captureされた領域のbedファイルにマッチするリードだけがbamに保存され、バリアントコールの対象になります。バリアントコールは行わずにキャプチャ領域のリードだけ取り出したい時にも使えます。

注意

出力されるVCFはほぼすべてのポジションを含んでいます。バリアントコールに使うには、この後でフィルタリングする必要があります。

 

elPrep 5の論文では、NA12878を例にelPrep 5の使い方についても説明しています。また、GATKのオリジナルの Best Practicesのフローと比較して、リソース使用量を議論しています。簡単にまとめると、elPrepを使うことでディスク使用量、計算時間などを大幅に抑えられ、サーバーのレンタルコストを節約できるという結果です。論文にアクセスしてみてください。

引用

Multithreaded variant calling in elPrep 5

Charlotte Herzeel  ,Pascal Costanza ,Dries Decap,Jan Fostier,Roel Wuyts,Wilfried Verachtert

PLOS ONE, Published: February 4, 2021

 

参考


関連

 

ロングリードを使ったSNVとSVのフェーシングを行う LongPhase

 

 ロングリード・フェーシングは、二倍体ゲノムの再構築、バリアント・コーリングの改善、メタゲノミクスにおける微生物株の解決などに用いられてきた。しかし、既存の手法では、大きな構造変化(Structural Variation: SV)によって位相差ブロックが破壊されてしまい、集団規模のフェーシングを行うには効率が満足できない。本論文では、超高速アルゴリズムLongPhaseを紹介する。このアルゴリズムは、ヒトゲノムの一塩基多型(SNP)とSVを同時に、最先端のWhatsHapやMarginの10倍の速さである約10-20分でフェーシングすることができる。特にLongPhaseは、ロングリード(N50=26Mbp)のみで、ほぼ染色体レベルのはるかに大きなフェーズブロックを生成する。LongPhaseとNanoporeの組み合わせは、追加のトリオ、染色体の構造、シングルセルのストランド・セクデータを必要とせずに、染色体レベルのフェージングを提供するコスト効率の高いアプローチであることを実証している。

 

インストール

htslibライブラリをビルドしてからGithubの手順に従ってビルドした。

依存

  • LongPhase relies on htslib for parsing BAMs. We package the htslib folder into our project. 

Github

リリース

git clone https://github.com/twolinin/LongPhase.git
cd LongPhase
autoreconf -i
./configure
make -j 4

> ./longPhase phase

phase: missing arguments. --ont or --pb

phase: missing SNP file.

phase: missing bam file.

 

Usage:  phase [OPTION] ... READSFILE

      --help                     display this help and exit.

      --dot                      each contig/chromosome will generate dot file. 

      --ont, --pb                ont: Oxford Nanopore genomic reads.

                                 pb: PacBio HiFi/CCS genomic reads.

      --sv-file=NAME             input SV vcf file.

      -s, --snp-file=NAME        input SNP vcf file.

      -b, --bam-file=NAME        input bam file.

      -o, --out-prefix=NAME      prefix of phasing result.

      -r, --reference=NAME       reference fasta.

      -t, --threads=Num          number of thread. 

      -d, --distance=Num         phasing two variant if distance less than threshold. default:300000

      -c, --crossBlock=Num       each block tries to connect with next N blocks. default:1

 

 

実行方法

ランにはbamとbam.bai、fastafasta.fai、そしてVCFファイルが必要。

#ONT
LongPhase phase \
-s SNP.vcf -b alignment.bam -r reference.fasta -t 8 \
-o output_prefix --ont

#PacBio
LongPhase phase \
-s SNP.vcf -b alignment.bam -r reference.fasta -t 8 \
-o output_prefix --pb

 

リリースからテストデータ(210GBほどある)が公開されている。

#SNV
LongPhase phase -s SNP_PEPPER.vcf -b HG002_ULR_60x.bam -r GRCh37.fa \
-t 20 -o output_prefix --ont

#SV and SNV
LongPhase phase -s SV_sniffles.vcf -b HG002_ULR_60x.bam -r GRCh37.fa \
-t 20 -o output_prefix --ont

出力ファイル

SNPとSVのco-phasingを行った場合、2つのVCF(SNP用とSV用)が出力される。GTフィールドにはフェーズされたSVが、PSフィールドにはブロックIDが格納される。Githubに例があるので興味がある方は確認して下さい。

 

引用

LongPhase: an ultra-fast chromosome-scale phasing algorithm for small and large variants

Jyun-Hong Lin, Liang-Chi Chen, Shu-Qi Yu, Yao-Ting Huang

bioRxiv, Posted September 11, 2021

 

関連


 

公開されているプラスチドゲノムのアノテーションとinverted repeatsを調べる airpg

2021 9/18 使い方を理解していなかったので一旦コマンドは消去

 

 ほとんどの顕花植物では、プラスチドのゲノムは、大小のシングルコピーと2つの逆方向反復配列領域からなる4分割構造をしている。近年、何千ものプラスチドのゲノムが配列決定され、公的な配列リポジトリに登録されている。その中でも特に、逆方向反復配列領域の長さや位置を指定するアノテーションの質には問題があることが知られている。しかし、多くの生物学的研究では、公開されているプラスチドのゲノムを額面通りに使用し、暗黙のうちにその配列アノテーションの正しさを前提としている。

 公表されているプラスチドゲノムの中で、逆方向反復配列のアノテーションが不完全または間違っている頻度を自動的に評価するPythonパッケージ、airpgを紹介する。具体的には、NCBI Nucleotideから様々な検索条件でプラスチドのゲノムを自動的に検索し、逆方向反復配列の長さと位置を調査し、ゲノム配列の自己比較によって逆方向反復配列のアノテーションを確認する。また、重複するゲノムレコードを自動的に識別・削除する機能や、純粋に逆方向反復配列を持たない分類を考慮する機能も備えている。airpgを用いて、2020年末までにNCBI Nucleotideに登録されたflowering plantsの全プラスチドゲノムにおける逆方向反復配列アノテーションの有無を調査し、レコードのメタデータとの関連性を統計的に分析した結果、ゲノムレコードのリリース年と出版状況が、完全長およびequal-lengthの逆方向反復配列のアノテーションの頻度に大きく影響することが明らかになった。

 NCBI Nucleotideに登録されているプラスチドゲノムの数は近年飛躍的に増加しており、今後10年間でさらに多くのゲノムが登録されると考えられる。 airpgを用いることで、これらのプラスチドゲノムの逆方向反復配列やその配列アノテーションに自動的にアクセスして評価することができ、公開されているプラスチドゲノムの信頼性向上に貢献することができる。このソフトウェアは、http://pypi.python.org/pypi/airpgPythonパッケージインデックスから自由に利用できる。

 

インストール

ubuntu18.04のpython3.9の仮想環境にpipで導入した(pip install git+https://github.com/michaelgruenstaeudl/airpg.git)。

Github

#pip(pypi)
pip install airpg

pip install git+https://github.com/michaelgruenstaeudl/airpg.git

> airpg_analyze.py -h

$ airpg_analyze.py -h

usage: airpg_analyze.py [-h] --infn INFN --outfn OUTFN --mail MAIL [--blocklist BLOCKLIST] [--query QUERY] [--recordsdir RECORDSDIR] [--datadir DATADIR] [--verbose]

 

Michael Gruenstaeudl <m.gruenstaeudl@fu-berlin.de>, Tilman Mehl <tilmanmehl@zedat.fu-berlin.de> -- Copyright (C) 2019-2021 Michael Gruenstaeudl and Tilman Mehl -- Retrieve the plastid genomes identified by the first

script and evaluate their inverted repeats -- 2021.03.05

 

optional arguments:

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

  --infn INFN, -i INFN  Path to input file; input is a summary table of NCBI accessions (tab-delimited, accession numbers in second column)

  --outfn OUTFN, -o OUTFN

                        Path to output file that contains information on IR positions and length

  --mail MAIL, -m MAIL  Mail address needed for Entrez search on NCBI PubMed (any valid mail address works)

  --blocklist BLOCKLIST, -b BLOCKLIST

                        (Optional) Path to blocklist file

  --query QUERY, -q QUERY

                        (Optional) Entrez string to query NCBI PubMed

  --recordsdir RECORDSDIR, -r RECORDSDIR

                        (Optional) Path to records directory

  --datadir DATADIR, -d DATADIR

                        (Optional) Path to data directory

  --verbose, -v         (Optional) Enable verbose logging

>  airpg_confirm.py -h

$ airpg_confirm.py -h

usage: airpg_confirm.py [-h] --infn INFN --outfn OUTFN --datadir DATADIR [--minlength MINLENGTH] [--maxlength MAXLENGTH] [--verbose]

 

Michael Gruenstaeudl <m.gruenstaeudl@fu-berlin.de>, Tilman Mehl <tilmanmehl@zedat.fu-berlin.de> -- Copyright (C) 2019-2021 Michael Gruenstaeudl and Tilman Mehl -- Retrieve the plastid genomes identified by the first

script and evaluate their inverted repeats -- 2021.05.19

 

optional arguments:

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

  --infn INFN, -i INFN  Path to input file; input is a summary table on reported IR positions and length (tab-delimited, accession numbers in first column)

  --outfn OUTFN, -o OUTFN

                        Path to output file that contains extended table IR positions and length

  --datadir DATADIR, -d DATADIR

                        Path to folder containing record-specific subfolders that store each record's complete sequence in FASTA format

  --minlength MINLENGTH, -n MINLENGTH

                        (Optional) Minimal length of IR for BLASTing

  --maxlength MAXLENGTH, -x MAXLENGTH

                        (Optional) Maximum length of IR for BLASTing

  --verbose, -v         (Optional) Enable verbose logging

> airpg_identify.py -h

$ airpg_identify.py  -h

usage: airpg_identify.py [-h] -o OUTFN [-q QUERY] [-b BLOCKLIST] [-u]

 

Michael Gruenstaeudl <m.gruenstaeudl@fu-berlin.de>, Tilman Mehl <tilmanmehl@zedat.fu-berlin.de> -- Copyright (C) 2019-2020 Michael Gruenstaeudl and Tilman Mehl -- Conduct a query of NCBI Nucleotide and identify plastid

genome records stored there -- 2020.12.17

 

optional arguments:

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

  -o OUTFN, --outfn OUTFN

                        Path to output file

  -q QUERY, --query QUERY

                        (Optional) Entrez query that will replace the standard query

  -b BLOCKLIST, --blocklist BLOCKLIST

                        (Optional) Path to file of blocklisted genera that will be removed from the retrieved plastid sequences

  -u, --update_only     (Optional) Only add entries with more recent creation date than the most recent existing entry

エラーが出たので対処方法を確認中。

=> 著者に修正していただいて動作するようになった (現在はGithubの方で対応)。

 

 

4つチュートリアルが用意されています。レポジトリにアクセスしてみて下さい。

引用

airpg: automatically accessing the inverted repeats of archived plastid genomes
Tilman Mehl & Michael Gruenstaeudl 
BMC Bioinformatics volume 22, Article number: 413 (2021) 

 

MMseqs2 コマンド其の4 分類群をアサインする mmseqs taxonomyコマンド

 

今年出た論文(*1)より

MMseqs2 taxonomyは、メタゲノムのコンティグに分類学上のラベルを付与する新しいツールである。各コンティグから可能性のある全てのタンパク質断片を抽出し、分類学的なアノテーションに貢献できるものを素早く取り出し、それらにロバストなラベルを付与し、weighted votingによってコンティグの分類学的な同一性を決定する。そのフラグメント抽出ステップは、生命の全てのドメインの分析に適している。MMseqs2 taxonomyは、最先端のツールと比較して2~18倍の速度で動作し、分類学的参照データベースの作成と操作、分類学的割り当ての報告と視覚化のための新しいモジュールも含まれている。
MMseqs2 taxonomyは、フリーのMMseqs2オープンソースソフトウェアパッケージの一部であり、LinuxmacOSWindowsで利用できる。

 

user guide

https://mmseqs.com/latest/userguide.pdf

tutorial

https://github.com/soedinglab/MMseqs2/wiki/Tutorials

 

 

f:id:kazumaxneo:20210916002624p:plain

MMseqs2 taxnomyのLCA (wiki) 推定アプローチ (2bLCA)。マニュアルより転載。

 

 

 インストール

以前の記事を参照

> mmseqs taxonomy -h

# mmseqs taxonomy -h

usage: mmseqs taxonomy <i:queryDB> <i:targetDB> <o:taxaDB> <tmpDir> [options]

 By Milot Mirdita <milot@mirdita.de> & Martin Steinegger <martin.steinegger@snu.ac.kr> & Eli Levy Karin <eli.levy.karin@gmail.com>

options: prefilter:                      

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

 --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 [nucl:nucleotide.out,aa:VTML80.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]

 --k-score INT                    k-mer threshold for generating similar k-mer lists [2147483647]

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

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

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

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

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

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

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

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

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

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

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

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

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

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

align:                          

 -a BOOL                          Add backtrace string (convert to alignments with mmseqs convertalis module) [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 [1]

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

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

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

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

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

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

 --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 FLOAT                      Pseudo count admixture strength [1.000]

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

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

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

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

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

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

 --orf-filter INT                 Prefilter query ORFs with non-selective search

                                  Only used during nucleotide-vs-protein classification

                                  NOTE: Consider disabling when classifying short reads [1]

 --orf-filter-e DOUBLE            E-value threshold used for query ORF prefiltering [1.000E+02]

 --orf-filter-s FLOAT             Sensitivity used for query ORF prefiltering [2.000]

 --lca-mode INT                   LCA Mode 1: single search LCA , 2/3: approximate 2bLCA, 4: top hit [3]

 --tax-output-mode INT            0: output LCA, 1: output alignment 2: output both [0]

 --majority FLOAT                 minimal fraction of agreement among taxonomically assigned sequences of a set [0.500]

 --vote-mode INT                  Mode of assigning weights to compute majority. 0: uniform, 1: minus log E-value, 2: score [1]

 --lca-ranks STR                  Add column with specified ranks (',' separated)

 --tax-lineage INT                0: don't show, 1: add all lineage names, 2: add all lineage taxids [0]

 --blacklist STR                  Comma separated list of ignored taxa in LCA computation [12908:unclassified sequences,28384:other sequences]

 --rescore-mode INT               Rescore diagonals with:

                                  0: Hamming distance

                                  1: local alignment (score only)

                                  2: local alignment

                                  3: global alignment

                                  4: longest alignment fulfilling window quality criterion [0]

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

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

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

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

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

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

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

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

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

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

                                  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]

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

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

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

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

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

 --start-sens FLOAT               Start sensitivity [4.000]

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

common:                         

 --compressed INT                 Write compressed output [0]

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

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

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

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

 

examples:

 # Download a sequence database with taxonomy information

 mmseqs databases UniProtKB/Swiss-Prot swissprotDB tmp

 

 # Assign taxonomy based on 2bLCA

 mmseqs taxonomy queryDB swissprotDB result tmp

 

 # Assign taxonomy based on top hit

 mmseqs taxonomy queryDB swissprotDB result tmp --lca-mode 4

 

 # Assign taxonomy without ORF prefilter

 # Classifies higher percentage for short nucleotide input (e.g. short reads) at the cost of speed

 mmseqs taxonomy queryNuclDB swissprotDB result tmp --orf-filter 0

 

 # Create a Krona report

 mmseqs taxonomyreport swissprotDB result report.html --report-mode 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)

 - Mirdita M, Steinegger M, Breitwieser F, Soding J, Levy Karin E: Fast and sensitive taxonomic assignment to metagenomic contigs. bioRxiv, 2020.11.27.401018 (2020)

 

 

 

データベースの準備

前回の記事に書きましたが、mmseqs databaseコマンドを使うことでビルド済みのデータベースをダウンロードすることができます。ここではメタゲノムアセンブリのコンティグにラベルを当てることを想定して、Unirefデータベースをダウンロードします(紹介した時の記事)。

mkdir tmp
mmseqs databases UniRef90 UniRef90_DB tmp

(UnirefはseqTaxDBに対応している)

 

mmseqs taxonomyのラン

1、クエリの配列(ここではメタゲノムのアセンブリ)をデータベースに変換する。

mmseqs createdb metagenome.fasta queryDB

 

2、一時ディレクトリtmpの作成(高いI/Oが必要なので大規模な検索ではSSDなどを使う。ただし、tmpには十分な空き容量がないといけない)。

 

3、 mmseqs taxonomyの実行

1で作成したクエリ配列のデータベース、準備したデータベース、出力prefix、作業ディレクトリの順番に指定する。

mmseqs taxonomy queryDB UniRef90_DB output_prefix tmp --majority 0.5 --lca-mode 3 --vote-mode 1 --tax-lineage 2 --orf-filter 1
  • --tax-lineage INT   0: don't show, 1: add all lineage names, 2: add all lineage taxids [0]
  •  --majority FLOAT   minimal fraction of agreement among taxonomically assigned sequences of a set [0.500]
  • --vote-mode INT   Mode of assigning weights to compute majority. 0: uniform, 1: minus log E-value, 2: score [1]
  • --lca-mode INT   LCA Mode 1: single search LCA , 2/3: approximate 2bLCA, 4: top hit [3]
  • --orf-filter INT   Prefilter query ORFs with non-selective search. Only used during nucleotide-vs-protein classification.  NOTE: Consider disabling when classifying short reads [1]

6フレームで含まれるタンパク質断片を抽出し、2bLCA法で検索する(--lca-mode 3)。次に、割り当てられたすべてのフラグメントの中から投票を行い、-log(E-value)の重みが50%以上の支持を得た(--majority 0.5)最も具体的な分類ラベルを選択する (--vote-mode 1)。パラメータ--tax-lineage 2は、NCBIのtaxidsとしての完全な系統情報を出力することを示している。

ラン中はデフォルトで利用可能なすべてのCPUが使用される。その間、他の計算は出来ないほど効率的にコアが使用されるので注意する。計算が終わると、output_prefixのファイルが複数生成される。

 

4、出力をTSVに変換する。

1で作成したクエリ配列のデータベース、3の出力prefix、出力タブ区切りファイルの順番に指定する。

mmseqs createtsv queryDB output_prefix result.tsv

f:id:kazumaxneo:20210915214301p:plain

 

easy-taxonomyコマンドを使うとより簡単に実行できます。

mmseqs easy-taxonomy QUERY.fasta targetDB output_prefix  tmp

 

引用

Fast and sensitive taxonomic assignment to metagenomic contigs 
M Mirdita, M Steinegger, F Breitwieser, J Söding, E Levy Karin
Bioinformatics, Published: 18 March 2021 

 

関連