Blast Output Results Analysis

[milo0615]

Hello,

I have successfully blasted 2500 merged COGS against 40 different assembly databases. However, I now have 40 large blast results which I need to analyze and select the assembly with the most hits. My questions are the following:

- What would be the best way or the best practice to analyze all of the blast results to check for the assembly with the most hits?

- Is there a free application that would help with the analysis?

Thank you in advance.

[atcghelix]

When you say you want the assembly with the most hits, do you mean the assembly that had the fewest number of "no hits found" (i.e. perfect score is 2500--each had at least one hit), or the assembly with the most hits (each of the 2500 COGS have multiple hits---perfect score would be way more than 2500 total hits).

[milo0615]

Originally Posted by atcghelix

When you say you want the assembly with the most hits, do you mean the assembly that had the fewest number of "no hits found" (i.e. perfect score is 2500--each had at least one hit), or the assembly with the most hits (each of the 2500 COGS have multiple hits---perfect score would be way more than 2500 total hits).

I would say the assembly that had the fewest number of "no hits found." Do you think that would be a better selection?

[atcghelix]

I'm not sure--it sort of depends on what you want to know. Short kmer values will probably have more hits overall, but the hits will be shorter. I often am trying to find the assembly kmer value that has the highest number of hits that fully span the length of the query sequence.

If you just want to see how many 'No hits found' there are, you can use:
grep -c 'No hits found' <filename>

[milo0615]

Originally Posted by atcghelix

I'm not sure--it sort of depends on what you want to know. Short kmer values will probably have more hits overall, but the hits will be shorter. I often am trying to find the assembly kmer value that has the highest number of hits that fully span the length of the query sequence.

If you just want to see how many 'No hits found' there are, you can use:
grep -c 'No hits found' <filename>

So you find the best kmer assembly based on the "highest number of hits?" I just want to know which is the best optimal kmer assembly by blasting it against COGS, or should I pick the best hits based on the e-value? How do you pick your best kmer?

[atcghelix]

I'm not sure what the particulars of your experiment are. I'm often dealing with targeted sequencing, where we are trying to sequence a subset of a few thousand regions of a genome. For something like this, evaluating the best assembly is a little tricky. For each assembly kmer value, I blast my target regions against the assembly. I want to find the assembly that:

1) matches as many target regions as possible (i.e. lowest number of "No hits found" in the blast report)
2) maximizes the number of target regions that have contigs that match along their entire length (if a target region is 300 basepairs long, I want assemblies that make contigs where the alignment between the target and the contig is 300 base pairs long, or as close as possible).

This is more involved than counting with grep. It involves going through each blast query, denoting how long the query sequence is, getting the length of the alignment for the longest hit, and comparing the two. I use Bio::SearchIO from the bioPerl package for this sort of thing.

I'm not sure what the standard is, or if there is an agreed upon protocol for evaluating assemblies where the goals are more complex than maximizing N50. It's an active area of research for sure (i.e. http://www.biomedcentral.com/1471-2164/14/465). It may be appropriate to merge assemblies from different kmer values as well then reevaluate--something else to look into.

[mike.t]

This is a great opportunity for you to learn to program in Python or Perl!