How to optimize the number of amplicons ordered for a PCR wet experiment with several genomic ROIs?

Question

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:

Assume we have a precomputed set of all necessary amplicons 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)

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?

Answer 1

I ended up implementing the following greedy algorithm that should work fairly good for any 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)