[Bioperl-l] standalone blast outformat

[dimitark at bii.a-star.edu.sg]

Hi Florent and Jason,
Thank you for your replies.

well i spent quite some time to figure that out. But still no results.
Yes i use BLAST+. So I tried BLASTN as it is, without bioperl. Just to  
see if i can do: blastn -outfmt "7 qcovs other_things". Well it works  
like this. But in bioperl can not pass do that. I dont know why tho.

After i tried to understand how tile_hsps() works. But no results.  
Somehow the documentation on that matter is not clear and not a simple  
example is present. So i could not figure out how to use it.

Now i am trying to figure out how to use FASTA or SSEARCH to solve my  
problem. I suppose i have to get the sequences of the HSPs  
representing my hit and then run ssearch like this:

    ssearch my_query_seq.fa my_hit_hsps.fa

Then parse the output. Just have to figure out what exactly do i need  
from the output :)

Thank you again
D.

Quoting Jason Stajich <jason.stajich at gmail.com>:

> Assuming 1 HSP you can do:
>
> my $qcov = $hsp->query->length / $query->length;
>
> Otherwise you have to add up the HSPs and their coverage (and remove  
> redundancy of overlapping hits) to determine what parts of the query  
> were aligned across all the HSPs that were found. You can try the  
> tile_hsps function to also get this. Or if you are using WUBLAST you  
> can get the -links option to grab the logical groups of  
> non-overlapping HSPs to tile across.
>
> Or use FASTA or SSEARCH which will provide a single longest path  
> through the alignment.
>
>
> On Oct 29, 2013, at 10:12 PM, Florent Angly <florent.angly at gmail.com> wrote:
>
>> Hi Dimitar,
>>
>> I suppose you are using BLAST+, not legacy BLAST. There is a HOWTO
>> here that can get you started:
>> http://www.bioperl.org/wiki/HOWTO:BlastPlus
>> Regarding your issue, I had a look at the blastn manual, and I doesn't
>> look like you can write '-outfmt "7 qcovs ..."' to me. If anything, it
>> would be '-outfmt "qcovs ..."'.
>> So, maybe try to get your command to work as intended in a terminal
>> before putting it in BioPerl.
>> Best,
>>
>> Florent
>>
>>
>> On Wed, Oct 30, 2013 at 1:45 PM,  <dimitark at bii.a-star.edu.sg> wrote:
>>> Hi guys,
>>> i was wondering how can i get the query coverage when blasting  
>>> with Bioperl.
>>> So i figure out one way is to specify the outformat to tabular(7) and the
>>> specify that i want the qcovs(query coverage per subject). From  
>>> BLAST manual
>>> i see one can do: -outfmt "7 qcovs ..." but in Bioperl i can only say:
>>> -outformat => 7. If i do: outformat => "7 qcovs" i get an error.
>>>
>>> Well another option i suppose is to get the lengths of all hsps from a hit
>>> and estimate the query coverage that way but im not exactly sure if this is
>>> correct way.
>>>
>>> Any idea how can i get the query coverage?
>>>
>>> Thank you very much for your time and any hints
>>>
>>> Regards
>>> Dimitar
>>>
>>> _______________________________________________
>>> Bioperl-l mailing list
>>> Bioperl-l at lists.open-bio.org
>>> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>> _______________________________________________
>> Bioperl-l mailing list
>> Bioperl-l at lists.open-bio.org
>> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>
> Jason Stajich
> jason.stajich at gmail.com
> jason at bioperl.org