【2023年中国高校大数据挑战赛 】赛题 B DNA 存储中的序列聚类与比对 Python实现

【2023年中国高校大数据挑战赛 】赛题 B DNA 存储中的序列聚类与比对 Python实现

更新时间:2023-12-29

1 题目

赛题 B DNA 存储中的序列聚类与比对

近年来,随着新互联网设备的大量涌入和对其服务需求的指数级增长,越来越多的数据信息被产生与收集。预计到 2021 年,数据中心内部的IP流量将达到14.7 ZB,数据中心之间的流量将达到 2.8 ZB。如何储存与运输如此庞大的数据已经成为了难题。DNA存储技术是一项着眼于未来的具有划时代意义存储技术,正成为应对数据爆炸的关键技术之一。DNA存储技术指的是使用人工合成的脱氧核糖核苷酸(DNA)作为介质进行信息存储的技术,其具有理论存储量大、维护方便的优点。具体来说,DNA存储将计算机的二进制信息转换为四种碱基(腺嘌呤A、胸腺嘧啶T、鸟嘌呤G和胞嘧啶C)组成的DNA序列(相当于转换为四进制),之后合成为DNA分子干粉。需要读取信息时,将DNA分子进行PCR扩增(这步将会使得原有DNA序列进行扩增复制),之后使用测序仪测出DNA信息。然而在合成、测序等阶段会存在一定的错误,有概率随机发生碱基删除、增添或者替换。下图是某个序列合成测序后的示意图,可以看出由于发生了碱基删除、增添和替换,进而将ATGCATGC变成了AGCAATTC:

在这里插入图片描述

因此,对于我们设计好的DNA序列,实际生产测序出来后的序列会存在以下差异:

  • 测序后的序列将比原始序列的数量多很多,因为原始序列会被随机扩增成很多条。

  • 测序后的序列相比于原始序列有可能存在错误,包括某个碱基缺失、替换、或添加了某个未知碱基,甚至会出现断链。

针对以上两个特点,目前往往需要对测序后的序列进行聚类与比对。其中聚类指的是将测序序列聚类以判断原始序列有多少条,聚类后相同类的序列定义为一个簇。比对则是指在聚类基础上对一个簇内的序列进行比对进而输出一条最有可能的正确序列。通过聚类与比对将会极大地恢复原始序列的信息,但需要注意由于DNA测序后序列众多,如何高效地进行聚类与比对则是在满足准确率基础上的另一大难点。

“train_reference.txt”是某次合成的目标序列,其中第一行为序号,第二行为序列内容。通过真实合成、测序后读取到的测序序列文件为“train_reads.txt”,我们已经对测序序列进行了分类,该文件第一行为目标序列的序号,第二行为序列内容。

基于赛题提供的数据,自主查阅资料,选择合适的方法完成如下任务:

**任务 1:**观察数据集“train_reads.txt”、“train_reference.txt”,针对这次合成任务,进行错误率(插入、删除、替换、断链)、拷贝数方面的分析。其中错误率定义为某个碱基发生错误的概率,需要对不同类型的错误率分别进行分析。拷贝数定义为原始序列复制的数量。

**任务 2:**设计开发一种模型用于对测序后的序列“train_reads.txt”进行聚类,并根据“train_reads.txt”的标签验证模型准确性。模型主要从两方面评估效果:

(1)聚类后准确性(包括簇的数量以及簇内纯度)、(2)聚类速度(以分钟为单位)。

任务 3: “test_reads.txt”是我们在另一种合成环境下合成的测序文件(与 “train_reads.txt”的目标序列不相同),请用任务 2 所开发的模型对其进行聚类,给出聚类耗时以及“test_reads.txt”的目标序列数量,给出拷贝数分布图。

任务 4: 聚类后能否通过比对恢复原始信息也是极为关键的,设计开发一种用于同簇序列的比对模型,该模型可以针对同簇的DNA序列进行比对并输出最有可能正确的目标序列。 请使用该工具对任务 3 中“test_reads.txt”的聚类后序列进行比对,并输出“test_reads.txt”最有可能的目标序列,并分析“test_reads.txt”的错误率。(请用一个“test_ref.txt”的文件记录“test_reads.txt”的目标序列,文件内序列的形式为:

AAAA……
AAAT……
AATA……
……
CCCC……

即序列只用回车间隔,不需要加其他符号,序列顺序按照从前到后,ATGC依次的顺序。此外,需要在论文中展示前十条目标序列的聚类结果。)

