Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-18T15:44:11.847Z
Answers.LastActivityDate: 2017-05-18T15:50:52.953Z
Answers.LastEditDate: 2017-05-18T15:50:52.953Z
Answers.UserDisplayName: Sam Nicholls
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-18T12:06:48.567Z
Answers.Text: An HSP is still allowed to contain gaps, though. So I'd still have to parse one of the `qseq` or `hseq` :(
Answers.UserDisplayName: Sam Nicholls
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-18T13:26:07.920Z
Answers.Text: 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_.
Answers.UserDisplayName: holmrenser
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-18T13:32:44.263Z
Answers.Text: 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).
Answers.UserDisplayName: Sam Nicholls
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-18T13:35:31.233Z
Answers.Text: 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.
Answers.UserDisplayName: holmrenser
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-18T13:38:37.557Z
Answers.Text: 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).
Answers.UserDisplayName: Sam Nicholls
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-05-18T10:05:47.703Z
Answers.LastActivityDate: 2017-05-18T10:05:47.703Z
Answers.UserDisplayName: holmrenser
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-18T09:31:27.710Z
Comments.Text: 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?
Comments.UserDisplayName: holmrenser
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-18T09:37:19.093Z
Comments.Text: @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...
Comments.UserDisplayName: Sam Nicholls
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-18T10:49:46.660Z
Comments.Text: 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.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-18T11:08:39.523Z
Comments.Text: @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).
Comments.UserDisplayName: Sam Nicholls
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-18T11:17:31.967Z
Comments.Text: 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.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-18T12:05:53.493Z
Comments.Text: @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!
Comments.UserDisplayName: Sam Nicholls
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-18T12:09:11.110Z
Comments.Text: 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 :).
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-05-18T12:11:54.580Z
Comments.Text: @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?
Comments.UserDisplayName: Sam Nicholls
Tags: <phylogenetics><blat><blast>
Title: Understanding the significance of BLAT score calculations for read fragment alignment
ThreadText: # 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:
  
  ```
    [...]
    <Hit_hsps>
      <Hsp>
        <Hsp_num>1</Hsp_num>
        <Hsp_bit-score>51.6</Hsp_bit-score>
        <Hsp_score>122</Hsp_score>
        <Hsp_evalue>9.8e-06</Hsp_evalue>
        <Hsp_query-from>2</Hsp_query-from>
        <Hsp_query-to>98</Hsp_query-to>
        <Hsp_hit-from>35</Hsp_hit-from>
        <Hsp_hit-to>65</Hsp_hit-to>
        <Hsp_query-frame>2</Hsp_query-frame>
        <Hsp_hit-frame>0</Hsp_hit-frame>
        <Hsp_identity>24</Hsp_identity>
        <Hsp_positive>27</Hsp_positive>
        <Hsp_gaps>1</Hsp_gaps>
        <Hsp_align-len>32</Hsp_align-len>
           <Hsp_qseq>ITAIGAGLQGPAGCEVIDAGGLLVMPGGIDTH</Hsp_qseq>
           <Hsp_hseq>IAAVGTGLE-PAGAEIIDAGGLLVMPGGIDVH</Hsp_hseq>
        <Hsp_midline>I A+G GL+ PAG E+IDAGGLLVMPGGID H</Hsp_midline>
      </Hsp>
    </Hit_hsps>
  </Hit>
  
  ```
  
  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      ++      <PROTEIN>   662     321     354     <READ>      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](http://genome.ucsc.edu/FAQ/FAQblat.html#blat4) 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**: @username_1 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**: @username_2 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**: @username_2 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**: @username_2 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 [@username_2](https://bioinformatics.stackexchange.com/users/298/username_2)) 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:
  
  ```
  <Hit>
    <Hit_num>1</Hit_num>
    <Hit_id>gnl|BL_ORD_ID|1515029</Hit_id>
    <Hit_def>43989.cce_0262 (Cyanothece ATCC 51142)</Hit_def>
    <Hit_accession>1515029</Hit_accession>
    <Hit_len>65</Hit_len>
    <Hit_hsps>
      <Hsp>
        <Hsp_num>1</Hsp_num>
        <Hsp_bit-score>40.0466</Hsp_bit-score>
        <Hsp_score>92</Hsp_score>
        <Hsp_evalue>0.00664016</Hsp_evalue>
        <Hsp_query-from>155</Hsp_query-from>
        <Hsp_query-to>253</Hsp_query-to>
        <Hsp_hit-from>12</Hsp_hit-from>
        <Hsp_hit-to>44</Hsp_hit-to>
        <Hsp_query-frame>-1</Hsp_query-frame>
        <Hsp_hit-frame>0</Hsp_hit-frame>
        <Hsp_identity>17</Hsp_identity>
        <Hsp_positive>27</Hsp_positive>
        <Hsp_gaps>0</Hsp_gaps>
        <Hsp_align-len>33</Hsp_align-len>
        <Hsp_qseq>LRGAICSMEHIEEALGKLKDWARKLIELLLGPR</Hsp_qseq>
        <Hsp_hseq>ITGAVCLMDYLEKVLEKLRELAQKLIETLLGPQ</Hsp_hseq>
        <Hsp_midline>+ GA+C M+++E+ L KL++ A+KLIE LLGP+</Hsp_midline>
      </Hsp>
    </Hit_hsps>
  </Hit>
  
  ```
  
  You should be able to get the position information from the following tags: `<Hsp_query-from>`,`<Hsp_query-to>`,`<Hsp_hit-from>`,`<Hsp_hit-to>`.
  
  To translate back to you original DNA sequence I can imagine you need the `<Hsp_query-frame>` tag.
  
  And then of course you get the evalue for free: `<Hsp_evalue>`
  
  > **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).  
  
  ---
  
  