<p>It would appear that <code>BLAT</code> and <a href="http://www.ebi.ac.uk/about/vertebrate-genomics/software/exonerate" rel="nofollow noreferrer"><code>exonerate</code></a> (as suggested by <a href="https://bioinformatics.stackexchange.com/users/298/terdon">@terdon</a>) are the only pieces of software that provide enough information for "HSPFragments" to be parsed from their outputs, at least according to <a href="http://biopython.org/DIST/docs/api/Bio.SearchIO._model-module.html" rel="nofollow noreferrer">BioPython's documentation</a>, anyway:</p>

<blockquote>
  <p>Most search programs only have HSPs with one HSPFragment in them,
  making these two objects inseparable. However, there are programs
  (e.g. BLAT and Exonerate) which may have more than one HSPFragment
  objects in any given HSP. If you are not using these programs, you can
  safely consider HSP and HSPFragment as a single union.</p>
</blockquote>

<p>It would seem the options are sticking with <code>BLAT</code> and generating scores afterwards with the <a href="http://genome.ucsc.edu/FAQ/FAQblat.html#blat4" rel="nofollow noreferrer">script described in their FAQ</a>, or using <code>exonerate</code> which provides a raw score as part of its default output. What these scores actually mean I will have to investigate, but for now I think the case is closed!</p>

CC BY-SA 3.0
2017-05-18T15:44:11.847Z
2017-05-18T15:50:52.953Z
2017-05-18T15:50:52.953Z
Sam Nicholls
<p>I would suggest using blastx and parsing the output XML. As an example, have a look at this: <a href="http://www.perlmonks.org/?node_id=1006367" rel="nofollow noreferrer">http://www.perlmonks.org/?node_id=1006367</a>.</p>

<p>The relevant part of the hit would be the following:</p>

<pre><code>&lt;Hit&gt;
  &lt;Hit_num&gt;1&lt;/Hit_num&gt;
  &lt;Hit_id&gt;gnl|BL_ORD_ID|1515029&lt;/Hit_id&gt;
  &lt;Hit_def&gt;43989.cce_0262 (Cyanothece ATCC 51142)&lt;/Hit_def&gt;
  &lt;Hit_accession&gt;1515029&lt;/Hit_accession&gt;
  &lt;Hit_len&gt;65&lt;/Hit_len&gt;
  &lt;Hit_hsps&gt;
    &lt;Hsp&gt;
      &lt;Hsp_num&gt;1&lt;/Hsp_num&gt;
      &lt;Hsp_bit-score&gt;40.0466&lt;/Hsp_bit-score&gt;
      &lt;Hsp_score&gt;92&lt;/Hsp_score&gt;
      &lt;Hsp_evalue&gt;0.00664016&lt;/Hsp_evalue&gt;
      &lt;Hsp_query-from&gt;155&lt;/Hsp_query-from&gt;
      &lt;Hsp_query-to&gt;253&lt;/Hsp_query-to&gt;
      &lt;Hsp_hit-from&gt;12&lt;/Hsp_hit-from&gt;
      &lt;Hsp_hit-to&gt;44&lt;/Hsp_hit-to&gt;
      &lt;Hsp_query-frame&gt;-1&lt;/Hsp_query-frame&gt;
      &lt;Hsp_hit-frame&gt;0&lt;/Hsp_hit-frame&gt;
      &lt;Hsp_identity&gt;17&lt;/Hsp_identity&gt;
      &lt;Hsp_positive&gt;27&lt;/Hsp_positive&gt;
      &lt;Hsp_gaps&gt;0&lt;/Hsp_gaps&gt;
      &lt;Hsp_align-len&gt;33&lt;/Hsp_align-len&gt;
      &lt;Hsp_qseq&gt;LRGAICSMEHIEEALGKLKDWARKLIELLLGPR&lt;/Hsp_qseq&gt;
      &lt;Hsp_hseq&gt;ITGAVCLMDYLEKVLEKLRELAQKLIETLLGPQ&lt;/Hsp_hseq&gt;
      &lt;Hsp_midline&gt;+ GA+C M+++E+ L KL++ A+KLIE LLGP+&lt;/Hsp_midline&gt;
    &lt;/Hsp&gt;
  &lt;/Hit_hsps&gt;
