Answers.Body: <p>I ended up implementing the following greedy algorithm that should work fairly good for any <a href="https://en.wikipedia.org/wiki/Binary_relation" rel="nofollow noreferrer">binary relation</a> even if it is not transitive nor symmetric.</p>
  <p>Assume you have the following binary relation (you may plug in any function you want):</p>
  <pre><code>def does_a_cover_b(a,b, min_dist=30, max_dist=130):
    return min_dist &lt; b.left_roi - a.left_end &lt; max_dist or min_dist &lt; a.right_end - b.right_roi &lt; max_dist
  </code></pre>
  <p>Using the auxiliary functions to define the binary relation, draw the graph, and build the graph:</p>
  <pre><code>import networkx as nx
  
  def does_a_cover_b(a,b, min_dist=30, max_dist=130):
    return min_dist &lt; b.left_roi - a.left_end &lt; max_dist or min_dist &lt; a.right_end - b.right_roi &lt; max_dist
  
  def draw(G):
    pos = nx.spring_layout(G,scale=2)
    nx.draw(G, pos, with_labels=True, font_weight='bold', node_size=750, width=2)
  
  def build_a_graph(primer_list):
    G = nx.DiGraph()
  
    for i,p_i in enumerate(primer_list):
      G.add_node(i, data=p_i)
  
    for i,p_i in enumerate(primer_list):
      for j,p_j in enumerate(primer_list):
        if i == j:
          continue
        if does_a_cover_b(p_i, p_j):
          G.add_edge(i,j)
    return G
  
  </code></pre>
  <p>we may implement the following algorithm:</p>
  <pre><code>p_set = set(primers)
  cover_set = set()
  covered_primers = set()
  while len(covered_primers)&lt;N and len(p_set)&gt;0:
    G = build_a_graph(p_set)
    nodes =sorted(G.out_degree, key=lambda item: item[1], reverse=True)
    node = nodes[0][0]
    absorber = G.nodes[node]['data']
    p_set.remove(absorber)
    cover_set.add(absorber)
    covered_primers.add(absorber)
  
    for edge in G.out_edges(node):
      primer = G.nodes[edge[1]]['data']
      p_set.remove(primer)
      covered_primers.add(primer)
  
  print(f&quot;Used a cover set of length {len(cover_set)} primers to cover the total primer set with {len(covered_primers)} primers&quot;)
  plot_primer_set(primers)
  plot_primer_set(cover_set)
  </code></pre>
  <p><a href="https://i.sstatic.net/ePf7n.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ePf7n.png" alt="enter image description here" /></a></p>
  <p><a href="https://i.sstatic.net/3qHvc.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/3qHvc.png" alt="enter image description here" /></a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-25T18:55:28.173Z
