Do you know a program/script to count how many sequences have mismatches on each position of my primers?

Question

I'm working with primer analysis and design, and so far I've performing it in Excel but: a) I don't like Excel, b) I do it manually, c) I don't like Excel, d) it gets tedious when working with thousands of sequences and tens of primers. My analysis consists of establishing how many sequences have mismatches on each position of my primers, and also count how many mismatches each sequence has with the corresponding primers (see attached image). This last analysis I've been able to (partially) solve it by using ncbi-blast+ tool. However, I will find very helpful a program/script where you input your sequences and primers fasta files and retrieve the number of mismatches at every primer position. Does anyone know something that could help me achieve all this? If it also has an output tab-separated file as the one attached that would be perfect for easy visualization of results. Suggestion/orientations are also grateful! I just want to know if there's something to help me out, otherwise I'll try writing a script to automatize my analysis. Thank you!

Answer 1

I'm amazed at your tenacity in getting Excel to do what you want; that looks like many hours of work to count mismatches.

The hammer I use at the moment for this sort of stuff is LAST, which will do local alignments of sequences to references. If the number of mismatches are not too large (as would be expected for primer sequences), a local alignment will extend to the end of one of the sequences, and I expect that should work fine for your purposes.

I'm going to demonstrate two options, one that can be done just using Excel and web services, and one that needs command-line access. Note that these are only two ways to fetch mismatch counts; you might want to do more than that (it's not clear from your question), and there are many other ways. Also, this may not work for your particular circumstances, but I think it should be quick enough that it's worth a try:

Option 1 - Web-based LAST

In this, I use the Web version of LAST to generate alignments for the primers [as "reference"] mapped to the sequences [as "query"].

1. To help create the mismatch (or in this case more correctly, "edit distance") statistic, in the Advanced Parameters section, set the a: gap existence cost to 1, and the b: gap extension cost to 1.
2. Click on Align!
3. Once the Alignment Results appear, click on the Alignments in tabular format link down the bottom, and save the file to your computer as last.tab.csv (the .csv at the end is the most important bit, to make sure that Excel recognises it as something it can open).
4. Open that file in Excel, set the delimiter to <tab>
5. The first column should be the alignment score, which in this case is the match score. This can be converted into a mismatch score by subtracting it from the primer sequence length (column 6).

Option 2 - Local LAST Installation

Installation

Install LAST and samtools.

Index Generation

The first step is to generate a reference index. LAST doesn't care too much about what's query and what's reference, but because I want all the genome sequences to appear on each line, the genome sequences will be the query, and the primer sequences will be the target. I use Martin Frith's lastdb mantra to do this:

$ lastdb -uNEAR -R01 primers.fa primers.fa

Mapping

The next step is to do the actual alignment with LAST. I set the match and mismatch to be both 1 (i.e. the default for the Web version), set the gap existence and extension to be 1, change the output format to tab, and save the result to last.tab.csv:

$ lastal -r 1 -q 1 -a 1 -b 1 primers.fa sequences.fa -f tab > last.tab.csv

Aggregating data

Finally, it's a matter of aggregating the count information. I use grep -v to exclude lines that begin with a #, use awk to subtract field 1 from field 6, then sort the results and use uniq -c to count the mismatches:

$ grep -v '^#' last.tab.csv | awk '{print ($6 - $1)}' | sort | uniq -c
2118 0
   5 1
   5 2