doi:10.1093/bioinformatics/bth147
“Primer design using genetic algorithm”
一、背景&问题什么是引物设计，引物设计很重要，是PCR实验（聚合酶链式反应，是一种让DNA序列以指数方式快速扩增的实验）的前置动作，引物的质量好坏影响PCR实验的成功与否。迄今为止，已经开发了许多基于经典算法的引物设计方法。
在进行聚合酶链式反应实验之前，需要一对引物来剪切目标DNA子序列。然而，这是一项繁琐的任务，因为需要满足太多的约束。已经提出了各种设计引物的方法，但是它们中的大多数在设计的引物上没有限制性位点，并且不满足特异性限制。
二、解决方案在给定的模板序列里，生成种群进行引物对的随机搜索，并计算对应的适应度，选择拥有优秀适应度的作为子代，模拟染色体的交叉、变异过程。待迭代次数足够且种群结果收敛后，得到最优的引物对结果，该引物对具有特异性和限制性位点。
三、总流程模板：TP53基因的CDS部分区域，长度为587
产物长度：80~150
交叉概率：0.6
变异概率：0.1
（一）初始化种群初始化过程随机产生群体P的初始个体，群体P有N个解Pt，一个解对应一对引物，解Pt被定义为向量：Pt=(Fs，Fe，Rs，Re)。但由于四个组成部分是依赖的（其中一个在交叉或变异过程中被改变，有时可能导致引物或产物违反长度约束）。
为了避免这个问题，将依赖性质转化为独立的形式，Pt=(Fs，α，β，γ)。其中α=Fe-Fs， β=Rs-Fe， γ=Re-Rs。
Fs:上游引物的起始位点，0≤Fs≤模板长度-最小产物长度
α:上游引物长度，18≤α≤26
β:上游引物结束位点与下游引物起始位点之间的距离，最小产物长度-α-γ≤β≤最大产物长度-α-γ and Fs+α+β+γ≤模板长度
γ:下游引物长度，18≤γ≤26
（二）适应度函数&选择适应度函数：
子函数满足要求置0，不满足要求置1
leng:长度（18~26）
lengd:长度差（≤3）
Tmd:Tm差（≤5）
퐺퐶푝:GC含量（0.4~0.6）
Term:引物3’末端是否为G或C
Uni:引物特异性（每条引物与模板序列之间只有一处退火位置）
Sc:引物不能有自身互补序列，否则产生发夹结构
PC:一个引物不应与另一个引物互补配对，否则产生引物二聚体
푅푡:是否有限制性位点,进行模式匹配，验证酶是否在其上。其中模式匹配长度：
文章所使用的限制性酶及其核酸序列：
选择：
对计算完适应度值的引物组进行筛选，其fitness越低，引物越符合要求。取fitness的倒数，将反比例关系转成正比例关系后，进行筛选（轮盘赌）。
每个引物对被选择作为下一代的机率：
（三）交叉&变异交叉：
产生小于16的十进制随机数，转成四位二进制数后，将位数为1的位置进行互换，得到新的引物对。对新得到的引物对Pt=(Fs，α，β，γ)进行检查，若同时满足产物长度合理和下游引物结束位点不越界，则进入下一步；否则重新开始交叉过程。
eg:
变异：
随机改变四个部分中的一个，按照初始化时的变量范围来赋值，然后检查变异后的个体是否符合约束。检查的方式同上。
eg:
流程图：
四、运行结果有提前出现适应度最优的引物即fitness=0的情况，且解有时不唯一，如下图所示（很可能是设定的产物长度明显小于模板长度的缘故）
五、代码核心部分fitness函数：
def each_fitness(pt):#适应度函数
return leng(pt)+3*lengd(pt)+3*Tmd(pt)+3*GCp(pt)+3*Term(pt)+50*Uni(pt)+10*Sc(pt)+10*PC(pt)+Rt(pt)
def leng(pt):#长度限制
if 18<=len(Bf(pt)) <= 26 and 18 <= len(Br(pt)) <= 26:
return 0
else:
return 1
def lengd(pt):#长度差限制
if abs(len(Bf(pt))-len(Br(pt)))<=3:
return 0
else:
return 1
def Tmd(pt):#Tm差限制
if abs(Tm(Bf(pt))-Tm(Br(pt)))<=5:
return 0
else:
return 1
def GCp(pt):#GC含量限制
if 0.4 <= GC(Bf(pt)) <= 0.6 and 0.4 <= GC(Br(pt)) <= 0.6:
return 0
else:
return 1
def Term(pt):#两引物3'末端情况
return Termf(pt)+Termr(pt)
def Termf(pt):#上游引物3'末端情况
f3=pt[0]+pt[1]-1
if GDC[f3] == 'g' or GDC[f3] == 'c':
return 0
else:
return 1
def Termr(pt):#下游引物3'末端情况
r3=pt[0]+pt[1]+pt[2]
if GD[r3] == 'g' or GD[r3] == 'c':
return 0
else:
return 1
def Uni(pt):#两引物特异性情况
return Unif(pt)+Unir(pt)
def Unif(pt):#上游引物特异性情况（模式匹配和GDC）
if GDC.count(Bf(pt)) == 1:
return 0
else:
return 1
def Unir(pt):#下游引物特异性情况（模式匹配和GD）
if GD.count(Br(pt)) == 1:
return 0
else:
return 1
def Sc(pt):#自互补
tempf=0
tempr=0
#上游引物情况
for i in range(0, len(Bf(pt))):
strf = Bf(pt)[i:i + 4]
if Bf(pt).find(reverse(complete(strf)),i+4)!=-1:#这里的4可以调，以调整发夹不配对圈的大小
tempf=1
break
#下游引物情况
for i in range(0, len(Bf(pt))):
strr = Br(pt)[i:i + 4]
if Bf(pt).find(reverse(complete(strr)),i+4)!=-1:#这里的4可以调，以调整发夹不配对圈的大小
tempr=1
break
if tempf==0 and tempr==0:
return 0
else:
return 1
def PC(pt):#间互补
if getNumofCommonSubstr(Bf(pt),complete(Br(pt))) <= 3 and getNumofCommonSubstr(Bf(pt),reverse(complete(Br(pt)))) <= 3:
return 0
else:
return 1
def Rt(pt):#搜索限制位点
temp=1
# 上游引物限制位点
for i in similar_list:
if Bf(pt).count(i)!=0:
temp=0
break
# 下游引物限制位点
for i in similar_list:
if reverse(Bf(pt)).count(i)!=0:
temp=0
break
if temp==0:
return 0
else:
return 1
def get_fitness(mat):#获得种群适应度值清单
fitness=np.zeros((pop_size,1))
for i in range(pop_size):
#print(mat[i])
fitness[i]=each_fitness(mat[i])
return fitness
初始化种群：
def init_mat():#初始化种群
mat=np.zeros((pop_size,vector_size))
for i in range(0,pop_size):
#初始化生成vector
Fs=np.random.randint(0,template_length-min_length+1)
alpha=np.random.randint(18,27)
gamma=np.random.randint(18,27)
if max_length-alpha-gamma+1
选择：
def select(pop, profitness):#选择最优（根据fitness的倒数）
idx = np.random.choice(np.arange(pop_size), size=pop_size, replace=True,
p=(profitness) / (profitness.sum()))
better_pop = pop[idx]
return better_pop
交叉&变异：
def crossover_and_mutation(pop):
new_pop = []
for father in pop: # 遍历种群中的每一个个体，将该个体作为父亲
# print(type(father)) ndarray
# print(father.shape) (48,)
child = father # 孩子先得到父亲的全部基因（这里我把一串二进制串的那些0，1称为基因）
# 产生子代时不是必然发生交叉，而是以一定的概率发生交叉 0~Pe
if np.random.rand() < Pe:
child=crossover(child,pop)
# 每个后代有一定的机率发生变异,以一定的概率发生变异 0~Pm
if np.random.rand() < Pm:
child = mutation(child)
new_pop.append(child)
return new_pop #list类型
def crossover(child,pop):
while True:
mother = pop[np.random.randint(pop_size)] # 再种群中选择另一个个体，并将该个体作为母亲,类型与father一样
seq=bina(np.random.randint(0,16))
for i in range(4):
if seq[i]==1:
temp=child[i]
child[i]=mother[i]
mother[i]=temp
if 80<=child[1]+child[2]+child[3]<=150 and child[0]+child[1]+child[2]+child[3]<=template_length:
break
return child
def mutation(child):
while True:
index=np.random.randint(0,4)
if index==0:
child[0]=np.random.randint(0,template_length-min_length+1)
if index==1:
child[1]=np.random.randint(18,27)
if index==2:
if max_length - child[1] - child[3] + 1 < template_length - child[0] - child[1] - child[3] + 1:
beta = np.random.randint(min_length - child[1] - child[3], max_length - child[1] - child[3] + 1)
else:
beta = np.random.randint(min_length - child[1] - child[3], template_length - child[0] - child[1] - child[3] + 1)
if index==3:
child[3]=np.random.randint(18,27)
if 80<=child[1]+child[2]+child[3]<=150 and child[0]+child[1]+child[2]+child[3]<=template_length:
break
return child 