附件 1:train_reference.txt train数据集的正确序列
附件 2:train_reads.txt train数据集的合成测序后序列
附件 3:test_reads.txt test数据集的合成测序后序列

参考文献:

  • Dong Y, Sun F, Ping Z, et al. DNA storage: research landscape and future prospects[J]. National Science Review, 2020, 7(6): 1092-1107.

  • Fu L, Niu B, Zhu Z, et al. CD-HIT: accelerated for clustering the next-generation sequencing data[J]. Bioinformatics, 2012, 28(23): 3150-3152.

2 问题分析

2.1 问题一

定义一个函数来比较两个字符串序列,可以自己写for循环去比较,也可以使用字符串比较工具SequenceMatcher。

2.2 问题二

DNA序列的聚类可以采用基于字符串相似度的聚类方法,比如Levenshtein、SMITH-WATERMAN、N-gram方法、或基于序列编码(如k-mer计数)的机器学习聚类方法。

2.3 问题三

在问题二的基础上,对train_reads.xt和test_reads文件和k-mer词频矩阵进行聚类分析,以判断原始序列有多少条。统计每个簇中的序列数量,得到拷贝数分布图。

2.4 问题四

(1)同簇的DNA序列比对方法:对每个簇中的序列进行多数投票,多数序列中出现的碱基将被选为最终序列的对应位置的碱基.

(2)对于每个聚类簇,进行列方向的比对,也就是对于序列的每个位置,从属于该簇的所有序列中选取每个位置上最常出现的碱基作为该位置的最终碱基。

(3)对多数投票的结果,进一步进行相似性评分,比较每个簇的共识序列(从投票中获得的序列)与引用序列库(理想的序列)中的序列。

(4)对于找到的共识序列,将其结果按照聚类簇的索引排序并输出,以方便与目标序列文件(“test_ref.txt”)进行比对,来确定错误位置和错误率。

(5)改进角度:使用更加复杂的比对算法,例如全局比对、局部比对算法、Smith-Waterman、Needleman-Wunsch算法,这些算法考虑了插入、删除和替换,并能够为每种类型的差错提供权重。

3 Python实现

3.1 问题一

import pandas as pd
from difflib import SequenceMatcher
from collections import Counter
from pyecharts.charts import Bar, Pie
from pyecharts import options as opts# 读取目标序列文件和测序序列文件
reference_seq_s = pd.read_csv('data/train_reference.txt',sep=' ',names=['ID','DNA_ref'])
reads = pd.read_csv('data/train_reads.txt',sep=' ',names=['ID','DNA'])
merged_df = pd.merge(reference_seq_s, reads, on='ID', how='inner')# 初始化统计变量
insertion_errors = 0
deletion_errors = 0
replacement_errors = 0
chain_breaks = 0
copy_numbers = Counter()# 定义一个函数来比较两个序列,并统计不同类型的错误
def analyze_sequence(ref_seq, test_seq):global insertion_errors, deletion_errors, replacement_errors, chain_breaks# 略for tag, i1, i2, j1, j2 in s.get_opcodes():if tag == 'replace':replacement_errors += max(i2 - i1, j2 - j1)elif tag == 'delete':deletion_errors += (i2 - i1)elif tag == 'insert':insertion_errors += (j2 - j1)elif tag == 'equal':pass  # No errorif len(ref_seq) != len(test_seq):chain_breaks += 1# 进行错误统计和拷贝数计算
for index, row in merged_df.iterrows():analyze_sequence(row['DNA_ref'], row['DNA'])copy_numbers[row['ID']] += 1

