2023年高校大数据挑战赛B题DNA存储中的序列聚类与比对求解全过程文档及程序

2023年高校大数据挑战赛

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

原题再现：

  近年来，随着新互联网设备的大量涌入和对其服务需求的指数级增长，越来越多的数据信息被产生与收集。预计到2021年，数据中心内部的IP流量将达到14.7 ZB，数据中心之间的流量将达到2.8ZB。如何储存与运输如此庞大的数据已经成为了难题。
  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数据集的合成测序后序列

整体求解过程概述(摘要)

  DNA存储技术以高密度和稳定性的DNA分子为基础，超越传统介质，存储容量大，长期稳定，通过转换和合成实现信息存储与读取，为应对巨大数据需求提供创新解决方案，本文通过一系列数学模型完成了基因序列聚类分析与序列对比。
  针对问题一，本文首先对原始数据进行预处理，包括分离行号和序列以及清理序列数据。随后，计算拷贝数并定义了可能出现的DNA序列的插入、删除、替换和断链错误。在这一基础上，本文采用贪心算法和BLAST算法，以确定这四类错误发生的概率，并通过可视化的方式展示了结果。进一步进行了错误类型的分析，结果显示插入错误和删除错误的发生概率约为0.5%-1%，替换错误的发生概率超过 1%，而断链错误的概率在0.1%-0.3%之间。
  针对问题2，本文首先从序列信息中提取了多个特征，包括碱基频率、碱基对频率和k-mer 频率等。随后，采用了基于粒子群算法改进的K-Means聚类方法对数据进行了聚类，同时以 Jaccard 相似度系数大于 0.75 作为验证聚类结果的指标。在评估模型时，本文综合考虑了准确率和聚类速度两个指标，通过归一化处理和方差倒数法确定了权重，以平衡这两个指标的重要性。通过模型的迭代优化，本文在迭代 50 次时取得了最佳表现，综合得分最高，准确率达到97.2%。值得一提的是，本文的模型在处理数据集时表现出了出色的效率，仅需3.05分钟完成聚类。
  针对问题三，本文在问题二模型的基础上重新定义K值，采用了肘部法则和轮廓系数相结合的方法确定最佳的K值。通过肘部法则观察SSE下降趋缓的位置，初步确定K值范围为1950-2000。然后利用轮廓系数进一步缩小选择范围，最终确定K值为1975，对应的轮廓系数高达0.93。在模型求解中，以K=1975作为聚类簇数，整个聚类过程耗时47 分钟。通过对各聚类中对应序列的拷贝数分布图进行分析，经过平滑曲线处理后得知各个聚类中心对应的基因序列数量没有明显的差异。
  针对问题四，本文利用隐马尔科夫模型对“test_reads.txt”中的基因序列进行了建模和分析。通过问题三中对每个类别的划分，为每类基因序列定义了隐藏状态，并计算了在各类别中基因序列状态转移的概率。随后，利用Baum-Welch算法学习模型参数，包括状态转移概率和观察概率，通过迭代估计获得了最终参数。接下来，应用Viterbi 算法进行序列解码，得到了所有聚类中心的目标序列。在正文中呈现了前十条目标序列的聚类结果，并对插入、删除、替换和断链错误率进行了分析。以此同时，对问题一和问题四中的四种错误发生概率进行了对比分析。结果显示，在两个数据集中，插入错误的概率接近，均为0.007。然而，在train_reads.txt数据集中，删除错误和断链错误的概率相对较高，而在test_reads.txt中，替换错误的概率较大。
  本文所构建的基因序列聚类模型不仅在实际应用中展现了较好的性能，而且从数学角度上得到了充分的合理性分析，为基因序列聚类问题提供了深刻的理论支持。

模型假设：

  1. 假设目标序列数据集和测序序列数据集的数据来源真实可靠；
  2. 假设每个测序序列至多存在一种错误类型；
  3. 假设不考虑PCR反应中可能出现的交叉污染或交叉扩增的影响。