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 !

  • $\begingroup$ What is the purpose of simulating sequences from a genome, rather than simply using subsequences from the genome? If this is purely a programming exercise, the question would be better asked on a different StackExchange site. If there is a biological investigation behind it, you may be ignoring genomic context that is necessary for the proper function of the DNA sequence (e.g. see here). $\endgroup$
    – gringer
    Commented Jun 15, 2019 at 1:55
  • 2
    $\begingroup$ @gringer I want to generate sequences that will minimize homology to the target genome at the scale of the kmer, so I don't want to just pick subsequences. Right, it might have been a good question for stack overflow as well ! $\endgroup$
    – cmdoret
    Commented Jun 15, 2019 at 19:01

1 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.