# 总的测序次数
total_reads = len(merged_df)# 绘制错误率和拷贝数统计图
def create_charts():# 错误率统计图error_bar = (Bar(init_opts=opts.InitOpts(width="700px", height="500px")).add_xaxis(['Insertion', 'Deletion', 'Replacement', 'Chain Breaks']).add_yaxis('Errors', [insertion_errors, deletion_errors, replacement_errors, chain_breaks]).set_global_opts(title_opts=opts.TitleOpts(title="DNA Sequence Errors")))# 拷贝数统计图copy_num_pie = (Pie(init_opts=opts.InitOpts(width="700px", height="500px")).add("",[list(z) for z in zip([str(id) for id in copy_numbers.keys()], copy_numbers.values())],radius=["40%", "75%"],).set_global_opts(title_opts=opts.TitleOpts(title="DNA Sequence Copy Numbers"),legend_opts=opts.LegendOpts(orient="vertical", pos_top="15%", pos_left="2%"),).set_series_opts(label_opts=opts.LabelOpts(formatter="{b}: {c}")))return error_bar, copy_num_pie# 创建和渲染图表
error_bar, copy_num_pie = create_charts()
error_bar.render("breakdown_of_errors.html")
copy_num_pie.render("dna_copy_numbers.html")

在这里插入图片描述
在这里插入图片描述

3.2 问题二

方法一:基于Levenshtein距离的聚类算法

import pandas as pd
from sklearn.cluster import AgglomerativeClustering
import Levenshtein
import time# 读取数据
reference_seq_s = pd.read_csv('data/train_reference.txt', sep=' ', names=['ID', 'DNA_ref'])
reads = pd.read_csv('data/train_reads.txt', sep=' ', names=['ID', 'DNA'])# 计算Levenshtein距离矩阵(由于计算量大,这里只示范计算前n个序列的距离矩阵)
n = len(reads)
distance_matrix = [[0] * n for _ in range(n)]
for i in range(n):for j in range(i+1, n):略。。。# 聚类
start_time = time.time()
clustering_model = AgglomerativeClustering(affinity='precomputed', linkage='complete', n_clusters=None, distance_threshold=1.0)
clustering_model.fit(distance_matrix)
duration = time.time() - start_time# 评估聚类结果,这里计算不同簇的数量
clusters = clustering_model.labels_
cluster_counts = pd.Series(clusters).value_counts()import numpy as np
import matplotlib.pyplot as plt
from scipy.cluster.hierarchy import dendrogram# 画出树状图
def plot_dendrogram(model, **kwargs):children = model.children_distance = np.arange(children.shape[0])no_of_observations = np.arange(2, children.shape[0]+2)linkage_matrix = np.column_stack([children, distance, no_of_observations]).astype(float)dendrogram(linkage_matrix, **kwargs)plt.figure(figsize=(15, 8))
plot_dendrogram(clustering_model, labels=range(len(reads)))
plt.ylabel("Distance")
plt.savefig('img/层次聚类.png',dpi=100)
plt.show()

在这里插入图片描述
方法二:基于SMITH-WATERMAN算法的聚类


import pandas as pd
import numpy as np
from sklearn.cluster import AgglomerativeClustering
import matplotlib.pyplot as plt
import itertools
# from Bio import pairwise2# 数据读取
reference_seq_s = pd.read_csv('data/train_reference.txt', sep=' ', names=['ID','DNA_ref'])
reads = pd.read_csv('data/train_reads.txt', sep=' ', names=['ID','DNA'])# SMITH-WATERMAN算法的实现
def smith_waterman_alignment(s1, s2, match_score=3, gap_cost=2):# 初始化得分矩阵A = np.zeros((len(s1) + 1, len(s2) + 1), int)for i, j in itertools.product(range(1, A.shape[0]), range(1, A.shape[1])):match = A[i - 1, j - 1] + (match_score if s1[i - 1] == s2[j - 1] else -match_score)delete = A[i - 1, j] - gap_costinsert = A[i, j - 1] - gap_costA[i, j] = max(match, delete, insert, 0)return np.max(A)# 编辑距离矩阵的计算
def compute_distance_matrix(reads):n_reads = len(reads)distance_matrix = np.zeros((n_reads, n_reads))for i in range(n_reads):for j in range(i+1, n_reads):alignment_score = smith_waterman_alignment(reads[i], reads[j])distance_matrix[i, j] = distance_matrix[j, i] = alignment_score # we use alignment score directly herereturn distance_matrix# Run SMITH-WATERMAN on the dataset
distance_matrix = compute_distance_matrix(reads['DNA'].values)# 聚类算法
def cluster_sequences(distance_matrix, n_clusters=2):# 使用层次聚类,可以使用其他聚类算法clustering = AgglomerativeClustering(n_clusters=n_clusters, affinity='precomputed', linkage='complete')# 使用 1 减 距离矩阵,是为了将距离转化为相似度clustering.fit(1 - distance_matrix)return clustering.labels_
# 聚类和评估
cluster_labels = cluster_sequences(distance_matrix)
reads['Cluster'] = cluster_labels# 评估簇的纯度
def evaluate_cluster_purity(cluster_labels, actual_labels):contingency_table = pd.crosstab(cluster_labels, actual_labels)purity = np.sum(np.max(contingency_table, axis=0)) / np.sum(contingency_table.sum())return purity# 可视化
def visualize_clustering(reads, cluster_labels):plt.figure(figsize=(12, 8))colors = ['r', 'g', 'b', 'y', 'c', 'm']for i in np.unique(cluster_labels):plt.plot(reads[reads['Cluster'] == i]['DNA'].index, [i] * sum(reads['Cluster'] == i), 'x', color=colors[i % len(colors)], label=f'Cluster {i}')plt.title('Clustering of DNA sequences')plt.xlabel('Sequence Index')plt.ylabel('Cluster ID')plt.legend()plt.show()visualize_clustering(reads, cluster_labels)