&lt;/Hit&gt;
</code></pre>

<p>You should be able to get the position information from the following tags: <code>&lt;Hsp_query-from&gt;</code>,<code>&lt;Hsp_query-to&gt;</code>,<code>&lt;Hsp_hit-from&gt;</code>,<code>&lt;Hsp_hit-to&gt;</code>.</p>

<p>To translate back to you original DNA sequence I can imagine you need the <code>&lt;Hsp_query-frame&gt;</code> tag.</p>

<p>And then of course you get the evalue for free: <code>&lt;Hsp_evalue&gt;</code></p>

CC BY-SA 3.0
2017-05-18T12:06:48.567Z
An HSP is still allowed to contain gaps, though. So I'd still have to parse one of the `qseq` or `hseq` :(
Sam Nicholls
CC BY-SA 3.0
2017-05-18T13:26:07.920Z
Correct me if I'm wrong, but the tags I mentioned (Hsp_query-from, Hsp_query-to, Hsp_hit-from, Hsp_hit-to) are coordinates on _the original sequence_.
holmrenser
CC BY-SA 3.0
2017-05-18T13:32:44.263Z
You are correct (I believe)! But the HSP can contain gaps, so I can't just take those co-ordinates as given to extract the DNA from the read between those positions (as I don't want the gapped nucleotides).
Sam Nicholls
CC BY-SA 3.0
2017-05-18T13:35:31.233Z
The fact that the read aligns to that position, and that the coordinates are on the original sequence means you actually can. Whether or not you want to do a gapped alignment is a different matter I would think.
holmrenser
CC BY-SA 3.0
2017-05-18T13:38:37.557Z
Not for what I need - long story short is I'm trying to "flatten" the reads into a pseudo-consensus. I need to know exactly where (what actual position) the gaps are, so I can handle them appropriately. I can't just use the starts and ends of the read as given, because taking that subsequence will include one (or more) nucleotides that should be removed to accommodate the gap(s).
Sam Nicholls
CC BY-SA 3.0
2017-05-18T10:05:47.703Z
2017-05-18T10:05:47.703Z
holmrenser
<p>I'm attempting to reconstruct the read fragments that form a hit to a target, excluding/skipping gaps. <em>BLAST</em>'s <code>outfmt 6</code>, and alternatives such as Diamond's <code>M8</code> format only give the number of gap opens on the alignment. One can't just use the <code>start</code> and <code>end</code> positions of the alignment as indices to slice the read DNA sequence, as this will include the nucleotides that were gapped by the alignment.</p>

<p><code>PSL</code> appeared to be the only format I could find that gave the actual start positions (and lengths) of each of the pieces that formed a hit (that is, each gap forces a new block).</p>

<p>The regular <em>BLAST</em> output is human readable and is a nightmare to parse. The <code>XML</code> format mentioned has a structure for each HSP, but an HSP can still include gaps:</p>

<pre><code>  [...]
  &lt;Hit_hsps&gt;
    &lt;Hsp&gt;
      &lt;Hsp_num&gt;1&lt;/Hsp_num&gt;
      &lt;Hsp_bit-score&gt;51.6&lt;/Hsp_bit-score&gt;
      &lt;Hsp_score&gt;122&lt;/Hsp_score&gt;
      &lt;Hsp_evalue&gt;9.8e-06&lt;/Hsp_evalue&gt;
      &lt;Hsp_query-from&gt;2&lt;/Hsp_query-from&gt;
      &lt;Hsp_query-to&gt;98&lt;/Hsp_query-to&gt;
      &lt;Hsp_hit-from&gt;35&lt;/Hsp_hit-from&gt;
      &lt;Hsp_hit-to&gt;65&lt;/Hsp_hit-to&gt;
      &lt;Hsp_query-frame&gt;2&lt;/Hsp_query-frame&gt;
      &lt;Hsp_hit-frame&gt;0&lt;/Hsp_hit-frame&gt;
      &lt;Hsp_identity&gt;24&lt;/Hsp_identity&gt;
      &lt;Hsp_positive&gt;27&lt;/Hsp_positive&gt;
      &lt;Hsp_gaps&gt;1&lt;/Hsp_gaps&gt;
      &lt;Hsp_align-len&gt;32&lt;/Hsp_align-len&gt;
         &lt;Hsp_qseq&gt;ITAIGAGLQGPAGCEVIDAGGLLVMPGGIDTH&lt;/Hsp_qseq&gt;
         &lt;Hsp_hseq&gt;IAAVGTGLE-PAGAEIIDAGGLLVMPGGIDVH&lt;/Hsp_hseq&gt;
      &lt;Hsp_midline&gt;I A+G GL+ PAG E+IDAGGLLVMPGGID H&lt;/Hsp_midline&gt;
    &lt;/Hsp&gt;
  &lt;/Hit_hsps&gt;
