macでインフォマティクス

macでインフォマティクス

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

UniProtKBデータベースをダウンロードしてBLAST検索する。

2022/07/10誤字修正、07/12誤字修正

2022/07/28 ツイート追記

2024/10/09 追記

 

Universal Protein Resource (UniProt)は、European Bioinformatics Institute (EBI) (*2)とSIB Swiss Institute Bioinformaticsが共同研究して構築している知識ベースである(*1)。
タンパク質の機能、分類、相互参照など、広範囲に渡って可能な限り多くのアノテーション情報がキュレーションされて付与されており、正確で一貫性のある豊富なアノテーションが付けられたタンパク質情報の知識ベースとなっている。UniProtは異なる用途に最適化された3つのデータベースから構成されている(*1, 7)。

  1. UniProt Reference Clusters (UniRef) は、近縁の配列を1つのレコードにまとめ、配列の類似性検索を高速化したものである(*4)(紹介)。3つの解像度で配列空間を完全にカバーすることにより、ユーザーの目的に合わせてより高速な類似性検索を実現している(*5)。
  2. UniProt Archive (UniParc) は、全タンパク質配列の包括的なリポジトリである。配列の一意な識別子と配列のみからなる(UPIに続く10個の16進数)。UniParcのタンパク質はデータベースの相互参照により、ソースデータベースとリンクされている。UniParcは、テキストベースと配列ベースの検索が可能である。UniParcに対して配列類似性検索を行うことは、UniParcにクロスリファレンスされている全てのデータベースに対して同じ検索を行うことと同じである(*6)。
  3. UniProt Knowledgebase (UniProtKB) は、2つの部分から構成されている(*7)。片方は完全に手動でアノテーションされたSwiss-Protセクション(UniProtKB/Swiss-Prot)で、専門家が文献|サブミットから手作業で抽出して収録している。もう片方はUniProtKB/Swiss-Proを補完するためにコンピュータアノテーションによって行われるTrEMBLセクション(UniProtKB/TrEMBL)で、EMBL/GenBank/DDBJ Nucleotide Sequence Databasesに存在する全てのコーディング配列(CDS)の翻訳と、文献から抽出したタンパク質配列、UniProtKB/Swiss-Protに提出されたタンパク質配列を収録している (*9)。Swiss-ProtよりもTrEMBLのエントリー数の方が圧倒的に多いが、TrEMBLデータベースのタンパク質は推定に過ぎず、アノテーションが不十分である。完全手動のアノテーション待ちのデータベースになっている(*7)。PIR-PSDリリースは2004年にこれらのセクションに完全に統合された。

UniProtを構成するSwiss-ProtとTrEMBLは元々別のプロジェクトで(*1, 4)、Swiss-ProtはSwiss Institute Bioinformaticsの知識ベースだった(*1, 7, 8)、TrEMBLはゲノムデータの爆発的増加に対応して1996年に導入された (*5, 7)。

UniProt Knowledgebaseは、長期休暇以外、8週間ごとに更新が追加される(*10)。また、UniProtKBのデータは一般公開日よりもかなり前に「凍結」されるため、関連データベースの更新とUniProtKBの一般公開版に関連データが現れるまでに8~16週間程度の遅れが生じる可能性がある(*10)。過去のリリースは少なくとも2年間はftpサイトに保存される。それ以降は、その年の最初のリリースであるYYYY_01アーカイブしている。そのため、あるリリースやサブセットを学習したツールを開発する場合、その年の最初のリリースを使用することが推奨されている(*10)。

 

 

UniProtKBのダウンロードとBLASTサーチ

UniProtKBのデータは、小さなデータであればこのリンク先から直接ダウンロードできますが、大きなデータはFTPサーバからダウンロードすることが推奨されています(* 11)。ヨーロッパ、中東、アフリカに住んでいる場合、イギリスまたはスイスにあるミラーサイトがより高速とされます(リンク)。

Downloads

https://www.uniprot.org/help/downloads

Current release (ftp.uniprot.org)

https://ftp.uniprot.org/pub/databases/uniprot/current_release/

(READMEに各ファイルの詳細が書かれている)

 

ここでは、Andrew Perryさんが公開されているスクリプトを使ってuniprot_sprotとuniprot_tremblをダウンロードし、makeblastdbでBLASTデータベースを作るまでを自動で実行します。更新が早いので自動化するスクリプトは役立ちます。このスクリプトは、uniprot_sprotとuniprot_tremblそれぞれのBLASTデータベースを作り、並行して、uniprot_sprotとuniprot_tremblの配列を結合したファイル(リリース年+その年のリリース番号がファイル名)からもBLASTデータベースを作ります。つまり、合計3つのデータベースを作る。ただ、uniprot_tremblは肥大しており重くなっています。冗長なヒットも多数出るので、このような方法でUniProtKBのデータベースをローカルに構築するよりも、クラスター化されたUniRefのデータベースの利用が一般的です。その点は注意して下さい(研究目的による)。また、検索を高速化するには、BLASTPではなく、mmseqs2紹介)やDIAMOND紹介)の利用も考えて下さい (DIAMONDがBLASTのデータベースを使えるにする取り組みも進行中)。また、ダウンロードせずともUniprotKBやUniRefを使えるオンラインサービスも存在します(紹介)。DIAMONDのデータベースを作る例(link)。

 