方法三:对测序序列进行k-mer编码。使用CountVectorizer把序列的k-mer列表转换成词频(term frequency)矩阵。使用K-means算法对k-mer词频矩阵进行聚类,聚类数设置为原始序列数。

import pandas as pd
import numpy as np
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.feature_extraction.text import CountVectorizer
import matplotlib.pyplot as plt
import time# 读取数据
reference_seq_s = pd.read_csv('data/train_reference.txt', sep=' ', names=['ID', 'DNA_ref'])
reads = pd.read_csv('data/train_reads.txt', sep=' ', names=['ID', 'DNA'])# k-mer计数函数
def get_kmers(sequence, k=3):return [sequence[x:x+k] for x in range(len(sequence) + 1 - k)]reads['kmers'] = reads['DNA'].apply(lambda x: get_kmers(x))# 将k-mer列表转换为字符串(以便进一步转换为向量)
reads['kmers_str'] = reads['kmers'].apply(lambda x: ' '.join(x))# 使用CountVectorizer生成k-mer的词频矩阵
vectorizer = CountVectorizer()
X = vectorizer.fit_transform(reads['kmers_str'])# PCA降维
pca = PCA(n_components=2)
X_reduced = pca.fit_transform(X.toarray())# KMeans聚类
# 确定簇的数量为原始序列数
n_clusters = len(reference_seq_s['ID'].unique())
kmeans = KMeans(n_clusters=n_clusters)start_time = time.time()# 训练模型
kmeans.fit(X)
end_time = time.time()# 计算总耗时
total_time = (end_time - start_time) / 60
print("聚类时间{:.2f} minutes.".format(total_time))labels = kmeans.labels_
reads['cluster'] = labels
# 可视化结果
plt.figure(figsize=(8, 6))
plt.scatter(X_reduced[:, 0], X_reduced[:, 1], c=labels, cmap='rainbow', alpha=0.6, edgecolors='w', s=50)
plt.savefig('img/k-cluster.png',dpi=100)
plt.show()

在这里插入图片描述

3.3 问题三

import pandas as pd
from sklearn.cluster import KMeans
from sklearn.feature_extraction.text import CountVectorizer
import pyecharts.options as opts
from pyecharts.charts import Bar
import time# k-mer计数函数
def get_kmers(sequence, k=3):return [sequence[x:x+k] for x in range(len(sequence) + 1 - k)]
# 读取数据
# reference_seq_s = pd.read_csv('data/train_reference.txt', sep=' ', names=['ID', 'DNA'])
reads = pd.read_csv('data/train_reads.txt', sep=' ', names=['ID', 'DNA'])
test_reads = pd.read_csv('data/test_reads.txt',names=['DNA'])reads['kmers'] = reads['DNA'].apply(lambda x: get_kmers(x))
# 将k-mer列表转换为字符串(以便进一步转换为向量)
reads['kmers_str'] = reads['kmers'].apply(lambda x: ' '.join(x))
# 应用k-mer处理
test_reads['kmers'] = test_reads['DNA'].apply(lambda x: get_kmers(x))
test_reads['kmers_str'] = test_reads['kmers'].apply(lambda x: ' '.join(x))# 使用CountVectorizer生成k-mer的词频矩阵
vectorizer = CountVectorizer()
# 先拟合训练数据
X_train = vectorizer.fit_transform(reads['kmers_str'])
# 再转换测试数据
X_test = vectorizer.transform(test_reads['kmers_str'])from sklearn.decomposition import PCA
# 用PCA降维以便可视化(仅用于降维和可视化,并不用于聚类)
pca = PCA(n_components=2)
X_reduced = pca.fit_transform(X_train.toarray())
# KMeans聚类
start_time = time.time()
n_clusters = len(reads['ID'].unique())  
kmeans = KMeans(n_clusters=n_clusters)
# 训练模型
kmeans.fit(X_train)
clusters = kmeans.fit_predict(X_test)
end_time = time.time()
# 输出聚类耗时
print(f"Clustering Time: {end_time - start_time}")# 统计每个簇的拷贝数
cluster_counts = pd.Series(clusters).value_counts().sort_index()

