【生物信息】ecDNA染色体外环状DNA的识别——AmpliconArchitect算法

写在前面

这个项目是我大三时期（2020年）做的SRTP，是与我们学校的生物信息专业的同学一起做的。生物上的东西我不是很了解，我只是帮我的同学找了一些工具并使用这些工具完成我们的SRTP。
SRTP题目是：染色体外DNA的识别与分析
我们主要是参考了这篇2019年发表在《nature》上的论文：Circular ecDNA promotes accessible chromatin and high oncogene expression
他们论文中提到了他们使用的AA算法，并提供了AA算法项目的开源地址，有了AA算法，我们就可以通过基因序列找到环状DNA，并把这些基因序列组合成环，输出成图片。
以下是我们在某个癌症细胞基因序列上提取到的ecDNA扩增子的重建结果:

实验环境

操作系统：deepin15.5 sp2
python: python3.8

原理简介

算法流程

AA采取许多步骤来预测扩增子（amplicons）的结构：
step1:间隔集确定：确定要重建的每个扩增子的间隔列表。
step2:SV检测：利用每个扩增子中的覆盖率和不一致的读段双簇（coverage and discordant read pairs）来检测拷贝数变化和结构变异。
step3:断点图构造：构造一个断点图，该断点图由序列egdes（基因组片段），断点边（连接的因组位置对）和源顶点（可选）组成，并预测所有边的拷贝数。
step4:循环分解：将断点图分解为简单的循环，以简单的方式表示预测的扩增子（amplicons）结构。
step5:交互式循环合并：提供一个界面，以交互方式合并或修改环视图（cycles）以探索候选结构。

实现细节

AmpliconArchitect实施了一系列步骤，从种子区间开始，最终重建了扩增子的完整结构，并为每个阶段提供了有益的结果：

区间搜索

在此步骤中，AmpliconArchitect从种子区间开始，迭代地标识属于扩增子的区间列表。首先创建一个存储种子区间的最大堆数据结构。
AmpliconArchitect重复以下步骤，直到最大堆为空或经过10次迭代。
在每次迭代中，AmpliconArchitect选择一个区间并以拷贝数敏感的方式确定不一致成对读段双簇，并选择双簇，其配偶读段比对位置位于先前看到的区间之外。然后，通过查询扩展部分是否被扩增，尝试扩展双簇。如果查询段的窗口（10kbp）至少有20％且覆盖深度>，或者小于20kbp且包含至少两个不一致边（双簇大小对应于拷贝数=2），则分类为扩增。然后，AmpliconArchitect迭代地将扩展部分的大小加倍，直到发现该扩展被扩增，然后通过迭代将扩展查询的大小减小一半，来有效地扩展查询的两部分。如果AmpliconArchitect能够成功扩展查询的双簇，则它将其进一步扩展100kbp，并记录扩展区间以用于将来的迭代。在AmpliconArchitect记录了所有放大的邻近区间之后，AmpliconArchitect标记所看到的区间，根据连接到先前看到的区间的不一致成对读段的数量更新了最大堆，并贪婪地选择堆顶部的区间以进行下一次迭代。AmpliconArchitect报告延伸步骤中所有扩增的区间。

区间重排，分区和可视化

AmpliconArchitect计算了所有覆盖范围的均值移位边界和相应段的初始拷贝数估计值。然后，它创建了不一致成对读段双簇列表，其双簇大小阈值由拷贝数估计值（或差异）确定。此外，对于没有找到匹配的不一致读段双簇的均值移位边界，它对双簇大小阈值仅为两个读段对的不一致成对读段执行敏感本地搜索。最终，它创建了所有重排的基因组位置集，并带有不一致成对读段双簇或均值移位边界（没有匹配的不一致成对读段）。使用这组位置，它将所有区间划分为序列边。对于第二阶段的输出，AmpliconArchitect创建了一个单独的图，称为结构变异视图，该图显示了带有覆盖深度直方图的区间集，均值移位段的初始拷贝数估计以及不一致成对读段双簇。

断点图的构造

AmpliconArchitect使用序列边分区构造了一个断点图。对于每个序列边，它创建了两个断点顶点，标记了基因组片段的起点和终点，并添加了连接两个顶点的序列边。AmpliconArchitect用特殊的源顶点扩展了顶点集。对于每个不一致成对读段双簇，如果两个簇都属于区间集，则会添加一个不一致断点边，该断点边连接相应序列边的各个端点。对于在区间集之外具有一个簇的所有双簇，它引入了一个源断点边，该边将源顶点连接到与扩增子区间内的簇相对应的断点顶点。它还添加了源断点边，这些源断点边将源顶点连接到与均值移位顶点相对应的断点顶点（没有对应的不一致双簇）和扩增子区间的端点。最后，AmpliconArchitect将对应于每个区间内连续序列边的断点对与断点边保持一致。对于每个序列边和断点边，AmpliconArchitect记录分别映射到该边的读段和读段对的数量。正确重建的断点图表示所有扩增子结构的叠加。每个环结构形成一个交替序列和断点边的环。端点与环外的基因组位置相连的线性结构可以表示为在源顶点处开始和结束的交替闭合步态。断点图的构建可能并不总是完整的，因为缺少边会导致结构的最终预测不准确。如果AmpliconArchitect无法检测到不一致断点边，但通过均值移位边界检测到重排的位置，则相应的扩增子结构将被表示为从源顶点开始和结束的一个或多个游走，从而部分缓解了此问题。此后，我们将断点图中的环定义限制为具有交替的序列边和断点边的闭合遍历，但该遍历可能包含最多两个连接到源的连续断点边。
使用平衡流优化的拷贝数估计：接下来，我们注意到每个边的拷贝数是每个扩增子结构的拷贝数的总和，其中边的每次遍历分别进行计数。结果，图中的拷贝数遵循平衡流动特性，其中序列边的拷贝数与连接到序列边的每个断点顶点的断点边的拷贝数的总和相匹配。AmpliconArchitect使用由拷贝数、边长度和序列覆盖深度参数确定的参数，将比对到每个边的读段片段的数量建模为Poisson分布。使用凸优化软件包Mosek获得平衡流的最佳解决方案。对于第三阶段的输出，AmpliconArchitect将图边及其拷贝计数报告为文本输出。

环分解

如上所述，线性或环扩增子结构可以在断点图中表示为一个或多个环。但是，即使使用正确的重构和对断点图的拷贝数分配，也无法始终明确地推断出环，尤其是在反复遍历边的情况下。相反，可能有两个或更多可能的环集和相关联的拷贝数，这样，将每个集合内的环组合在一起，最终可以得到具有相同拷贝数分配的相同断点图。在这里，环的组合仅意味着对每个环中每个图形边的拷贝数求和。枚举所有可能的扩增子结构并不总是可行的，因为可能的结构数目可能成倍增加。为