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

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 问题二

2.3 问题三

2.4 问题四

3 Python实现

3.1 问题一

python">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 error
    if 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

python">
# 总的测序次数
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")