在这里插入图片描述
在这里插入图片描述

3.4 问题四

(1)方法一

from sklearn.decomposition import PCA
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.feature_extraction.text import CountVectorizer
import time# k-mer计数函数
def get_kmers(sequence, k=3):return [sequence[x:x+k] for x in range(len(sequence) + 1 - k)]
# 读取数据
# reference_seq_s = pd.read_csv('data/train_reference.txt', sep=' ', names=['ID', 'DNA'])
reads = pd.read_csv('data/train_reads.txt', sep=' ', names=['ID', 'DNA'])
test_reads = pd.read_csv('data/test_reads.txt',names=['DNA'])reads['kmers'] = reads['DNA'].apply(lambda x: get_kmers(x))
# 将k-mer列表转换为字符串(以便进一步转换为向量)
reads['kmers_str'] = reads['kmers'].apply(lambda x: ' '.join(x))
# 应用k-mer处理
test_reads['kmers'] = test_reads['DNA'].apply(lambda x: get_kmers(x))
test_reads['kmers_str'] = test_reads['kmers'].apply(lambda x: ' '.join(x))
# 使用CountVectorizer生成k-mer的词频矩阵
vectorizer = CountVectorizer()
# 先拟合训练数据
X_train = vectorizer.fit_transform(reads['kmers_str'])
# 再转换测试数据
X_test = vectorizer.transform(test_reads['kmers_str'])
# 用PCA降维以便可视化(仅用于降维和可视化,并不用于聚类)
pca = PCA(n_components=2)
X_reduced = pca.fit_transform(X_train.toarray())
# KMeans聚类
start_time = time.time()
n_clusters = len(reads['ID'].unique())  
kmeans = KMeans(n_clusters=n_clusters)
# 训练模型
kmeans.fit(X_train)
clusters = kmeans.fit_predict(X_test)# 比对模型的Python代码实现
import numpy as np
from collections import Counter
from typing import List# 函数来计算多数投票后确定的序列
def consensus_sequence(seqs: List[str]) -> str:"""采取多数投票法,返回一个列表中最可能正确的目标序列。:param seqs: 需要进行多数投票的一系列序列:return: 最可能正确的目标序列"""# 将序列转置,以方便进行列方向投票transposed_seqs = list(zip(*seqs))consensus_seq = []# 对于每个位置,计算最常见的元素for column in transposed_seqs:counter = Counter(column)most_common = counter.most_common(1)[0][0]consensus_seq.append(most_common)return ''.join(consensus_seq)# 根据聚类结果对序列进行聚类
clustered_seqs = {}  # 存储每个原始序列ID对应的所有序列
# 对测试数据聚类
for idx, cluster_id in enumerate(clusters):if cluster_id not in clustered_seqs:clustered_seqs[cluster_id] = []clustered_seqs[cluster_id].append(test_reads['DNA'][idx])# 对于每个聚类,进行比对,并确定共识序列
consensus_seqs = {}
for cluster_id, seqs in clustered_seqs.items():consensus = consensus_sequence(seqs)consensus_seqs[cluster_id] = consensus
# 评估聚类质量和恢复的序列质量
reference_seqs = pd.read_csv('data/train_reference.txt', sep=' ', names=['ID', 'DNA'])# 评估聚类质量和恢复的序列质量
reference_seqs = pd.read_csv('data/train_reference.txt', sep=' ', names=['ID', 'DNA'])
# 计算共识序列与目标序列的错误率
def calculate_error_rate(original_seq: str, new_seq: str) -> float:"""计算恢复的序列与目标序列之间的错误率。:param original_seq: 原序列:param new_seq: 恢复的序列:return: 错误率"""errors = sum(1 for orig, new in zip(original_seq, new_seq) if orig != new)return errors / len(original_seq)# 错误率列表
error_rates = []
# 输出最可能正确的序列并计算错误率
for cluster_id, cons_seq in sorted(consensus_seqs.items()):original_seq = reference_seqs.loc[cluster_id,'DNA']error_rate = calculate_error_rate(original_seq, cons_seq)error_rates.append(error_rate)print(f"Cluster {cluster_id} Consensus: {cons_seq}, Error Rate: {error_rate}")# 分析总体错误率
overall_error_rate = np.mean(error_rates)
print(f"总体错误率: {overall_error_rate}")

