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

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:
for c in cover:
if does_a_cover_b(c,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?

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):

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):
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)