Understanding the significance of BLAT score calculations for read fragment alignment

Question

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 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?

Answer 1

It would appear that BLAT and exonerate (as suggested by @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, 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, 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 2

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>