I would like some advice on potential strategies to address the following problem.
I want to write a program that will generate DNA sequences that are optimized on two constraints based on an input genome:
- The GC content should be as close to the genome as possible.
- When taking all k-mers (of a predefined length) in the sequence, the sum of their frequencies in the genome should be as small as possible (i.e. it should be made of rare k-mers).
I currently score sequences using : $S = w_{GC} * (1 - |GC_{genome} - GC_{seq}|) + w_{Kmer} * (1 - \frac{mean(f(Kmer_{seq})}{max(f(Kmer_{genome}))}$
Where $w_{GC}$ and $w_{Kmer}$ are user-defined weights defining the relative importance of GC content and k-mer profile) and $f(Kmer_{seq})$ is the frequency of k-mers in the sequence.
Of course I could just use bruteforce and generate many random sequences of appropriate GC content, score them all and pick the very best. But 1) that's no fun and 2) it's slow.
The first thing I attempted was to extend the sequence one basepair at a time. Each iteration has 4 candidates (with a new A, C, T or G added). I pick the candidate using their score ($S$) as probability weight. This "blind" method is too simple and does not yield satisfactory results.
I thought a dynamic approach could work well (filling a table and backtracking), but I am not sure this is applicable since the cost of extending the sequence would be dependent on all new k-mers generated and not be a fixed value (such as in sequence alignment).
I also thought something like reinforcement learning could be appropriate here, but I don't know much about it, so I am not sure.
I know all of this is probably overkill, but this is mostly for the challenge / learning :)
Which strategy would be most appropriate for this ? Did I miss something simple ? All advice / suggestions are welcome !