Find protein homologs with BLASTp

Question

I'm trying to find homologs of a set of proteins using BLASTp. I'm working with custom databases.

I'm using evalue of 0.00001 as threshold.

I would like to filter queries having hits with >90% identities. Because BLASTp output is based on HSPs, I cannot filter by %identities/query, only by HSP.

I would like to know how to do this and also if I'm following a reasonable strategy.

Here is an alignment example: qcovs=100 but qcovhsp lower.

qseqid           sseqid         pident length mismatch gapopen qstart   qend  sstart    send   evalue bitscore qcovs qcovhsp
HPNK_01698      HAPS_0519       81.88   596     75      5       630     1225    615     1177    0.0      889    100     49
HPNK_01698      HAPS_0519       49.17   301     115     8       84      366     201     481     2e-56    214    100     23
HPNK_01698      HAPS_0519       53.64   261     61      6       436     684     616     828     6e-49    191    100     20
HPNK_01698      HAPS_0519       46.61   251     62      3       332     510     584     834     6e-46    181    100     15
HPNK_01698      HAPS_0519       53.27   214     79      4       1       194     1       213     1e-45    180    100     16
HPNK_01698      HAPS_0519       55.96   218     60      8       550     764     643     827     1e-40    164    100     18
HPNK_01698      HAPS_0519       51.56   225     61      7       516     731     642     827     1e-38    157    100     18
HPNK_01698      HAPS_0519       49.57   230     77      6       484     713     643     833     1e-37    154    100     19
HPNK_01698      HAPS_0519       57.89   76      26      1       364     433     760     835     1e-13   76.3    100     6

Code used

Make database

makeblastdb -in $Hparasuisfastadatabase -out H_parasuis_strains_gb_ALL.fna_databaseBLAST -dbtype prot -parse_seqids 

Run BLAST

blastp -db H_parasuis_strains_gb_ALL.fna_databaseBLAST -query 'out_2.fasta' -out HPNK_selected_vs_H_parasuis_strainss.tblastn -evalue 0.00001 -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore qcovs qcovhsp" -max_target_seqs 50

Thanks, Bernardo

Answer 1

First of all, if you want 90% identity, you can discard this hit. None of the HSPs pass that threshold. What's more, since you're working with proteins, there are no splicing issues involved and you should be able to get a single HSP spanning most of the query and subject sequences. Assuming, of course, you have a true homolog.

In your output, I see many small, overlapping HSPs, most of which have low identity. I can't be sure without seeing the sequence but it's a safe bet that what you have there are low complexity/repetitive regions and that's why you have so many separate HSPs. The only half way decent one starts at position 630 of the query sequence and is only 595 residues long, less than half your query protein. Either you have a very diverged N-terminal region or your HSP is just a conserved domain. Again, I would need to see the actual sequence alignment to be sure but that doesn't look like a true homolog (assuming your species are reasonably close which they must be if you're using a 90% identity threshold).

So, always assuming that your species are close enough to expect decent homologs, I would simply ignore the shorter HSPs and deal with those that represent more than, say, 80% of my query's length at >=90% identity. Shorter hits will most often be conserved domains or repetitive/low complexity regions. The thresholds you choose depend on the species you are studying.

If your species are not that close, don't use BLASTP at all. Instead, you can use something like hmmer. Collect a set of homologs from various species for each of your query proteins, build a matrix using these and use that matrix to search your database. You could also use Selenoprofiles which uses a similar approach.