how to identify blast hits with only 1 HSP (not limit the number of HSP)

[Kennels]

Hi,

Is there a way to identify blast hits where there was only one HSP to the database?

i.e. I don't want to limit the number of HSP's displayed, but I want to identify all queries which only had one HSP to the database to begin with in the blast result.

Maybe there is an accessory script somewhere that parses a blast report and identify the queries with this property?

thanks

[whataBamBam]

If I understand what your after correctly then I do the same using this..

cat BlastOutputFile | grep -B 3 -A 1 "# 1 hits found"

When BlastOutputFile is in outformat 7. Obviously if your using Xml output or any other format it'll be different.

[Kennels]

If I understand what your after correctly then I do the same using this..

cat BlastOutputFile | grep -B 3 -A 1 "# 1 hits found"

When BlastOutputFile is in outformat 7. Obviously if your using Xml output or any other format it'll be different.

Thanks bambam, but that's not what I'm after.

Yes that will pull out those with 1 match, but what I'm after is those that have only have one alignment 'section' (i.e. high-scoring segment pair).

e.g. an isoform could match to only one target in the database, but half of it could be aligning to to the 5' end, then a gap, then the other half the 3' end of the target.

I'm after hits which aligns without any interruptions, i.e. only 1 HSP. Blastn options don't seem to do what I want (e.g. -ungapped, -culling_limit, -max_hsps_per_subject just mask the information).

[sphil]

I'm after hits which aligns without any interruptions, i.e. only 1 HSP. Blastn options don't seem to do what I want (e.g. -ungapped, -culling_limit, -max_hsps_per_subject just mask the information).

So you mean you want to parse for full-length hsp?

[Kennels]

Yes and no.
Also looking for partial alignment, but only those with single hsp.

[atcghelix]

I usually use BioPerl for something like this. The following will get you close (the formatting of the output blast file generated by Bio::SearchIO::Writer::TextResultWriter is slightly different, but basically the same).

Usage would be:
perl script.pl --blast inputBlastFile.txt --out reducedBlastFile.txt

#!/usr/bin/perl

use strict;
use warnings;
use Getopt::Long;
use Bio::SearchIO;
use Bio::SearchIO::Writer::TextResultWriter;

my $blastFile;
my $outFile;

GetOptions("blast=s" => \$blastFile,
           "out=s"   => \$outFile) || die "Couldn't get options with GetOpt::Long.\n";

my $blastIO = Bio::SearchIO->new(-file => $blastFile,
                                 -format => 'blast'); #you might need to change this

my $searchIOWriter = Bio::SearchIO::Writer::TextResultWriter->new();
my $out = Bio::SearchIO->new(-writer => $searchIOWriter,
                             -file => '>outtest.txt');

while (my $result = $blastIO->next_result) {
    if ($result->num_hits() == 1) {
        my $hit = $result->next_hit();
        if ($hit->num_hsps() == 1) {
            print "A single HSP for a single hit found for " . $result->query_name() . "\n";
            $out->write_result($result);
        }
    }
}

[Kennels]

Thanks I'll give it a try.

[lindenb]

if you used the XML output of blast, you could use the following XSLT stylesheet:

<?xml version='1.0'  encoding="ISO-8859-1"?>
<xsl:stylesheet xmlns:xsl='http://www.w3.org/1999/XSL/Transform' version='1.0' >
<xsl:output method="text"/>
<xsl:template match="/">
<xsl:apply-templates select="BlastOutput/BlastOutput_iterations/Iteration/Iteration_hits/Hit"/>
</xsl:template>
<xsl:template match="Hit">
<xsl:if test="count(Hit_hsps/Hsp) = 1">
<xsl:value-of select="Hit_id"/><xsl:text>	</xsl:text>
<xsl:value-of select="Hit_def"/><xsl:text>	</xsl:text>
<xsl:value-of select="Hit_accession"/><xsl:text>
</xsl:text>
</xsl:if>
</xsl:template>
</xsl:stylesheet>

usage:

$ xsltproc --novalid  stylesheet.xsl  blast.xml

[whataBamBam]

Thanks bambam, but that's not what I'm after.

Yes that will pull out those with 1 match, but what I'm after is those that have only have one alignment 'section' (i.e. high-scoring segment pair).

e.g. an isoform could match to only one target in the database, but half of it could be aligning to to the 5' end, then a gap, then the other half the 3' end of the target.

I'm after hits which aligns without any interruptions, i.e. only 1 HSP. Blastn options don't seem to do what I want (e.g. -ungapped, -culling_limit, -max_hsps_per_subject just mask the information).

Ah sorry I get you. Yes that would be really useful. I'll have a look at the script posted and see what that's doing.

I think Gmap will do something like that.. it isn't really designed to match pairs of transcript length sequences (it's designed to map a transcript length sequence to a genome) but I actually do use it to map transcripts to ESTs and it seems to work really well. I then filter the output based on length of the alignment and can find the matches where the EST aligned along it's whole length without any significant gaps. Which is probably quite similar to what your trying to do.

Possibly I really shouldn't be using Gmap to do this so if someone can tell me that then great because I actually do it all the time.

[Kennels]

I usually use BioPerl for something like this. The following will get you close (the formatting of the output blast file generated by Bio::SearchIO::Writer::TextResultWriter is slightly different, but basically the same).

Just letting you know the script did what I wanted.
I just modified the script a little:

#!/usr/bin/perl

use strict;
use warnings;
use Getopt::Long;
use Bio::SearchIO;
use Bio::SearchIO::Writer::TextResultWriter;

my $blastFile = $ARGV[0];     ## added one argument
my $outFile;

GetOptions("blast=s" => \$blastFile,
           "out=s"   => \$outFile) || die "Couldn't get options with GetOpt::Long.\n";

my $blastIO = Bio::SearchIO->new(-file => $blastFile,
                                 -format => 'blasttable'); ### changed from 'blast'

my $searchIOWriter = Bio::SearchIO::Writer::TextResultWriter->new();
my $out = Bio::SearchIO->new(-writer => $searchIOWriter,
                             -file => '>outtest.txt');

while (my $result = $blastIO->next_result) {
    if ($result->num_hits() == 1) {
        my $hit = $result->next_hit();
        if ($hit->num_hsps() == 1) {
            
	print $result->query_name() . "\n";   ## just print the result ID
            $out->write_result($result);
        }
    }
}

Since my output was in -outfmt 6 from blastn, I changed the format to 'blasttable'
Also I just wanted the ID of the query for single HSP, so I just modified the stdout and redirect to a file to get my list of IDs.
I ran as:

parse_blastreport_single_hsp.pl Result.blastn6 > singleHSPs.txt

For others, also need to make sure the output has '-max_target_seqs' in blastn is set to more than 1.

Many thanks!