Github
pansapiens/download-uniprot-release.sh

Download the latest Uniprot (Swissprot+Trembl) release, make a BLAST database · GitHub

#ディレクトリを作って実行。ローカルにデータベースを作る。
mkdir uniprotDB && cd uniprotDB/
./download-uniprot-release.sh

出力

2022_02/  (古い順)

上の方にあるのがダウンロードしたfasta形式のファイル。2022/07/11アクセス時は、uniprot_sprot.fasta(完全なUniProtKB/Swiss-Protデータセット)には567,483配列、uniprot_trembl.fasta(完全なUniProtKB/TrEMBLデータセット)には231,354,261配列含まれていました。2007年当時TrEMBLは4,672,908エントリーだったと記載(*14)があるので、15年で爆発的に増えたことが分かります。 ディレクトリ2022_02/のサイズは505GBほどでした。

 

 

BLASTPを実行します。ここではuniprot_sprotとuniprot_tremblを結合したデータベースを使いますが、Swiss-Protの配列のみ使いたいなら、uniprot_sprotを使います。-outfmt 7ならスコア順に並びます。(追記;最大ヒット数を増やす-max_target_seqsがデフォルト500なので必要なら増やす(multi-fastaなら各配列ごとの最大ヒット))

#uniprot_sprot & uniprot_trembl(2022-02は更新バージョンによって変わる)
blastp -db 2022_02/uniprot_2022_02 -query protein.faa -out result.txt -evalue 1e-1 -num_threads 20 -outfmt 7 -max_target_seqs 100000

#uniprot_sprot
blastp -db 2022_02/uniprot_sprot -query protein.faa -out result.txt -evalue 1e-1 -num_threads 20 -outfmt 7 -max_target_seqs 100000

$top hitを取り出すなら
echo qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore |sed s/' '/$'\t'/g > best_hits.txt
#then, extract best hit
cat result.txt | awk '/hits found/{getline;print}' | grep -v "#" >> best_hits.txt

クエリ数が多いとBLASTサーチには時間がかかります。1クエリだけの検索では、20スレッド指定で1分かかりました(配列とパラメータによる)。

 

BLAST結果からタンパク質情報を取り出す

ヒットした一意な識別子は、UniProt Retrieve/ID mappingサービスで、機能的なアノテーションや分類群などの情報を得たり(ID mappingと呼ぶ)、該当する識別子の配列をFASTA形式でダウンロードすることができます。また、NCBI refseq protein IDなどへの変換も可能です。Retrieve/ID mappingサービスは最近インターフェイスが更新されました。

補足;UniProt Retrieve/ID mappingサービスは、クエリが冗長な場合、冗長なIDは除いて結果を出すので注意する。例えば、blastp でヒットした10000配列のIDをクエリとしてRetrieve/ID mappingでタンパク質情報を取り出したが、2000ヒットしかなかったなら、10000配列に冗長なIDが多数あったと推定される。

 

コメント

ツイッターで教えてもらいましたが、最近NCBIもClustered nr(ClusteredNR)という90%の同一性と90%の長さの範囲内でクラスタ化されたnrを公開していますね(現在Experimental)(*15)。90%の同一性と90%の長さということは、少しパラメータが異なりますがUniRef90 に近い圧縮率でしょうか。この機能が構築された動機は、BLAST検索結果がほぼ同一の、”ゲノムのみが出た”研究がされていない生物からの過剰なヒットでBLAST結果が埋め尽くされる事が多いのが一因かと思います(refseq representative genomesを選ぶと減らすことはできる)。ある距離範囲で代表配列(よく注釈されたもの)だけ残して他の冗長な配列は捨てることで、より均一な分類学的距離でサンプリングされた配列セットを作れます。そうやって作られたデータベースを使うと、同一のタンパク質が複数の種に存在するとか、汚染であるとかなども調べやすくなります。詳細はNCBIの記事を読んでみて下さい(*15)。

興味深いのは、MMseqs2ソフトウェアが使われているという点でしょうか。実際Clustered nrを試してみましたが、リソースの割り当てが少ないのか現在は検索時間が相当長めでした。結果も微妙に気になりましたが、実験中のサービスなのでまずは試してみて下さい。

 

補足;DeepMind社がEMBL-EBIと提携し、AlphaFoldタンパク質構造データベースを公開しています(*16)。データベースの最初のリリースでは、ヒトプロテオームに含まれる2万個のタンパク質全てと、いくつかのモデル生物のプロテオームが含まれています(長いタンパク質はスプリットされている)。2022/07/11現在のv2アップデートでは、48のモデル生物とSwiss-Protの542,380配列まで拡張されています(リンク)(*A)