总体错误率:0.509

(2)方法二

from sklearn.decomposition import PCA
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.feature_extraction.text import CountVectorizer
import time# k-mer计数函数
def get_kmers(sequence, k=3):return [sequence[x:x+k] for x in range(len(sequence) + 1 - k)]
# 读取数据
# reference_seq_s = pd.read_csv('data/train_reference.txt', sep=' ', names=['ID', 'DNA'])
reads = pd.read_csv('data/train_reads.txt', sep=' ', names=['ID', 'DNA'])
test_reads = pd.read_csv('data/test_reads.txt',names=['DNA'])reads['kmers'] = reads['DNA'].apply(lambda x: get_kmers(x))
# 将k-mer列表转换为字符串(以便进一步转换为向量)
reads['kmers_str'] = reads['kmers'].apply(lambda x: ' '.join(x))
# 应用k-mer处理
test_reads['kmers'] = test_reads['DNA'].apply(lambda x: get_kmers(x))
test_reads['kmers_str'] = test_reads['kmers'].apply(lambda x: ' '.join(x))
# 使用CountVectorizer生成k-mer的词频矩阵
vectorizer = CountVectorizer()
# 先拟合训练数据
X_train = vectorizer.fit_transform(reads['kmers_str'])
# 再转换测试数据
X_test = vectorizer.transform(test_reads['kmers_str'])
# 用PCA降维以便可视化(仅用于降维和可视化,并不用于聚类)
pca = PCA(n_components=2)
X_reduced = pca.fit_transform(X_train.toarray())
# KMeans聚类
start_time = time.time()
n_clusters = len(reads['ID'].unique())  
kmeans = KMeans(n_clusters=n_clusters)
# 训练模型
kmeans.fit(X_train)
clusters = kmeans.fit_predict(X_test)import numpy as np
import pandas as pd
from collections import Counter# Needleman-Wunsch算法实现
def needleman_wunsch(seq1, seq2, match_score=1, gap_cost=1, mismatch_cost=1):n = len(seq1)m = len(seq2)score_matrix = np.zeros((n+1, m+1))# Initialize score matrix and traceback pathsfor i in range(n+1):score_matrix[i][0] = -i * gap_costfor j in range(m+1):score_matrix[0][j] = -j * gap_cost# Fill in score matrixfor i in range(1, n+1):for j in range(1, m+1):if seq1[i-1] == seq2[j-1]:match = score_matrix[i-1][j-1] + match_scoreelse:match = score_matrix[i-1][j-1] - mismatch_costdelete = score_matrix[i-1][j] - gap_costinsert = score_matrix[i][j-1] - gap_costscore_matrix[i][j] = max(match, delete, insert)# Traceback to compute the alignmentalign1 = ""align2 = ""i = nj = mwhile i > 0 and j > 0:score_current = score_matrix[i][j]score_diagonal = score_matrix[i-1][j-1]score_up = score_matrix[i][j-1]score_left = score_matrix[i-1][j]if score_current == score_diagonal + (match_score if seq1[i-1] == seq2[j-1] else -mismatch_cost):align1 += seq1[i-1]align2 += seq2[j-1]i -= 1j -= 1elif score_current == score_left - gap_cost:align1 += seq1[i-1]align2 += "-"i -= 1elif score_current == score_up - gap_cost:align1 += "-"align2 += seq2[j-1]j -= 1while i > 0:align1 += seq1[i-1]align2 += "-"i -= 1while j > 0:align1 += "-"align2 += seq2[j-1]j -= 1return align1[::-1], align2[::-1]# 从聚类结果中恢复出最可能的序列
def recover_sequence(cluster_seqs):# 序列长度可能不同,先找到最长的序列长度return consensus_sequencefrom functools import reduce
# 使用先前完成的KMeans结果clusters
# 假设clusters为序列的聚类结果,test_reads为相应的序列数据
cluster_dict = {i: [] for i in range(n_clusters)}
for i, cluster in enumerate(clusters):cluster_dict[cluster].append(test_reads['DNA'][i])# 对每个簇进行比对,并且输出最可能正确的序列
consensus_sequences = {}
for cluster_id, seqs in cluster_dict.items():if len(seqs) > 1:# 使用reduce函数将同簇序列两两比对consensus = reduce(lambda x, y: recover_sequence([x, y]), seqs)else:# 如果簇内只有一个序列,则将其作为最可能的序列consensus = seqs[0]consensus_sequences[cluster_id] = consensus# 将得到的“最可能正确的序列”写入到文件
with open('data/test_ref.txt', 'w') as f_out:for seq in consensus_sequences.values():f_out.write(seq + '\n')

