6
$\begingroup$

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

$\endgroup$
  • $\begingroup$ If I'm not wrong, qcovhsp is a perc_identity for ONE HSP. I would have to calculate manually for each hit as an average for all HSPs perc_identity for one hit. $\endgroup$ – biotech Nov 12 '14 at 13:53
  • $\begingroup$ I even don't like hit_perc_identity because of its reductive information assumption. Average is normally influenced by extreme values. $\endgroup$ – biotech Nov 12 '14 at 13:55
  • $\begingroup$ What species are you comparing? How distant are they? If close enough, "real" homologs should be able to form a single HSP spanning most of the target sequences. $\endgroup$ – terdon Nov 12 '14 at 14:06
  • $\begingroup$ You can calculate the average HSP coverage for the entire sequence (from all alignments). You can also calculate average score per residue. The other option is to do an end-to-end alignment with tools such as stretcher. $\endgroup$ – WYSIWYG Nov 12 '14 at 14:19
  • $\begingroup$ Check out my new question related to BLAST: biology.stackexchange.com/questions/23958/… $\endgroup$ – biotech Nov 12 '14 at 18:53
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ They are strains of the same bacterial specie. I agree with the thresholds you suggest. Thanks $\endgroup$ – biotech Nov 12 '14 at 18:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.