基于动态规划算法的DNA序列比对函数，给出两条序列（v和w）的打分矩阵
一.什么是动态规划算法
1.1总体思想
·动态规划算法与分治法类似，基本思想也是将待求解的问题分成若干个子问题
·经过分解得到的子问题往往不是互相独立的，有些子问题被重复计算多次
·如果能够保存已解决的子问题答案，在需要时再找出来已求得的答案，就可以避免大量重复计算，从而得到多项式时间算法（备忘录）
1.2使用动态规划求解的问题需要具备的基本要素
1）重复子问题
·递归算法求解问题时，每次产生的子问题并不总是新问题，有些子问题被反复计算多次，这种性质被称为子问题的重叠性质
·动态规划算法，对每一个子问题只解一次，而后将其解保存在一个表格中，当再次需要解此子问题时，只是简单地用常数时间查看一下结果。
·通常不同的子问题个数随问题的大小呈多项式增长，用动态规划算法只需要多项式时间，从而获得较高的解题效率。
2）最优子结构
·一个问题的最优解包含着其子问题的最优解，这种性质称为最优子结构性质
·分析问题的最优子结构性质，首先假设由问题的最优解导出的子问题的解不是最优，然后再设法说明在这个假设下可构造出比原问题最优解更好的解，从而导致矛盾
·利用问题的最优子结构性质，以自底向上的方式递归地从子问题的最优解逐步构造出整个问题的最优解
·最优子结构是一个问题能用动态规划算法求解的前提
1.3动态规划求解的基本步骤
1）找出最优解的性质，并刻画其结构特征
2）递归地定义最优质
3）以自底向上的方式计算出最优值
4）根据计算最优值时得到的信息，构造最优解
二.打分矩阵代码
import random # 1. 生成一个指定长度的随机DNA序列
def generate_dna(length): dna_bases = 'ACGT' # DNA的四个碱基 sequence = '' # 初始化空序列 for _ in range(length): # 循环length次 sequence += random.choice(dna_bases) # 每次从四个碱基中随机选一个添加到序列中 return sequence # 返回生成的DNA序列 # 2. 在DNA序列中随机位置插入一个可能突变的motif
def insert_motif(dna, motif, mutation_rate): mutated_motif = '' # 初始化突变的motif为空字符串 for base in motif: # 遍历motif中的每个碱基 if random.random() < mutation_rate: # 如果随机数小于突变率 mutated_motif += random.choice('ACGT') # 则该位置碱基随机突变 else: mutated_motif += base # 否则保持原样 insert_point = random.randint(0, len(dna)) # 在DNA序列中随机选择一个插入点 return dna[:insert_point] + mutated_motif + dna[insert_point:] # 插入突变的motif # 3. 生成多条带有随机插入motif的DNA序列
def generate_sequences(num_sequences, dna_length, motif_length, mutation_rate): sequences = [] # 初始化一个空列表来存储生成的序列 motif = generate_dna(motif_length) # 生成一个随机的motif for _ in range(num_sequences): # 循环生成指定数量的序列 dna = generate_dna(dna_length) # 生成一个DNA序列 dna_with_motif = insert_motif(dna, motif, mutation_rate) # 插入motif sequences.append(dna_with_motif) # 将序列添加到列表中 return sequences # 返回生成的序列列表 # 使用函数生成序列并打印
sequences = generate_sequences(5, 20, 5, 0.1)
for seq in sequences: print(seq)三.编写一个函数，生成m条DNA序列，每条序列长度为k，然后对每条序列随机插入一个长度为L的motif，motif的突变率为n
import random # 1. 生成一个指定长度的随机DNA序列
def generate_dna(length): dna_bases = 'ACGT' # DNA的四个碱基 sequence = '' # 初始化空序列 for _ in range(length): # 循环length次 sequence += random.choice(dna_bases) # 每次从四个碱基中随机选一个添加到序列中 return sequence # 返回生成的DNA序列 # 2. 在DNA序列中随机位置插入一个可能突变的motif
def insert_motif(dna, motif, mutation_rate): mutated_motif = '' # 初始化突变的motif为空字符串 for base in motif: # 遍历motif中的每个碱基 if random.random() < mutation_rate: # 如果随机数小于突变率 mutated_motif += random.choice('ACGT') # 则该位置碱基随机突变 else: mutated_motif += base # 否则保持原样 insert_point = random.randint(0, len(dna)) # 在DNA序列中随机选择一个插入点 return dna[:insert_point] + mutated_motif + dna[insert_point:] # 插入突变的motif # 3. 生成多条带有随机插入motif的DNA序列
def generate_sequences(num_sequences, dna_length, motif_length, mutation_rate): sequences = [] # 初始化一个空列表来存储生成的序列 motif = generate_dna(motif_length) # 生成一个随机的motif for _ in range(num_sequences): # 循环生成指定数量的序列 dna = generate_dna(dna_length) # 生成一个DNA序列 dna_with_motif = insert_motif(dna, motif, mutation_rate) # 插入motif sequences.append(dna_with_motif) # 将序列添加到列表中 return sequences # 返回生成的序列列表 # 使用函数生成序列并打印
sequences = generate_sequences(5, 20, 5, 0.1)
for seq in sequences: print(seq)四.对生成的已插入突变motif的序列集合，编写一套函数，寻找其中的motif,可指定motif长度;
# 定义一个名为 find 的函数，它接受两个参数：
# sequence_set 是一个包含多个 DNA 序列的列表
# motif_length 是我们想要查找的子序列（也称为 motif）的长度
def find(sequence_set, motif_length): # 初始化一个空集合 motifs，用于存储找到的所有唯一的 motif motifs = set() # 遍历 sequence_set 中的每一条序列 for sequence in sequence_set: # 对于每一条序列，我们从其第一个碱基开始，直到剩下的碱基数量少于 motif_length 为止 # 这样确保我们可以从序列中提取出完整长度的 motif for i in range(len(sequence) - motif_length + 1): # 从当前位置 i 开始，提取长度为 motif_length 的子序列 motif = sequence[i:i+motif_length] # 将提取到的 motif 添加到 motifs 集合中 # 由于 motifs 是一个集合，所以重复的 motif 会被自动去除 motifs.add(motif) # 函数返回包含所有唯一 motif 的集合 return motifs # 测试数据：一个包含三条 DNA 序列的列表和一个 motif 长度值
sequence_set = ['ACGTTAGC', 'GTATCGAG', 'CGTACGTA']
motif_length = 4 # 调用 find 函数，传入测试数据，并将返回的结果存储在变量 motifs 中
motifs = find(sequence_set, motif_length)
# 打印出找到的所有唯一 motif
print(motifs)五.对生成的已插入突变motif的序列集合，编写一套函数，基于分支界定法寻找指定长度的motif,并与遍历法比较计算效率
import itertools
import time # 计算motif在序列中的得分
def calculate_score(motif, sequences): score = 0 for seq in sequences: min_distance = float('inf') for i in range(len(seq) - len(motif) + 1): distance = sum(motif[j] != seq[i + j] for j in range(len(motif))) min_distance = min(min_distance, distance) score += min_distance return score # 分支界定法
def find_motif_branch_bound(sequences, motif_length): best_motif = None best_score = float('inf') def search(motif, depth): nonlocal best_motif, best_score if depth == motif_length: score = calculate_score(motif, sequences) if score < best_score: best_score = score best_motif = motif return for base in 'ACGT': search(motif + base, depth + 1) start_time = time.time() search('', 0) end_time = time.time() return best_motif, end_time - start_time # 遍历法
def find_motif_brute_force(sequences, motif_length): best_motif = None best_score = float('inf') start_time = time.time() for motif in itertools.product('ACGT', repeat=motif_length): motif = ''.join(motif) score = calculate_score(motif, sequences) if score < best_score: best_score = score best_motif = motif end_time = time.time() return best_motif, end_time - start_time # 示例序列集合和motif长度
sequences = ['ACGTAGCTAG', 'ACGGATCGTA', 'TAGCTAGCTA', 'TCGATCGATT']
motif_length = 3 # 使用分支界定法寻找motif
motif_bb, time_bb = find_motif_branch_bound(sequences, motif_length)
print(f"Branch and Bound Motif: {motif_bb}, Time: {time_bb:.4f}s") # 使用遍历法寻找motif
motif_bf, time_bf = find_motif_brute_force(sequences, motif_length)
print(f"Brute Force Motif: {motif_bf}, Time: {time_bf:.4f}s")