1

I have a sequence of chars that I multialign with the Matlab function 'multialign'. The result is a char matrix with the multialigned sequences:

e.g. with only 3 seqences.

----GC
AT--GC
ATGCGC

Next I score every column of the alignment on the base of the similarity between the sequences. In the previous example: in the fist column the 'A' is the most frequent symbol and it appears 2 out of 3 times, so the score is 2/3, in the last column the letter 'C' appears every time so I score it 1, and so on. The final score is the mean of the single scores.

Now the real question: I create these sequences with a function that uses a threshold to decide whether to add a char or not, so I can have shorter or longer sequences. This is another example with a different threshold that I score in the same way:

-----ATATGGCGC
AT-ATGCA-G-C--
ATG-TGC--G-C--

I wanted to use fminsearch to search for the best threshold but my problem is that varying the threshold only a little the score doesn't change so the algorithm used by fminsearch doesn't work (for example if you start from 10 the value chosen for the next step is something like 9.75...).

This is a pseudo-Matlab code for what I do:

[bestthr, bestscore] = fminsearch(@(x) fcnthr(data,x),[10]);

function score = fcnthr(data,thr)

sequences = generateSequencesFromData(data,thr);
multialignmatrix = multialign(sequences);
score = scorealignments(multialignmatrix);

score = 1/score; %I want to find the maximum score so fminsearch searches the minimum of 1/score.

end

Can someone help me?

Marco
  • 1,454
  • 1
  • 16
  • 30
  • genetic algorithms (ga) may help with its integer programming support, or even brute-force for-loop may be useful – Dmitry Galchinsky Oct 05 '13 at 18:01
  • Thank you for your answer. Brute force has been my first solution but it's time consuming and I need to choose an upper bound to stop searching so in these days I have been searching for a more performing solution. I also thought of using GAs but they tend to be slow, however maybe I can give them a try :) – Marco Oct 05 '13 at 18:41

0 Answers0