在这里插入图片描述

4 完整代码

请看名片扣我

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/230289.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

【日积月累】Java Lambda 表达式

目录 【日积月累】Java Lambda 表达式 1.前言2.语法3.应用场景3.1简化匿名内部类的编写3.1简化匿名内部类的编写3.2简化集合类中的操作3.3实现函数式接口3.4简化多个方法的调用3.5简化异步编程 4.总结5.参考 文章所属专区 日积月累 1.前言 Lambda表达式是一个匿名函数&#…

快速打通 Vue 3(二):响应式对象基础

很激动进入了 Vue 3 的学习,作为一个已经上线了三年多的框架,很多项目都开始使用 Vue 3 来编写了 这一组文章主要聚焦于 Vue 3 的新技术和新特性 如果想要学习基础的 Vue 语法可以看我专栏中的其他博客 Vue(一):Vue 入…

2023年.AI域名销售额达550万美元 2024还要继续涨

根据域名投资专家Elliot Silver的最新文章,2023年公开报道的.AI域名销售额已经达到了550万美元,而2022年和2021年分别为90万美元和120万美元。 Silver观察到过去几年.AI域名销售额呈现逐年增长的趋势,尤其是2023年的销售额相较前两年有了显著…

【计算机毕业设计】SSM二手交易网站

项目介绍 该项目分为前后台,前台普通用户角色,后台管理员角色。 管理员主要功能如下: 登陆,商品分类管理,商品管理,商品订单管理,用户管理等功能。 用户角色主要功能如下: 包含以下功能:查看所有商品,用户登陆注册…

骨传导耳机不踩坑推荐指南,南卡/韶音/墨觉实测告诉你答案!

你知道怎么选骨传导耳机吗?作为一个音响测评博主,我用过不下10款骨传导耳机,有的是普通款式,有的是专业运动款式,甚至为了维修也拆过一些骨传导耳机。可以说,骨传导耳机的选购绝不是表面看起来那么简单&…

单位转换工具类

