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:

      <Hsp_midline>I A+G GL+ PAG E+IDAGGLLVMPGGID H</Hsp_midline>

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?

  • 1
    $\begingroup$ 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? $\endgroup$
    – holmrenser
    May 18 '17 at 9:31
  • $\begingroup$ @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... $\endgroup$ May 18 '17 at 9:37
  • $\begingroup$ 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. $\endgroup$
    – terdon
    May 18 '17 at 10:49
  • $\begingroup$ @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. $\endgroup$ May 18 '17 at 11:08
  • 1
    $\begingroup$ 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 :). $\endgroup$
    – terdon
    May 18 '17 at 12:09

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!


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_def>43989.cce_0262 (Cyanothece ATCC 51142)</Hit_def>
      <Hsp_midline>+ GA+C M+++E+ L KL++ A+KLIE LLGP+</Hsp_midline>

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>

  • $\begingroup$ An HSP is still allowed to contain gaps, though. So I'd still have to parse one of the qseq or hseq :( $\endgroup$ May 18 '17 at 12:06
  • 1
    $\begingroup$ 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. $\endgroup$
    – holmrenser
    May 18 '17 at 13:26
  • $\begingroup$ 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). $\endgroup$ May 18 '17 at 13:32
  • $\begingroup$ 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. $\endgroup$
    – holmrenser
    May 18 '17 at 13:35
  • $\begingroup$ 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). $\endgroup$ May 18 '17 at 13:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.