Generating DNA sequences with constraints

Question

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 !

Answer 1

You can use a $l$-order Markov chain. Here is the procedure:

1. Count $l$-mers in your genome. For small genomes, you can do that in Python. For large genomes, you may need jellyfish or KMC3.

2. Draw a $l$-mer randomly based on the distribution of $l$-mers. This generates the first $l$ bases.

3. Let $s$ be the last $(l-1)$ subsequence from the generated sequence. There will be at most four $l$-mers prefixed with $s$. It gives you a distribution of the next A/C/G/T base conditional on $s$. Randomly draw a base from this distribution. You get the next base.

4. Repeat step 3 until you generate a sequence long enough.

The generated sequence will have the same $l$-mer distribution to the original sequence. The GC content will be the same, too. In your case, if $k$ is small enough, you can let $l=k$. If $k$ is large, you have to choose a smaller $l$.

The procedure above doesn't consider the strand symmetry of DNA, but probably this doesn't matter much. It is also possible to adapt the algorithm to address this.

As to implementation, you can use python, but a high performance language like C will be tens of times faster at least.