<p>You can use a <span class="math-container">$l$</span>-order Markov chain. Here is the procedure:</p>

<ol>
<li><p>Count <span class="math-container">$l$</span>-mers in your genome. For small genomes, you can do that in Python. For large genomes, you may need jellyfish or KMC3.</p></li>
<li><p>Draw a <span class="math-container">$l$</span>-mer randomly based on the distribution of <span class="math-container">$l$</span>-mers. This generates the first <span class="math-container">$l$</span> bases.</p></li>
<li><p>Let <span class="math-container">$s$</span> be the last <span class="math-container">$(l-1)$</span> subsequence from the generated sequence. There will be at most four <span class="math-container">$l$</span>-mers prefixed with <span class="math-container">$s$</span>. It gives you a distribution of the next A/C/G/T base conditional on <span class="math-container">$s$</span>. Randomly draw a base from this distribution. You get the next base.</p></li>
<li><p>Repeat step 3 until you generate a sequence long enough.</p></li>
</ol>

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

<p>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.</p>

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

CC BY-SA 4.0
2019-06-16T16:37:25.837Z
2019-06-16T16:43:54.023Z
2019-06-16T16:43:54.023Z
user172818
<p>I would like some advice on potential strategies to address the following problem.</p>

<p>I want to write a program that will generate DNA sequences that are optimized on two constraints based on an input genome:</p>

<ul>
<li>The GC content should be as close to the genome as possible.</li>
<li>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).</li>
</ul>

<p>I currently score sequences using  :
<span class="math-container">$S = w_{GC} * (1 - |GC_{genome} - GC_{seq}|) + w_{Kmer} * (1 - \frac{mean(f(Kmer_{seq})}{max(f(Kmer_{genome}))}$</span></p>

<p>Where <span class="math-container">$w_{GC}$</span> and <span class="math-container">$w_{Kmer}$</span> are user-defined weights defining the relative importance of GC content and k-mer profile) and <span class="math-container">$f(Kmer_{seq})$</span> is the frequency of k-mers in the sequence.</p>

<p>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.</p>

<p>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 (<span class="math-container">$S$</span>) as probability weight. This "blind" method is too simple and does not yield satisfactory results.</p>

<p>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).</p>

<p>I also thought something like reinforcement learning could be appropriate here, but I don't know much about it, so I am not sure.</p>

<p>I know all of this is probably overkill, but this is mostly for the challenge / learning :)</p>

<p>Which strategy would be most appropriate for this ? Did I miss something simple ? All advice / suggestions are welcome !</p>

CC BY-SA 4.0
2019-06-15T01:55:43.117Z
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](https://twitter.com/gringene_bio/status/993821037513994242)).
gringer
CC BY-SA 4.0
2019-06-15T19:01:15.987Z
@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 !
cmdoret
<sequence-analysis><genomics><k-mer>
Generating DNA sequences with constraints
# Generating DNA sequences with constraints

(Asked by: username_0 on 2019-06-15)

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 !

**Tags:** sequence-analysis, genomics, k-mer
> **Comments:**
> **username_1**: 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](https://twitter.com/gringene_bio/status/993821037513994242)).  
> **username_0**: @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 !  

---

## 1 answers

### Answer by username_2 on 2019-06-16. Score: 4 | Accepted answer

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.


---

