This is a test version of Biostars. For the public version, visit https://www.biostars.org.
BLAST run parameters and parsing advice

I have a set of transposon element(TE) sequences and their genomic coordinates.

For each TE, I have concatenated a mosaic of their flanking 5' and 3' sequences, up to 50nt on each end, so up to 100nt total mosaic length. Let's call this the "pre-insertion" sequence.

Now, for each mosaic as query, I need to report coordinates and copy numbers for regions in the genome with matches at >= x% identity, over at least y% of query length (excluding the one self-match).

My questions are:

A. how best to run the BLAST search, and (outfmt 6 or 11 or what?)

B. Importantly, how best to parse the BLAST search results to report copy number and genomic coordinates for each mosaic query.

C. To reiterate, I do not want to report the gapped self-match (where the gap is the TE sequence). But I realize there could be other gapped matches, where the intervening sequence is not the original TE sequence, but some other sequence. In principle, by retaining coords of query in some manner and comparing them to BLAST results, I could filter out self-match cases, right?

Is mine too specific a parsing requirement, or will any off-the-shelf tools like blast_formatter or bp_search2gff.pl work out for me? I've never had to parse BLASTn output before, so this is new to me. Looking forward to your advice and suggestions. Thanks all!

In the example visualization, you can see on chr 1 and 4, there are instances of pre-insertion sequences - I want to report each of them with their start-stop coords.

On chr 1, there is another match with a black box in between, these matches will be reported by blastN. But those matches will be separated by a non-matching, non-TE sequence. Therefore, I will not consider this a pre-insertion sequence, and want to ignore such cases.

On Chr 2, on the minus strand, and on Chr 3, on the plus strand, are matches to the flanking sequences, separated by original TE sequence. If you imagine the entire region duplicated from Chr 2 to Chr 3, this can be possible. Here again, these do not represent pre-insertion sequence, and I want to ignore such cases.

Insertional-Mosaic-BLASTn-Vs-Genome

blast blastn blastparser
• 1,661 views • link
written by Anand Rao Anand Rao

Hello Anand Rao!

upon your request we closed this question (not suited for your needs)

• 1 views • link
written by lieven.sterck lieven.sterck

1 answer

I'm not quite understanding your set-up here (I'll come back to that) but to answer some of your questions already:

A . output format 6 or 7 (the tabular ones) are the preferred choice here since those allow easy parsing afterwards

B. This will require some very basic scripting (perl? python?) but also awk for instance will get you a long way. Are you familiar with any of those.

C. That is a correct approach indeed, but might be a bit cumbersome though. Depending on the size of the TEs you started with I would be surprised to see any of those gapped self-matches reported (given that you will only have a short alignable part and a very large gap, blastn will likely never report this as a valid match). You should in any way also consider this for the 'inverse' cases (== those you would like to se reported), so count on potential substantial false negative rate in this analysis.

As I mentioned in the beginning: I'm not sure what you want to achieve here. Could you better explain what you mean with those mosaic sequences? Do you take 50n upstream of the start of the TE and 50n downstream of the end and concatenate those?

• 1 views • link
written by lieven.sterck lieven.sterck

Thanks Lieven. I've added an image and some lines in italics to help me visualize match cases better. Look fwd to your response.

• 1 views • link
written by Anand Rao Anand Rao

OK, very clear picture to illustrate your question/setup.

This might in theory work indeed, however, I seem to believe that the insertion site determination happens on a much smaller region than the 50n (100n) you use. Moreover, there are many other factors playing for TE insertion sites than only the actual sequence (eg. methylation state)

• 1 views • link
written by lieven.sterck lieven.sterck

I realized local alignment tools (such as BLASTn that use the Smith Waterman algo) are not the way to go, I need a global end-to-end alignment tool (like VMATCH or Bowtie2 that use the Needleman Wunsch algo) So I have to retract my question, since it's premise of using BLASTn is itself wrong! ha :}

• 1 views • link
written by Anand Rao Anand Rao

OK, fair enough.

I will close your question then (not delete it), you never know it might be useful to someone.

and one last comment on your blast cmdline: -perc_identity and -qcov_hsp_perc are real numbers , so in your example you filter on 0.8% identity (and not 80 !)

• 1 views • link
written by lieven.sterck lieven.sterck

Log in to answer this question.

More posts like this