&lt;/Hit&gt;
</code></pre>

<p>I'm trying to be lazy, I don't want to have to parse the <code>Hsp_qseq</code> and <code>Hsp_hseq</code> to determine exactly where my gaps are. I just want co-ordinates of where my read hits a protein. The <code>PSL</code> format (below) tells me this (<code>tStarts</code>):</p>

<pre><code>match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block   blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
[...]
21      4       0       0       1       8       1       24      ++      &lt;PROTEIN&gt;   662     321     354     &lt;READ&gt;      101     2       101     2       8,17,   321,337,        2,50,
</code></pre>

<p>As an alternative, I've executed the <code>BLAT</code> aligner as a means to find alignments of translated DNA reads against a small protein database. I selected <code>BLAT</code> as it provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing). That is, the PSL has a comma delimited field that enumerates each of the starting positions of a query to a target.</p>

<p>However, <code>BLAT</code> outputs a significantly greater number of alignments than alternative tools. This is likely due to its identity percentage threshold defaulting to 25% when emulating a <code>blastx</code> style query.
I can cope with this (apologies to my disks for the I/O), but for filtering purposes, the output <code>PSL</code> format appears to leave me nothing but raw match counts to work with.</p>

<p>The <code>PSL</code> gives me what I want, except now I don't have information on its significance. I know that the <a href="http://genome.ucsc.edu/FAQ/FAQblat.html#blat4" rel="nofollow noreferrer">BLAT FAQ</a> describes a Perl script for adding "web-based identity and score calculations". I can actually get identity percentages and scores from this script, but I cannot tell what the score actually represents. Is it comparable to bit-score? How can I recover something like an e-value from a <code>PSL</code> output file?</p>