单位转换工具类 1. 工具类转换- 定义装换枚举转换类型- 创建转换工具类,1. 通过反射去除字段,2.对照传入map标记的字段需要转换的类型转换3. 重新赋值 2. 注解转换- 定义注解- 解析注解 1. 工具类转换 - 定义装换枚举转换类型 public enum UnitConvertType {/*** 精确度*/ACC…

break,continue跳出指定循环小案例

某一天,你犯了一个错误,你老婆罚你做5天家务,每天去洗碗,洗碗到第三天心软了,原谅你了只有第三太不用洗碗 public class BreakDemo {public static void main(String[] args) {//某一天,你犯了一个错误&am…

视频监控可视化云平台EasyCVR智能视频技术优势分析

TSINGSEE青犀视频安防视频管理系统EasyCVR视频智能融合共享平台,是一个支持Windows/Linux(CentOS ubuntu)/国产化系统的视频管理平台。平台可以支持多协议接入,通过视频应用引擎将多种格式的视频数据转换为统一的视频流数据,支持无插件H5直播…

现在的人们如何看待数据隐私?

PrimiHub一款由密码学专家团队打造的开源隐私计算平台,专注于分享数据安全、密码学、联邦学习、同态加密等隐私计算领域的技术和内容。 在当前时代,每一次点击、触摸或按键都留下了数字痕迹。但是我们对自己的个人数据几乎没有控制的权限,这让…

主流桌面浏览器Chrome,FireFox和Edge等如何禁用弹出式窗口阻止程序,这里有详细步骤

为什么你想知道如何禁用浏览器中的弹出式窗口阻止程序?毕竟,弹出式窗口是网络的祸害:显示烦人的广告、虚假的安全消息和其他刺激,会分散你的浏览注意力,甚至可能包含恶意代码。 所有主要的桌面浏览器现在都默认阻止弹出式窗口,那么你到底为什么要取消阻止这些害虫呢?事…

c++的三大特性之关于继承

目录 继承的概念及定义 基类和派生类对象赋值转换 继承中的作用域 派生类的默认成员函数 继承与友元,静态成员 继承的概念及定义 概念: 继承(inheritance)机制是面向对象程序设计使代码可以复用的最重要的手段,它允许程序员在保持原有类…

gitee创建仓库

描述 本文章记录了怎么在gitee上创建项目,以及使用vscode提代码到远程呢个仓库,如何创建一个新分支,并将新分支提交到远程仓库。 1、创建远程仓库 在创建远程仓库之前要先进行ssh密钥的设置 (1)打开黑窗口&#xff…

iptables 防火墙(二)

目录 1. SNAT 策略及应用 1.1 SNAT策略概述 1. 只开启路由转发,未设置地址转换的情况 2. 开启路由转发,并设置SNAT转换的情况 1.2 SNAT策略的应用 1. 2.1 共享固定IP上网 (1)打开网关的路由转发 (2)…

[C#]C# OpenVINO部署yolov8目标检测模型

【官方框架地址】 https://github.com/ultralytics/ultralytics.git 【算法介绍】 YOLOv8 抛弃了前几代模型的 Anchor-Base。 YOLO 是一种基于图像全局信息进行预测的目标检测系统。自 2015 年 Joseph Redmon、Ali Farhadi 等人提出初代模型以来,领域内的研究者们…

大数据Doris(四十九):Doris数据导出介绍

文章目录 Doris数据导出介绍 一、​​​​​​​使用示例

静态网页设计——个人简介网站

前言 使用经典前端三件套HTMLCSSJavascript编写了一个关于个人简介的静态网页,可以根据自己的需要,十分简单的进行修改。 首页 首页由上方的菜单栏以及菜单栏下面的轮播图组成,再往下走,是关于自己的兴趣爱好的部分&#xff0c…

django学习:页面渲染与请求和响应

1.请求过程 2.页面渲染 在app中新建一个目录(Directory),文件名命名为templates。该文件名命名是固定的,不可命名出错,如若后续步骤出错,该目录文件名是一个检查的重点项目。在该目录下新建一个html文件&a…

深度学习:大规模模型分布式训练框架DeepSpeed

深度学习:大规模模型分布式训练框架DeepSpeed DeepSpeed简介DeepSpeed核心特点DeepSpeed如何工作?DeepSpeed如何使用?参考文献 DeepSpeed简介 随着机器学习模型变得越来越复杂和庞大,训练这些模型所需的计算资源也在不断增加。特别…

如何查找iPhone中所有的应用程序

​ ​ Apple 的 App Store 共有约 200 万个适用于 iPhone 和 iPad 的应用程序。如果您像我们一样,您的 iOS 或 iPadOS 设备上可能有数十个应用程序,但没有机会将它们全部整理好。您很容易忘记主屏幕上应用程序图标的位置。 幸运的是,iPhone…

向日葵远程工具的使用Mysql5.7的安装与配置

目录 一、向日葵远程安装与使用 二、Mysql 5.7 安装与配置 2.1 安装 2.2 Navicat Premium 12 测试连接 本机测试连接 外部访问MySQL测试连接 三、思维导图 一、向日葵远程安装与使用 简介: 向日葵远程控制是一款用于对远程PC进行管理和服务的软件,拥有5秒快速…