How to blast UniProt/RefSeq database from local clusters?

Question

I want to identify all members of several protein families and study the evolutionary relationships between them.

I tried to do exhaustive search for proteins homologous to a sequence via blast (against UniProt or RefSeq database). I figured out three ways to do that:

1. Webpage blast tools. But they have restrictions on parameters like memory usage and number of hits retrieved.
2. Download all proteomes and makeblastdb locally, as demonstrated here: https://dmnfarrell.github.io/bioinformatics/local-refseq-db. But downloading all fasta files can be slow for me.
3. Use the webservice clients provided by EMBL (https://www.ebi.ac.uk/jdispatcher/docs/). It provides two tools for blast, but I failed to find a thorough documentation for their usage or Q&A about them (is it a popular tool, or out-dated?). Direct troubleshooting by looking into the code may be feasible but overwhelming to me.

I wonder whether you have better ways to do that or suggestions/resources for above plans.

---

• "Why BLAST?" -- I thought algorithms like PSI-BLAST may help find remote homology, which can be un-annotated or annotated as other domains. In other words, my initial wish is to find something (maybe) evolutionarily related to a domain based on sequence/site/structure conservation but are not similar enough to be classified as "same domain".

• "How many is several?" -- I do case studies. Around ten domains.

• "given that BLAST doesn't / can't show all matches?" -- It's a problem I didn't expect. Local blast+ software seems only return 500 subjects. I thought it was unlimited.

Answer 1

The standard tool for this is NCBI's blast+ package. It is installable on most systems, and should be in the repositories of your distribution if you're using Linux.

Once installed, you can use it to run blast queries against all of the NCBI databases, including RefSeq. For example, using this human tP53 protein sequence saved as p53.pep:

>NP_000537.3 cellular tumor antigen p53 isoform a [Homo sapiens]
MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGPDEAPRMPEAA
PPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAKSVTCTYSPALNKMFCQLAKT
CPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHERCSDSDGLAPPQHLIRVEGNLRVEYLDDRN
TFRHSVVVPYEPPEVGSDCTTIHYNYMCNSSCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGR
DRRTEEENLRKKGEPHHELPPGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALEL
KDAQAGKEPGGSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD

You can blast it against the translated nr database with:

blastp -query p53.pep -db nr  -remote > blast.out

The output consists of several header lines with the blast references, and query details, then a list of matches:

uery= NP_000537.3 cellular tumor antigen p53 isoform a [Homo sapiens]

Length=393

RID: 4H2ECMAT013
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

XP_063556659.1 cellular tumor antigen p53 isoform X1 [Gorilla gor...  814     0.0   
7XZZ_K Chain K, Cellular tumor antigen p53 [Homo sapiens]             814     0.0   
8R1F_C Chain C, Cellular tumor antigen p53 [Homo sapiens]             814     0.0   
NP_000537.3 cellular tumor antigen p53 isoform a [Homo sapiens]       813     0.0   
AYE20617.1 mutant tumor protein p53 [Homo sapiens]                    813     0.0   
AYE20613.1 mutant tumor protein p53 [Homo sapiens]                    812     0.0   
AYE20623.1 mutant tumor protein p53 [Homo sapiens]                    812     0.0   
AAA61212.1 p53 antigen [Homo sapiens]                                 812     0.0   
XP_004058559.3 cellular tumor antigen p53 isoform X2 [Gorilla gor...  812     0.0   
CAA42633.1 p53 transformation suppressor [Homo sapiens]               812     0.0   
EAW90142.1 tumor protein p53 (Li-Fraumeni syndrome), isoform CRA_...  812     0.0   
AFN61606.1 tumor suppressor p53 [Homo sapiens]                        812     0.0   
UNB14053.1 p53 [synthetic construct]                                  811     0.0   
WOX59026.1 mutant tumor protein p53 transcript variant 1 [Homo sa...  811     0.0   
AYE20618.1 mutant tumor protein p53 [Homo sapiens]                    811     0.0
[...]   

Followed by the high-scoring segment pairs (HSPs):

>XP_063556659.1 cellular tumor antigen p53 isoform X1 [Gorilla gorilla gorilla]
Length=408

 Score = 814 bits (2103),  Expect = 0.0, Method: Compositional matrix adjust.
 Identities = 392/393 (99%), Positives = 392/393 (99%), Gaps = 0/393 (0%)

Query  1    MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP  60
            MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP
Sbjct  16   MEEPQSDPSVEPPLSQETFSDLWKLLPENNVLSPLPSQAMDDLMLSPDDIEQWFTEDPGP  75

Query  61   DEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAK  120
            DEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAK
Sbjct  76   DEAPRMPEAAPPVAPAPAAPTPAAPAPAPSWPLSSSVPSQKTYQGSYGFRLGFLHSGTAK  135

Query  121  SVTCTYSPALNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHE  180
            SVTCTYSP LNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHE
Sbjct  136  SVTCTYSPTLNKMFCQLAKTCPVQLWVDSTPPPGTRVRAMAIYKQSQHMTEVVRRCPHHE  195

Query  181  RCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNS  240
            RCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNS
Sbjct  196  RCSDSDGLAPPQHLIRVEGNLRVEYLDDRNTFRHSVVVPYEPPEVGSDCTTIHYNYMCNS  255

Query  241  SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELP  300
            SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELP
Sbjct  256  SCMGGMNRRPILTIITLEDSSGNLLGRNSFEVRVCACPGRDRRTEEENLRKKGEPHHELP  315

Query  301  PGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPG  360
            PGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPG
Sbjct  316  PGSTKRALPNNTSSSPQPKKKPLDGEYFTLQIRGRERFEMFRELNEALELKDAQAGKEPG  375

Query  361  GSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD  393
            GSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD
Sbjct  376  GSRAHSSHLKSKKGQSTSRHKKLMFKTEGPDSD  408

Now, depending on what you want to do, you are likely to want to use tblastn or perhaps psi-blast instead. Protein homology is about far more than simple sequence similarity so you should use tools that take structure and/or protein domains into account. HMMR is another good option. But that's a whole different question and we would require the exact biological question you are trying to answer to be able to guide you there. This tool, however, is the most common and easiest way I know to blast against the various standard DBs from your local machine.

Answer 2

If an intermediate requirement is "do thousands of BLAST searches", the question is almost always an XY problem, i.e. asking about the attempted intermediate solution rather than your actual problem. That appears to be the case here.

So I'm not going to answer your question about exhaustive BLAST searches, because you have not convinced me that it would help your actual problem. I'm also not sure what you mean by "exhaustive search for proteins", because the thresholds for identity depends on the protein family and how deep you want to go back into the phylogenetic history - some proteins within a family would still be considered related at 60% identity (where precise sequences are less important), and others would be considered different at 95% identity (e.g. essential housekeeping proteins). There's no universal rule or threshold, which makes an exhaustive search a lesson in futility.

This is why we ask for context. For a specific protein family, there may be an appropriate tool, but such a tool is unlikely to exist in the general case.

With respect to BLAST, it is a basic local alignment and search tool, and best used for an initial idea of the origin of individual sequences. It's a low-throughput tool that searches a comprehensive database of reference sequences.

I use it for sanity spot checks when doing real-time sequencing, to make sure that the reads the sequencer is producing are similar to what I expect. There have been a few times when I've stopped a sequencing run quickly because I've noticed from the BLAST results that the output is nonsense, or there is too much contamination, or the thing I think I'm sequencing is something different from what is actually being sequenced.

Some things to remember about BLAST (and its related family of BLAST-like tools):

• It is basic - BLAST makes assumptions about the nature of the sequence data and the desired output. In an attempt to be fast, it will output the first hits it finds that satisfy the specified parameters, not the best hits.

• It carries out a local alignment - for long, chimeric, or otherwise disordered inputs, BLAST will output good substring matches, up until the match starts being less good. It will not attempt to align entire sequences. If you have a protein with multiple subunits, you might get a different impression of phylogeny by these local matches when compared to a global whole-sequence alignment algorithm.

• It is slow - due to the size of the database and the way it carries out its searches, BLAST works well enough for short individual sequences, but has less useful processing time scales when doing hundreds to thousands of searches. There are other similar tools that are faster (e.g. Diamond), but they're less accessible. NCBI has done an excellent job at making their web BLAST tool easy to use, but painful enough for computationally taxing tasks that most people don't do that.

I may add more to this explanation later, but that'll do for now.