Answers.LastActivityDate: 2021-09-26T09:36:25.863Z
Answers.LastEditDate: 2021-09-26T09:36:25.863Z
Answers.UserDisplayName: 0x90
Body: <p>Imagine we have an experiment for which we would like to minimize the number of purchased primers. Let's look at the layout of a single primer below:</p>
  <p><a href="https://i.sstatic.net/M4Tyz.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/M4Tyz.png" alt="enter image description here" /></a></p>
  <p>Assume we have a precomputed set of all necessary <a href="https://en.wikipedia.org/wiki/Amplicon" rel="nofollow noreferrer">amplicons</a> and we are trying to reduce the costs and find the minimum set of amplicons that tile all regions of interest (ROIs) with some flanking.</p>
  <p>The algorithm seems to work, but I am not sure it is ideal, so I would like to know how to assert it finds the minimal covering set of amplicons and how to improve the runtime complexity.</p>
  <p>First, we define the underline representation of a primer/amplicon (I  use the terms interchangeably), and an auxiliary function for maintaining the same colors among different plotting sessions.</p>
  <pre><code>from dataclasses import dataclass
  import numpy as np
  
  def rand_color():
     return &quot;#&quot; + &quot;%06x&quot; % np.random.randint(0, 0xFFFFFF)
  
  @dataclass(eq=True,frozen=True)
  class Primer:
      left_end:  int      
      right_end: int
      left_roi: int
      right_roi: int
      color: str
      level: int
  </code></pre>
  <p>We then create a random set of primers:</p>
  <pre><code>N = 25; max_coor = 2_000
  min_primer_len, max_primer_len = 200, 500
  min_required_flank_len, max_required_flank_len = 10, 50
  roi_len = 20
  primers = []
  for i in range(N):
    left_end = np.random.randint(0,max_coor-max_primer_len)
    primer_length = np.random.randint(min_primer_len,max_primer_len)
    required_roi_flank = np.random.randint(min_required_flank_len,max_required_flank_len)
    left_roi = np.random.randint(left_end+required_roi_flank,left_end+primer_length-required_roi_flank-roi_len-1)
    right_roi = left_roi + roi_len
    primers.append(Primer(left_end, left_end+primer_length, left_roi, right_roi, color=rand_color(), level=i))
  
  </code></pre>
  <p>The algorithm for finding the minimal primer covering set is as follows:</p>
  <pre><code>primers.sort(key=lambda p:p.left_end)
  
  def does_a_cover_b(a,b, flank=0):
    return a.left_end&lt; b.left_roi-flank and a.right_end&gt; b.right_roi+flank  
  
  p_tmp = primers.copy()
  
  cover = set()
  for p in p_tmp:
    add_to_cover_set = True
    for c in cover:
      if does_a_cover_b(c,p):
        add_to_cover_set = False
    
    if add_to_cover_set:
      cover.add(p)
  </code></pre>
  <p>Now <code>cover</code> is supposed to be a minimal subset of primers that covers all ROIs. To visualize, we use <code>plot_primer_set</code> twice:</p>
  <pre><code>import matplotlib.pyplot as plt
  from shapely.geometry import Polygon
  
  def plot_primer_set(primer_set):
    plt.figure(figsize=(4,4), dpi=120)
    for c in primer_set:
      p_all = Polygon([(c.left_end,c.level),(c.left_roi,c.level),(c.right_roi,c.level),(c.right_end,c.level)])
      x,y = p_all.exterior.xy
      plt.plot(x,y, markersize=7, marker='o', linewidth=5, color=c.color)
      plt.yticks(np.arange(0, N, 1))
      axes = plt.axes()
      axes.set_xlim([-100, max_coor+100])
  
  plot_primer_set(primers)
  plot_primer_set(cover)
  
  </code></pre>
  <p><a href="https://i.sstatic.net/Kjzjy.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/Kjzjy.png" alt="enter image description here" /></a>
  <a href="https://i.sstatic.net/1cnJn.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/1cnJn.png" alt="enter image description here" /></a></p>
  <ol>
  <li>Is this the optimal and correct minimal set?</li>
  <li>My complexity is <span class="math-container">$O(N^2)$</span>, but I expected it to be <span class="math-container">$O(N \log{N})$</span>, what am I missing?</li>
  </ol>
  