今後、早期にUniRef90データベースの1億個のタンパク質の大部分をカバーするようにデータベースを拡張する予定とアナウンスされています( FTPリンク)。2022年7月現在、UniRef100は306,293,517配列、UniRef90は147,407,377配列、UniRef50は52,523,202配列含まれています。1億ということは、UniRefにもある短い配列を除いて、文字通りUniRef90のほぼ全ての配列の推定立体構造情報を作ろうとしているのだと思われます。これを企業主導で行っているのは違和感がありますが、SWISS-PROTも最初は個人の努力で進んでいたわけで(*8)、今後研究者に重要性が認知されて重要性が増せば、自然に公的なサービスに変わっていく感じでしょうか。

ちなみに、クエリ配列もColabFold (論文)で立体構造予測し、AlphaFold DBをダウンロードすれば、ローカルで類似構造予測ができます。その場合、ローカルにインストールしたFoldseek (preprint)やDALI(*B)が使えます。Swiss-Protとの類似構造予測結果が、アミノ酸配列レベルの検索やサポートベクターマシンの分類器(論文)などの結果と違いが出るのか、twilight zone(論文)が変わるのかなど気になるところです。

*A 上のuniprot_sprot.fastaは567,483配列あるので少しだけ少ない。Swiss-Protには最短2a.aと非常に短いアミノ酸配列も含まれておりそれが除かれたと思われる。Swiss-Protになぜこのような配列があるのかよく分からない。

*B DALIオンライン版でもalphafoldDBを対象に検索できるものの、alphafoldDBのSwiss-Protは選べない。

 

追記

こちらを使うと、今回行った、データベースのダウンロード、diamond blastx 、ID mappingの工程を全て自動で実行できます。

 

引用

wiki

*1 Swiss-Prot - Wikipedia

*2 欧州バイオインフォマティクス研究所 - Wikipedia

*3 ExPASy - Wikipedia

*4 DDBJ サービスの UniProt への対応

https://www.ddbj.nig.ac.jp/news/ja/2005-03-24.html

*5 UniRef

https://www.uniprot.org/help/uniref

*6 UniParc

https://www.uniprot.org/help/uniparc

*7 uniprot current release README (FTP)

https://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/README

*8 SWISS-PROTは いかにして生まれたか?

https://www.jstage.jst.go.jp/article/kagakutoseibutsu1962/39/4/39_4_239/_pdf/-char/ja

Uniprot Q&A

*9 Why is UniProtKB composed of 2 sections, UniProtKB/Swiss-Prot and UniProtKB/TrEMBL? 

https://www.uniprot.org/help/uniprotkb_sections 

*10 How frequently is UniProt released?

https://www.uniprot.org/help/synchronization

*11 Downloaded data seems incomplete or corrupted - how can I get help with download problems?

https://www.uniprot.org/help/metalink

*12 Mapping between UniProtKB and NCBI resources (GeneID, RefSeq): how does it work?

https://www.uniprot.org/help/ncbi_mappings

*13 ExPASy

https://www.expasy.org/

SIB Swiss Institute of Bioinformaticsが運営するオンライン・バイオインフォマティクス・リソース。160以上のデータベースとソフトウェアツールへのアクセスを提供し、ゲノム、プロテオミクス、構造生物学から進化・系統、システム生物学、医化学まで、生命科学と臨床研究の幅広い分野をサポートする、統合的なポータルサイト。ExPASyは1993年8月1日に設置され、ライフサイエンス分野の最初のウェブサイトであり、世界で最初の150番目のウェブサイトであった(*3, 6)。Expasyの現行バージョンは、2020年10月にリリースされた(wikiより)。

*14 UniProtKB/TrEMBL

http://www.bioinfo.pte.hu/more/TrEMBL.htm

*15

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

*16

https://www.deepmind.com/open-source/alphafold-protein-structure-database

*17

Swiss-Protには手動でキュレーションされた配列が含まれているが、含まれている配列の中にはソースとして微妙なものもあるように見える。例えばQ87JM4はRND型輸送体のMacBをコードしているが、ソースをみるとこの配列が決定された菌のゲノムプロジェクトとなっている。従っておそらく実験的に証明されたものではなく、配列同一性や保存されたモチーフの共通性からMacBとみなされたのだと想像されるが、この菌はもう1つMacBを持っていて、もしろそちらのMacBの方が近縁種ではよく保存されている。そうするとこのタンパク質は本当にMacBなのか、それとも輸送体スーパーファミリーとして起源は同じでも異なる基質を運ぶように進化したものではないかということも想像されてくる。しかし、多剤輸送体ならアノテーションとしては間違いはないとも言える(引用はほかにされていない)。「Swiss-Protの配列なので機能は100%完正しいだろう」とは限らないかもしれない事は少しだけ頭の中に入れておいたほうがよい。また時間が許せばアノテーションソースも少し見る癖をつけておくのが望ましい。

 

参考

*12について

UniProt Retrieve/ID mapping serviceで変換できる。

Uniprot channnel

Inside expert biocuration in UniProtKB/Swiss-Prot

 

関連

オンラインならSANSparallelが利用できます。


2022/07/28

思ったより早く来ました。想定より数が多いのも興味深いところですね。