CC BY-SA 3.0
2017-05-18T09:31:27.710Z
Just out of curiosity: can you elaborate on how BLAT 'provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing)'? What about blast HSPs?
holmrenser
CC BY-SA 3.0
2017-05-18T09:37:19.093Z
@holmrenser Don't you only get the HSPs by parsing the XML output? This is some laze on my part, as the HSPs in my XML described the hits from the perspective of the translated DNA, but I want to know where the alignments on the original sequence are without too much of my own error-prone interference. BLAT just gives the offset. It might be that I've overlooked something in the HSP data though...
Sam Nicholls
CC BY-SA 3.0
2017-05-18T10:49:46.660Z
No! Not at all! The default output format for blast is a simple text file and there are flags to get one HSP per line (no sequence, only the scores etc). You *can* get XML as well but that's just a way of making your life harder IMO. Please [edit] your question and clarify what the original (scientific) question you are asking is. Are you looking for exact matches? Distant homologs? Finally, it would be better if you could include an example of your PSL file if that's really what you want parsed.
terdon
CC BY-SA 3.0
2017-05-18T11:08:39.523Z
@terdon Despite my dislike of XML parsing, it seems a more structured way to get out the components of a hit compared to trying to infer structure from a sequence of lines. For what it's worth `BLAT` [doesn't have flags for outputting an HSP per line](https://genome.ucsc.edu/goldenpath/help/blatSpec.html#blatUsage).
Sam Nicholls
CC BY-SA 3.0
2017-05-18T11:17:31.967Z
BLAT might not, I admit I've only ever used BLAT online for simple things. BLAST, however (both NCBI and wublast), certainly does. Look at the "tabular" option for the `-outfm` flag in NCBI blast for example. I've spent more time than I care to remember parsing blast output, if you clarify the exact question you have, I should be able to help.
terdon
CC BY-SA 3.0
2017-05-18T12:05:53.493Z
@terdon Thanks, I know what options `BLAST` has, but I'm asking about something else. I have updated my question to try and clarify though!
Sam Nicholls
CC BY-SA 3.0
2017-05-18T12:09:11.110Z
Fair enough (and sorry, didn't mean to appear pushy). I was just wondering if using another tool (like BLAST which would give you the e-values in a parseable table format) wouldn't be a simpler approach. I realize that "redo the analysis with this tool instead" is not a perfect solution though :).
terdon
CC BY-SA 3.0
2017-05-18T12:11:54.580Z
@terdon Hah, no worries. A lot of bioinformatics is "damn, I need to run this with another tool". I'd happily move this to chat if you have any ideas. For what it's worth, `BLAST` is too slow (for the size of the data) and still doesn't give the granularity on gap information that I'm looking for :( Any thoughts?
Sam Nicholls
<phylogenetics><blat><blast>
Understanding the significance of BLAT score calculations for read fragment alignment
# Understanding the significance of BLAT score calculations for read fragment alignment

(Asked by: username_0 on 2017-05-18)

I'm attempting to reconstruct the read fragments that form a hit to a target, excluding/skipping gaps. BLAST's outfmt 6, and alternatives such as Diamond's M8 format only give the number of gap opens on the alignment. One can't just use the start and end positions of the alignment as indices to slice the read DNA sequence, as this will include the nucleotides that were gapped by the alignment.
PSL appeared to be the only format I could find that gave the actual start positions (and lengths) of each of the pieces that formed a hit (that is, each gap forces a new block).
The regular BLAST output is human readable and is a nightmare to parse. The XML format mentioned has a structure for each HSP, but an HSP can still include gaps:
  [...]
  
    
      1
      51.6
      122
      9.8e-06
      2
      98
      35
      65
      2
      0
      24
      27
      1
      32
         ITAIGAGLQGPAGCEVIDAGGLLVMPGGIDTH
         IAAVGTGLE-PAGAEIIDAGGLLVMPGGIDVH
      I A+G GL+ PAG E+IDAGGLLVMPGGID H
    
  


I'm trying to be lazy, I don't want to have to parse the Hsp_qseq and Hsp_hseq to determine exactly where my gaps are. I just want co-ordinates of where my read hits a protein. The PSL format (below) tells me this (tStarts):
match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block   blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
[...]
21      4       0       0       1       8       1       24      ++         662     321     354           101     2       101     2       8,17,   321,337,        2,50,

As an alternative, I've executed the BLAT aligner as a means to find alignments of translated DNA reads against a small protein database. I selected BLAT as it provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing). That is, the PSL has a comma delimited field that enumerates each of the starting positions of a query to a target.
However, BLAT outputs a significantly greater number of alignments than alternative tools. This is likely due to its identity percentage threshold defaulting to 25% when emulating a blastx style query.
I can cope with this (apologies to my disks for the I/O), but for filtering purposes, the output PSL format appears to leave me nothing but raw match counts to work with.
The PSL gives me what I want, except now I don't have information on its significance. I know that the BLAT FAQ describes a Perl script for adding "web-based identity and score calculations". I can actually get identity percentages and scores from this script, but I cannot tell what the score actually represents. Is it comparable to bit-score? How can I recover something like an e-value from a PSL output file?

**Tags:** phylogenetics, blat, blast
> **Comments:**
> **username_1**: Just out of curiosity: can you elaborate on how BLAT 'provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing)'? What about blast HSPs?  
> **username_0**: @holmrenser Don't you only get the HSPs by parsing the XML output? This is some laze on my part, as the HSPs in my XML described the hits from the perspective of the translated DNA, but I want to know where the alignments on the original sequence are without too much of my own error-prone interference. BLAT just gives the offset. It might be that I've overlooked something in the HSP data though...  
> **username_2**: No! Not at all! The default output format for blast is a simple text file and there are flags to get one HSP per line (no sequence, only the scores etc). You *can* get XML as well but that's just a way of making your life harder IMO. Please [edit] your question and clarify what the original (scientific) question you are asking is. Are you looking for exact matches? Distant homologs? Finally, it would be better if you could include an example of your PSL file if that's really what you want parsed.  
> **username_0**: @terdon Despite my dislike of XML parsing, it seems a more structured way to get out the components of a hit compared to trying to infer structure from a sequence of lines. For what it's worth `BLAT` [doesn't have flags for outputting an HSP per line](https://genome.ucsc.edu/goldenpath/help/blatSpec.html#blatUsage).  
> **username_2**: BLAT might not, I admit I've only ever used BLAT online for simple things. BLAST, however (both NCBI and wublast), certainly does. Look at the "tabular" option for the `-outfm` flag in NCBI blast for example. I've spent more time than I care to remember parsing blast output, if you clarify the exact question you have, I should be able to help.  
> **username_0**: @terdon Thanks, I know what options `BLAST` has, but I'm asking about something else. I have updated my question to try and clarify though!  
> **username_2**: Fair enough (and sorry, didn't mean to appear pushy). I was just wondering if using another tool (like BLAST which would give you the e-values in a parseable table format) wouldn't be a simpler approach. I realize that "redo the analysis with this tool instead" is not a perfect solution though :).  
> **username_0**: @terdon Hah, no worries. A lot of bioinformatics is "damn, I need to run this with another tool". I'd happily move this to chat if you have any ideas. For what it's worth, `BLAST` is too slow (for the size of the data) and still doesn't give the granularity on gap information that I'm looking for :( Any thoughts?  

---

## 2 answers

### Answer by username_0 on 2017-05-18. Score: 2 | Accepted answer

It would appear that `BLAT` and [`exonerate`](http://www.ebi.ac.uk/about/vertebrate-genomics/software/exonerate) (as suggested by [@terdon](https://bioinformatics.stackexchange.com/users/298/terdon)) are the only pieces of software that provide enough information for "HSPFragments" to be parsed from their outputs, at least according to [BioPython's documentation](http://biopython.org/DIST/docs/api/Bio.SearchIO._model-module.html), anyway:

> Most search programs only have HSPs with one HSPFragment in them,
> making these two objects inseparable. However, there are programs
> (e.g. BLAT and Exonerate) which may have more than one HSPFragment
> objects in any given HSP. If you are not using these programs, you can
> safely consider HSP and HSPFragment as a single union.

It would seem the options are sticking with `BLAT` and generating scores afterwards with the [script described in their FAQ](http://genome.ucsc.edu/FAQ/FAQblat.html#blat4), or using `exonerate` which provides a raw score as part of its default output. What these scores actually mean I will have to investigate, but for now I think the case is closed!


---

### Answer by username_1 on 2017-05-18. Score: 1

I would suggest using blastx and parsing the output XML. As an example, have a look at this: http://www.perlmonks.org/?node_id=1006367.
The relevant part of the hit would be the following:

  1
  gnl|BL_ORD_ID|1515029
  43989.cce_0262 (Cyanothece ATCC 51142)
  1515029
  65
  
    
      1
      40.0466
      92
      0.00664016
      155
      253
      12
      44
      -1
      0
      17
      27
      0
      33
      LRGAICSMEHIEEALGKLKDWARKLIELLLGPR
      ITGAVCLMDYLEKVLEKLRELAQKLIETLLGPQ
      + GA+C M+++E+ L KL++ A+KLIE LLGP+
    
  


You should be able to get the position information from the following tags: ,,,.
To translate back to you original DNA sequence I can imagine you need the  tag.
And then of course you get the evalue for free:

> **Comments:**
> **username_0**: An HSP is still allowed to contain gaps, though. So I'd still have to parse one of the `qseq` or `hseq` :(  
> **username_1**: Correct me if I'm wrong, but the tags I mentioned (Hsp_query-from, Hsp_query-to, Hsp_hit-from, Hsp_hit-to) are coordinates on _the original sequence_.  
> **username_0**: You are correct (I believe)! But the HSP can contain gaps, so I can't just take those co-ordinates as given to extract the DNA from the read between those positions (as I don't want the gapped nucleotides).  
> **username_1**: The fact that the read aligns to that position, and that the coordinates are on the original sequence means you actually can. Whether or not you want to do a gapped alignment is a different matter I would think.  
> **username_0**: Not for what I need - long story short is I'm trying to "flatten" the reads into a pseudo-consensus. I need to know exactly where (what actual position) the gaps are, so I can handle them appropriately. I can't just use the starts and ends of the read as given, because taking that subsequence will include one (or more) nucleotides that should be removed to accommodate the gap(s).  

---