Tags: <python><algorithms><pcr>
Title: How to optimize the number of amplicons ordered for a PCR wet experiment with several genomic ROIs?
ThreadText: # How to optimize the number of amplicons ordered for a PCR wet experiment with several genomic ROIs?
  
  (Asked by: username_0 on 2021-09-22)
  
  Imagine we have an experiment for which we would like to minimize the number of purchased primers. Let's look at the layout of a single primer below:
  
  [![enter image description here](https://i.sstatic.net/M4Tyz.png)](https://i.sstatic.net/M4Tyz.png)
  
  Assume we have a precomputed set of all necessary [amplicons](https://en.wikipedia.org/wiki/Amplicon) and we are trying to reduce the costs and find the minimum set of amplicons that tile all regions of interest (ROIs) with some flanking.
  
  The algorithm seems to work, but I am not sure it is ideal, so I would like to know how to assert it finds the minimal covering set of amplicons and how to improve the runtime complexity.
  
  First, we define the underline representation of a primer/amplicon (I use the terms interchangeably), and an auxiliary function for maintaining the same colors among different plotting sessions.
  
  ```
  from dataclasses import dataclass
  import numpy as np
  
  def rand_color():
     return "#" + "%06x" % np.random.randint(0, 0xFFFFFF)
  
  @dataclass(eq=True,frozen=True)
  class Primer:
      left_end:  int      
      right_end: int
      left_roi: int
      right_roi: int
      color: str
      level: int
  
  ```
  
  We then create a random set of primers:
  
  ```
  N = 25; max_coor = 2_000
  min_primer_len, max_primer_len = 200, 500
  min_required_flank_len, max_required_flank_len = 10, 50
  roi_len = 20
  primers = []
  for i in range(N):
    left_end = np.random.randint(0,max_coor-max_primer_len)
    primer_length = np.random.randint(min_primer_len,max_primer_len)
    required_roi_flank = np.random.randint(min_required_flank_len,max_required_flank_len)
    left_roi = np.random.randint(left_end+required_roi_flank,left_end+primer_length-required_roi_flank-roi_len-1)
    right_roi = left_roi + roi_len
    primers.append(Primer(left_end, left_end+primer_length, left_roi, right_roi, color=rand_color(), level=i))
  
  ```
  
  The algorithm for finding the minimal primer covering set is as follows:
  
  ```
  primers.sort(key=lambda p:p.left_end)
  
  def does_a_cover_b(a,b, flank=0):
    return a.left_end< b.left_roi-flank and a.right_end> b.right_roi+flank  
  
  p_tmp = primers.copy()
  
  cover = set()
  for p in p_tmp:
    add_to_cover_set = True
    for c in cover:
      if does_a_cover_b(c,p):
        add_to_cover_set = False
    
    if add_to_cover_set:
      cover.add(p)
  
  ```
  
  Now `cover` is supposed to be a minimal subset of primers that covers all ROIs. To visualize, we use `plot_primer_set` twice:
  
  ```
  import matplotlib.pyplot as plt
  from shapely.geometry import Polygon
  
  def plot_primer_set(primer_set):
    plt.figure(figsize=(4,4), dpi=120)
    for c in primer_set:
      p_all = Polygon([(c.left_end,c.level),(c.left_roi,c.level),(c.right_roi,c.level),(c.right_end,c.level)])
      x,y = p_all.exterior.xy
      plt.plot(x,y, markersize=7, marker='o', linewidth=5, color=c.color)
      plt.yticks(np.arange(0, N, 1))
      axes = plt.axes()
      axes.set_xlim([-100, max_coor+100])
  
  plot_primer_set(primers)
  plot_primer_set(cover)
  
  ```
  
  [![enter image description here](https://i.sstatic.net/Kjzjy.png)](https://i.sstatic.net/Kjzjy.png)
  [![enter image description here](https://i.sstatic.net/1cnJn.png)](https://i.sstatic.net/1cnJn.png)
  
  1. Is this the optimal and correct minimal set?
  2. My complexity is $O(N^2)$, but I expected it to be $O(N \log{N})$, what am I missing?
  
  **Tags:** python, algorithms, pcr
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_0 on 2021-09-25. Score: 1 | Accepted answer
  
  I ended up implementing the following greedy algorithm that should work fairly good for any [binary relation](https://en.wikipedia.org/wiki/Binary_relation) even if it is not transitive nor symmetric.
  
  Assume you have the following binary relation (you may plug in any function you want):
  
  ```
  def does_a_cover_b(a,b, min_dist=30, max_dist=130):
    return min_dist < b.left_roi - a.left_end < max_dist or min_dist < a.right_end - b.right_roi < max_dist
  
  ```
  
  Using the auxiliary functions to define the binary relation, draw the graph, and build the graph:
  
  ```
  import networkx as nx
  
  def does_a_cover_b(a,b, min_dist=30, max_dist=130):
    return min_dist < b.left_roi - a.left_end < max_dist or min_dist < a.right_end - b.right_roi < max_dist
  
  def draw(G):
    pos = nx.spring_layout(G,scale=2)
    nx.draw(G, pos, with_labels=True, font_weight='bold', node_size=750, width=2)
  
  def build_a_graph(primer_list):
    G = nx.DiGraph()
  
    for i,p_i in enumerate(primer_list):
      G.add_node(i, data=p_i)
  
    for i,p_i in enumerate(primer_list):
      for j,p_j in enumerate(primer_list):
        if i == j:
          continue
        if does_a_cover_b(p_i, p_j):
          G.add_edge(i,j)
    return G
  
  ```
  
  we may implement the following algorithm:
  
  ```
  p_set = set(primers)
  cover_set = set()
  covered_primers = set()
  while len(covered_primers)<N and len(p_set)>0:
    G = build_a_graph(p_set)
    nodes =sorted(G.out_degree, key=lambda item: item[1], reverse=True)
    node = nodes[0][0]
    absorber = G.nodes[node]['data']
    p_set.remove(absorber)
    cover_set.add(absorber)
    covered_primers.add(absorber)
  
    for edge in G.out_edges(node):
      primer = G.nodes[edge[1]]['data']
      p_set.remove(primer)
      covered_primers.add(primer)
  
  print(f"Used a cover set of length {len(cover_set)} primers to cover the total primer set with {len(covered_primers)} primers")
  plot_primer_set(primers)
  plot_primer_set(cover_set)
  
  ```
  
  [![enter image description here](https://i.sstatic.net/ePf7n.png)](https://i.sstatic.net/ePf7n.png)
  
  [![enter image description here](https://i.sstatic.net/3qHvc.png)](https://i.sstatic.net/3qHvc.png)
  
  
  ---
  
  