Hierholzer 算法 —有向图求解欧拉路径(环路)

  • 欧拉路径
    如果在一张图中,可以从一点出发遍历所有的边,那么遍历过程中的这条路径就叫做欧拉路径。如果这条路径是闭合的,那就称为欧拉回路
    简单地说,如果一张图可以被“一笔画”,那么“一笔画”的那个轨迹就叫做欧拉路径。
  • 欧拉图判定定理
    含有欧拉回路的图称为欧拉图,含有欧拉路径的图称为半欧拉图。
    无向图中,如果所有顶点的度数都为偶数,则为欧拉图;如果有两个顶点的度数为奇数,其他的为偶数,则为半欧拉图。
    有向图中,如果所有顶点的入度等于出度,那么就是欧拉图。判定半欧拉图的方法也很简单,大家可以自行推理。

Hierholzer 算法

问题简述:现给出一个有向图,且为欧拉图。求欧拉回路。

Hierholzer 算法过程

  1. 选择任一顶点为起点,遍历所有相邻边。
  2. 深度搜索,访问相邻顶点。将经过的边都删除。
  3. 如果当前顶点没有相邻边,则将顶点入栈。
  4. 栈中的顶点倒序输出,就是从起点出发的欧拉回路。

从一个可能的起点出发,进行深度优先搜索,但是每次沿着辅助边从某个顶点移动到另外一个顶点的时候,都需要删除这个辅助边。如果没有可移动的路径,则将所在结点加入到栈中,并返回。

dfs(node, trace){
	while(!node.adj.isEmpty()){
		Node next = node.adj.removeLast();
		dfs(next, trace);
	}
	trace.addLast(node);
}

最后得到的栈中保存的就是整个欧拉闭迹中的顶点。(要恢复我们需要不断出栈,因此如果你用列表来存欧垃迹的话需要反转一次)。

重新安排行程:给你一份航线列表 tickets ,其中 tickets[i] = [fromi, toi] 表示飞机出发和降落的机场地点。请你对该行程进行重新规划排序。所有这些机票都属于一个从 JFK(肯尼迪国际机场)出发的先生,所以该行程必须从 JFK 开始。如果存在多种有效的行程,请你按字典排序返回最小的行程组合。例如,行程 [“JFK”, “LGA”] 与 [“JFK”, “LGB”] 相比就更小,排序更靠前。
假定所有机票至少存在一种合理的行程。且所有的机票 必须都用一次且只能用一次。

在这里插入图片描述
class Solution:
    ''' Hierholzer 算法
    '''
    def findItinerary(self, tickets):
        def dfs(cur, graph, res):
            while graph[cur]: # 遍历所有边
                dfs(graph[cur].pop(), graph, res) # 访问并删除
            if not graph[cur]: # 将顶点入栈
                res.append(cur)
        import collections
        # 建图
        graph = collections.defaultdict(list)
        for start, end in sorted(tickets)[::-1]:
            graph[start].append(end)
        res = []
        dfs("JFK", graph, res)
        return res[::-1] # 倒序输出

AlphaFold2论文

https://www.nature.com/articles/s41586-021-03819-2

https://github.com/deepmind/alphafold

沐神AlphaFold讲解

论文补充材料:https://www.biorxiv.org/content/10.1101/2021.10.04.463034v1

摘要

蛋白质对生命至关重要,了解它们的结构可以促进对其功能的系统理解。通过大量的实验,已经确定了大约 100,000 种独特的蛋白质结构,但这仅代表了数十亿已知蛋白质序列中的一小部分。

蛋白质结构是指蛋白质分子的空间结构。作为一类重要的生物大分子,蛋白质主要由化学元素组成。所有蛋白质都是由20种不同的L型α氨基酸连接形成的多聚体,在形成蛋白质后,这些氨基酸又被称为残基。
蛋白质一级结构:组成蛋白质多肽链的线性氨基酸序列。
蛋白质二级结构:依靠不同氨基酸之间的C=O和N-H基团间的氢键形成的稳定结构,主要为α螺旋和β折叠。
蛋白质三级结构:通过多个二级结构元素在三维空间的排列所形成的一个蛋白质分子的三维结构。
蛋白质四级结构:用于描述由不同多肽链(亚基)间相互作用形成具有功能的蛋白质复合物分子。

仅根据其氨基酸序列预测蛋白质的三级结构,这是“蛋白质折叠问题”, 50 多年来一直是一个重要的开放性研究问题。尽管最近取得了进展,但现有方法仍远未达到原子级准确度,尤其是当没有可用的同源结构时。

同源结构,那些不同物种因来自共同祖先而具有的相似性结构。 例如现代马经较长时间的修饰成为具有一个趾,鼹鼠及其他洞穴动物成为瘤状肢体,大象的肢体成为柱状,这样功能各异的前肢有一个共同来源,他们都来自原始陆生脊椎动物五趾型的肢体。

在这里,我们提供了第一种计算方法,即使在不知道相似结构的情况下,它也可以以原子精度定期预测蛋白质结构。我们在具有挑战性的第 14 次蛋白质结构预测关键评估 (CASP14) 中验证了我们基于神经网络模型的完全重新设计的AlphaFold,在大多数情况下表现出与实验相媲美的准确性,并且大大优于其他方法。支持最新版本的 AlphaFold 是一种新颖的机器学习方法,它将关于蛋白质结构的物理和生物学知识,利用多序列比对,融入深度学习算法的设计中。

序列比对指将两个或多个序列排列在一起,标明其相似之处。序列中可以插入间隔(通常用短横线“-”表示)。对应的相同或相似的符号(在核酸中是A, T(或U), C, G,在蛋白质中是氨基酸残基的单字母表示)排列在同一列上。
tcctctgcctctgccatcat—caaccccaaagt
|||| ||| ||||| ||||| ||||||||||||
tcctgtgcatctgcaatcatgggcaaccccaaagt
多序列比对是成对比对的延伸,是为了在一次比对里面处理多于两条的的序列。多序列比对方法试图比对一个指定序列集合里面的所有序列,这可以帮助确定这些序列的共同区段。进行多序列比对有几种方法,最常用的一种是Clustal程序集,它使用渐进多序列比对算法。Clustal在cladistics中被用来建立进化树,在PSI-BLAST和Hidden Markov model (HMM)中用来建立序列档案以在序列数据库中搜索更远的同源序列。

从蛋白质序列预测蛋白质3D结构的计算方法的发展沿着两条互补的路径前进,分别关注物理相互作用或进化历史。 物理相互作用方案将我们认知的分子驱动力(molecular driving forces)整合到物理热力学或动力学模拟或统计模型逼近中。 虽然理论上非常吸引人,但由于分子模拟的计算难度、蛋白质稳定性的上下文依赖性以及难以产生足够准确的蛋白质物理学模型,这种方法已被证明对即使是中等大小的蛋白质也极具挑战性。 近年来,进化方案提供了一种替代方案,其中蛋白质结构约束来自蛋白质进化历史的生物信息学分析、与已解决结构的同源性和成对进化相关性。

这种生物信息学方法极大地受益于蛋白质数据库 (PDB) 中存储的实验蛋白质结构的稳定增长、基因组测序的爆炸式增长以及用于解释这些相关性的深度学习技术的快速发展。 尽管取得了这些进展,基于物理和进化历史的方法产生的预测远低于实验准确性。

我们(DeepMind团队)开发了第一个能够在大多数情况下预测蛋白质结构接近实验准确性的计算方法。 我们开发的神经网络 AlphaFold 已进入 CASP14 评估(2020 年 5 月至 7 月)。 CASP 评估每两年进行一次,使用未在 PDB 中存放或公开披露的最近解决的结构,因此它是对参与方法的盲测,长期以来一直作为结构预测准确性的金标准评估。

蛋白质结构预测 (CASP) 实验的批判性评估旨在建立蛋白质结构预测的当前技术水平,确定已取得的进展,并突出未来可能最有成效的工作重点。

… 省略了部分内容

在图 2a 中证明,CASP14 中展示AlphaFold的高准确率扩展到最近大量PDB结构样本,其中所有结构在我们的训练数据截止后都存储在PDB中,并作为完整链进行分析。此外,当主链预测准确时,观察到高侧链准确性(图 2b),并且表明我们预测的局部距离差异测试(pLDDT)置信度可靠地预测了 Ca 局部距离差异测试(lDDT-Cα)准确度相应的预测(图 2c)。我们还发现可以准确估计全局叠加度量模板建模分数 (TM-score)(图 2d)。总体而言,这些验证了 AlphaFold 在 CASP14 蛋白上的高精度和可靠性也可以迁移到最近 PDB 提交的未经整理的数据集中。

Fig. 1 | AlphaFold produces highly accurate structures.

PDB蛋白质结构数据库(Protein Data Bank,简称PDB)是美国Brookhaven国家实验室于1971年创建的,由结构生物信息学研究合作组织(Research Collaboratory for Structural Bioinformatics,简称RCSB)维护。和核酸序列数据库一样,可以通过网络直接向PDB数据库提交数据。

AlphaFold网络

AlphaFold 通过结合基于蛋白质结构的进化、物理和几何约束的新型神经网络架构和训练程序,大大提高了结构预测的准确性。特别是,我们展示了一种联合嵌入多序列比对 (MSA) 和成对特征的新架构、一种新的输出表示和相关损失函数,可实现准确的端到端结构预测、新的等变注意架构、使用中间损失函数来实现预测的迭代改进,屏蔽 MSA 损失与结构联合训练,使用自蒸馏从未标记的蛋白质序列中学习,以及自我估计准确率。

蒸馏,就是知识蒸馏,将教师网络(teacher network)的知识迁移到学生网络(student network)上,使得学生网络的性能表现如教师网络一般;或者大型模型迁移到小型模型中,小型模型的参数规模小,运行速度快,但性能与大型模型参不多。

AlphaFold 网络使用一级氨基酸序列和同源物的比对序列作为输入,直接预测给定蛋白质的所有重原子的 3-D 坐标(图 1e,请参阅方法了解输入的详细信息,包括数据库、MSA 构建和使用的模板)。 最重要的方法和组件的描述如下在附件中提供了完整的网络架构和训练过程在附件的方法章节中。

该网络包括两个主要阶段。 首先,网络的主干通过我们称为 Evoformer 的新型神经网络块的重复层处理输入,以生成 Nseq × Nres 数组(Nseq:序列数,Nres:残基数),表示已处理的 MSA 和 Nres × Nres 数组,表示残基对。 MSA 表示是用原始 MSA 初始化的,但请参阅 附件-方法 1.2.7 了解处理非常深的 MSA 的详细信息。 Evoformer 块包含许多新颖的基于注意力和非基于注意力的组件。 我们在“可解释性”部分展示了证据,表明在 Evoformer 块中早期出现了具体的结构假设并不断完善。 Evoformer 模块的关键创新是在 MSA 内交换信息的新机制和允许直接推理空间和进化关系的配对表示。

网络的主干之后是结构模块,该模块以蛋白质的每个残基(全局刚体框架)的旋转和平移的形式引入了明确的 3-D 结构。这些表示在简单的状态下初始化,所有旋转设置为一致,所有位置设置为原点,但快速发展和完善,具有精确原子细节的高度准确的蛋白质结构。网络这一部分的关键创新包括打破链原子结构以允许同时对结构的所有部分进行局部细化,一种新颖的等变变换器允许网络隐式推理未表示的侧链原子,以及一个损失项代替残基的方向正确性的重要权重。在结构模块和整个网络中,我们通过反复将最终损失函数应用于输出,然后将输出递归地提供给相同的模块来强化迭代细化的概念。使用整个网络的迭代细化(我们称之为“循环”) 对准确性有显着贡献,而额外的训练时间很少(有关详细信息,请参阅附件-方法-1.8)。

Evoformer模块

名为 Evoformer(图 1e 和 3a)的网络构建块的关键原理是将蛋白质结构预测视为 3-D 空间中的图推理问题,其中图的边缘由邻近的残基定义。 配对表示的元素编码有关残基之间关系的信息(图 3b)。 MSA 表示的列编码输入序列的各个残基,而行表示这些残基出现的序列。 在这个框架内,我们定义了许多更新操作,这些更新操作应用于每个块中,其中不同的更新操作被串联应用。

Fig. 3 | Architectural details.

MSA 表示通过在 MSA 序列维度上求和的逐元素外积更新配对表示。 与之前的工作不同,此操作应用于每个块中,而不是在网络中应用一次,这使得从不断发展的 MSA 表示到配对表示的连续通信成为可能。

在配对表示中,有两种不同的更新模式。两者都受到配对表示一致性必要性的启发——为了将氨基酸的配对描述表示为单个 3-D 结构,必须满足许多约束,包括距离上的三角不等式。基于这种直觉,我们根据涉及三个不同节点的边三角形来安排对表示的更新操作(图 3c)。特别是,我们向轴向注意力添加了一个额外的 logit 偏差,以包括三角形的“缺失边”,并且我们定义了一个非注意力更新操作“三角形乘法更新”,它使用两条边来更新缺失的第三条边(参见 附件-方法-1.6.5 了解详情)。三角形乘法更新最初是作为注意力更对称且更便宜的替代品而开发的,仅使用注意力或乘法更新的网络都能够产生高精度结构。然而,两个更新的组合更准确。

我们还在 MSA 表示中使用了一种轴向注意力的变体。 在 MSA 中的 per-sequence attention 期间,我们从 pair stack 中投射额外的 logits 以偏置 MSA attention。 这通过提供从配对表示返回到 MSA 表示的信息流来关闭循环,确保整个 Evoformer 模块能够完全混合对和 MSA 表示之间的信息,并为结构模块中的结构生成做好准备。

端到端结构预测

结构模块(图 3d)使用对表示和来自主干的 MSA 表示的原始序列行(“single representation”,“单一表示”)在具体的 3-D 主干结构上运行。 3-D 主干结构表示为 Nres 独立的旋转和平移,每个旋转和平移相对于全局框架(residue gas,“残余气”?,图 3e)。这些旋转和平移,代表 N-Cα-C 原子的几何形状,优先考虑蛋白质骨架的方向,以便每个残基的侧链位置在该框架内受到高度限制。相反,肽键几何形状完全不受约束,并且在应用结构模块期间观察到网络经常违反链约束,因为打破此约束允许对链的所有部分进行局部细化,而无需解决复杂的闭环问题。在微调期间通过违规损失项鼓励满足肽键几何形状。只有在 Amber力场中的梯度下降结构的预测后松弛中才能实现肽键几何形状的精确执行。根据经验,这种最终松弛不会提高模型的准确性,如通过全局距离测试 (GDT) 或 IDDT-Cα34 测量的,但确实消除了分散注意力的立体化学违规而不会损失准确性。

AMBER力场是在生物大分子的模拟计算领域有着广泛应用的一个分子力场。AMBER力场的优势在于对生物大分子的计算,其对小分子体系的计算结果常常不能令人满意。

residue gas表示分两个阶段迭代更新(图 3d)。首先,我们称为不变点注意力 ( Point Attention) 的新型几何感知注意操作用于更新 Nres 神经激活集(single representation,“单一表示”)而不改变 3-D 位置,然后对residue gas使用更新的激活。不变点注意力通过在每个残基的局部框架中产生的 3-D 点来增强每个通常的注意力查询、键和值,这样最终值对全局旋转和平移是不变的(参见方法“不变点注意力(IPA)”了解详情)。 3-D 查询和键也对注意力施加了强烈的空间/局部性偏差,这非常适合蛋白质结构的迭代细化。在每个注意力操作和逐元素转换块之后,该模块计算每个主干帧的旋转和平移的更新。这些更新在每个残差的局部框架内的应用使得整体注意力和更新块成为对residue gas的等变操作。

侧链 chi 角的预测以及结构的最终每个残基精度 (pLDDT) 是在网络末端的最终激活上使用小的每个残基网络计算的。 TM 分数 (pTM) 的估计值是从成对错误预测中获得的,该预测被计算为最终对表示的线性投影。最后的损失(我们称之为帧对齐点误差(FAPE)(图 3f))将预测的原子位置与许多不同对齐下的真实位置进行比较。对于每个对齐,通过将预测帧 (Rk,tk) 对齐到相应的真实帧来定义,我们计算所有预测原子位置 xi 与真实原子位置的距离。由此产生的 Nframes × Natoms 距离受到限制的 L1 损失的惩罚。这对原子相对于每个残基的局部框架是正确的产生了强烈的偏见,因此在其侧链相互作用方面是正确的,并为 AlphaFold 提供了手性的主要来源(增刊-方法 1.9.3 和增刊-图 9)

使用标记和未标记数据进行训练

AlphaFold 架构能够仅使用对 PDB 数据的监督学习来训练到高精度,但我们能够使用类似于noisy-student自我蒸馏的方法来提高准确性(见图 4a)。 在这个过程中,我们使用一个训练有素的网络来预测来自 Uniclust30的约 350,000 个不同序列的结构,并将预测结构的新数据集过滤为高置信度子集。 然后,我们使用 PDB 和这个新的预测结构数据集的混合作为训练数据从头开始训练相同的架构,其中各种训练数据增强(例如裁剪和 MSA 子采样)使网络难以重现先前预测的结构。 这种自蒸馏过程有效地利用了未标记的序列数据,并显着提高了所得网络的准确性。

Self-training是最简单的半监督方法之一,其主要思想是找到一种方法,用未标记的数据集来扩充已标记的数据集。算法流程如下:
(1)首先,利用已标记的数据来训练一个好的模型,然后使用这个模型对未标记的数据进行标记。
(2)然后,进行伪标签的生成,因为我们知道,已训练好的模型对未标记数据的所有预测都不可能都是好的,因此对于经典的Self-training,通常是使用分数阈值过滤部分预测,以选择出未标记数据的预测标签的一个子集。
(3)其次,将生成的伪标签与原始的标记数据相结合,并在合并后数据上进行联合训练。
(4)整个过程可以重复n次,直到达到收敛。

此外,我们随机屏蔽或突变 MSA 中的单个残基,并具有来自 Transformers (BERT) 式目标的双向编码器表示来预测 MSA 序列的屏蔽元素。 这个目标鼓励网络学习解释系统发育和协变关系,而无需将特定的相关统计量硬编码到特征中。 与最近的独立工作相比,BERT 目标是在相同的训练示例上与正常的 PDB 结构损失联合训练的,并且没有进行预训练。

解释神经网络

为了了解 AlphaFold 如何预测蛋白质结构,我们为网络中的 48 个 Evoformer 块中的每一个训练了一个单独的结构模块,同时保持主网络的所有参数保持不变(补充方法 1.14)。包括我们的回收阶段,这提供了 192 个中间结构的轨迹,每个完整的 Evoformer 块一个,其中每个中间体代表网络对该块最可能结构的信念。在前几个块之后产生的轨迹出奇地平滑,表明 AlphaFold 对结构进行了不断的增量改进,直到它不能再改进为止(见图 4b 的准确度轨迹)。这些轨迹也说明了网络深度的作用。对于非常具有挑战性的蛋白质,如 SARS-CoV-2 Orf8 (T1064),网络搜索并重新排列多层的二级结构元素,然后再确定一个好的结构。对于 LmrP (T1024) 等其他蛋白质,网络会在前几层内找到最终结构。请参阅增刊-视频 1-4 的 CASP14 目标 T1024、T1044、T1064 和 T1091 的结构轨迹显示了一系列蛋白质大小和难度的清晰迭代构建过程。在补充-方法 1.16 和增刊-图12-13,我们解释了 AlphaFold 层产生的注意力图。

Fig. 4 | Interpreting the neural network

图 4a 包含 AlphaFold 组件的详细消融,表明各种不同的机制有助于 AlphaFold 的准确性。 请参阅增刊-方法 1.13 详细描述了每个消融模型、它们的训练细节、消融结果的扩展讨论以及 MSA 深度对每次消融的影响(补充-图 10)。

MSA 深度和跨链接触

虽然 AlphaFold 在绝大多数沉积的 PDB 结构中具有很高的准确性,但我们注意到仍然存在影响准确性或限制模型适用性的因素。当平均比对深度小于约 30 个序列时,该模型使用多个序列比对,准确度大幅下降(详见图 5a)。我们观察到阈值效应,其中 MSA 深度超过约 100 个序列的改进导致小增益。我们假设需要 MSA 信息在网络的早期阶段粗略地找到正确的结构,但是将该预测细化为高精度模型并不关键取决于 MSA 信息。我们观察到的另一个实质性限制是,与异型接触的数量相比,AlphaFold 对于链内或同型接触很少的蛋白质要弱得多。这通常发生在较大复合物中的桥接结构域,其中蛋白质的形状几乎完全由与复合物中其他链的相互作用产生。相反,AlphaFold 通常能够为同聚体提供高精度预测,即使链基本上交织在一起(例如图 5b)。我们希望 AlphaFold 的想法很容易适用于预测未来系统中的完整异质复合物,并且这将消除具有大量异质接触的蛋白质链的困难。

Fig. 5 | Effect of MSA depth and cross-chain contacts.

相关工作

蛋白质结构预测有一个漫长而多样的发展,在许多优秀的评论中得到了广泛的介绍。尽管将神经网络应用于结构预测的历史悠久,但它们最近才开始改进结构预测。这些方法通过将蛋白质结构预测问题处理为将进化耦合的“图像” 转换为蛋白质距离矩阵的“图像”,然后将距离预测集成到启发式系统中,从而有效地利用了计算机视觉系统的快速改进。产生最终的 3-D 坐标预测。最近开发了一些工作来直接预测 3-D 坐标 ,但这些方法的准确性与传统的手工结构预测管道不匹配 。同时,基于注意力的语言处理网络的成功 和最近的计算机视觉激发了对解释蛋白质序列的基于注意力的方法的探索。

讨论

我们在设计 AlphaFold 时采用的方法是生物信息学和物理方法的结合:我们使用物理和几何归纳偏差来构建组件,这些组件可以从 PDB 数据中学习,并最大限度地减少手工制作的特征(例如,AlphaFold 在没有氢的情况下有效地构建氢键债券评分函数)。这使得网络从 PDB 中的有限数据中更有效地学习,但能够应对结构数据的复杂性和多样性。特别是,AlphaFold 能够处理缺失的物理环境,并在具有挑战性的情况下生成准确的模型,例如交织的同源异构体或仅在未知血红素组存在时才会折叠的蛋白质。处理未指定结构条件的能力对于从 PDB 结构中学习至关重要,因为 PDB 代表了结构已被求解的所有条件。一般来说,AlphaFold 被训练产生最有可能作为 PDB 结构的一部分出现的蛋白质结构。在特定化学计量或配体/离子可单独从序列预测的情况下,AlphaFold 可能会产生一种隐式遵守这些约束的结构。

AlphaFold 已经向实验界展示了它的实用性,包括分子置换和解释低温电子显微镜 (cryo-EM) 图 。 此外,由于 AlphaFold 直接输出蛋白质坐标,因此 AlphaFold 会根据蛋白质序列的长度以图形处理单元 (GPU) 分钟到 GPU 小时生成预测(例如,对于 384 个残基,每个模型大约 1 个 GPU 分钟,请参阅方法了解详细信息 )。 这开辟了在蛋白质组规模及其他范围内预测结构的令人兴奋的可能性。

可用基因组测序技术和数据的爆炸式增长彻底改变了生物信息学,但实验结构测定的内在挑战阻止了我们结构知识的类似扩展。 通过开发准确的蛋白质结构预测算法,再加上由实验社区组装的现有大型且精心准备的结构和序列数据库,我们希望加速结构生物信息学的进步,以跟上基因组学革命的步伐。 我们希望 AlphaFold 以及将其技术应用于其他生物物理问题的计算方法将成为现代生物学的重要工具。

在线内容

任何方法、附加参考资料、自然研究报告摘要、源数据、扩展数据、补充信息、致谢、同行评审信息; 作者贡献和竞争利益的详细信息; 数据和代码可用性声明可在 Highly accurate protein structure prediction with AlphaFold – Nature 上获得。

 Highly accurate protein structure prediction with AlphaFold – Nature

dijkstra算法求最短路径

基本思想

dijkstra的算法思想是从以上最短距离数组中每次选择一个最近的点,将其作为下一个点,然后重新计算从起始点经过该点到其他所有点的距离,更新最短距离数据。已经选取过的点就是确定了最短路径的点,不再参与下一次计算。

  1. 通过Dijkstra计算图G中的最短路径时,需要指定起点s(即从顶点s开始计算)。
  2. 此外,引进两个集合S和U。S的作用是记录已求出最短路径的顶点(以及相应的最短路径长度),而U则是记录还未求出最短路径的顶点(以及该顶点到起点s的距离)。
  3. 初始时,S中只有起点s;U中是除s之外的顶点,并且U中顶点的路径是”起点s到该顶点的路径”。然后,从U中找出路径最短的顶点,并将其加入到S中;接着,更新U中的顶点和顶点对应的路径。 然后,再从U中找出路径最短的顶点,并将其加入到S中;接着,更新U中的顶点和顶点对应的路径。 … 重复该操作,直到遍历完所有顶点。

操作步骤

  1. 初始时,S只包含起点s;U包含除s外的其他顶点,且U中顶点的距离为”起点s到该顶点的距离”[例如,U中顶点v的距离为(s,v)的长度,然后s和v不相邻,则v的距离为∞]。
  2. 从U中选出”距离最短的顶点k”,并将顶点k加入到S中;同时,从U中移除顶点k。
  3. 更新U中各个顶点到起点s的距离。之所以更新U中顶点的距离,是由于上一步中确定了k是求出最短路径的顶点,从而可以利用k来更新其它顶点的距离;例如,(s,v)的距离可能大于(s,k)+(k,v)的距离。
  4. 重复步骤(2)和(3),直到遍历完所有顶点。

第1步:初始化距离,其实指与D直接连接的点的距离。dis[c]代表D到C点的最短距离,因而初始dis[C]=3,dis[E]=4,dis[D]=0,其余为无穷大。设置集合S用来表示已经找到的最短路径。此时,S={D}。现在得到D到各点距离{D(0),C(3),E(4),F(*),G(*),B(*),A(*)},其中*代表未知数也可以说是无穷大,括号里面的数值代表D点到该点的最短距离。

第2步:不考虑集合S中的值,因为dis[C]=3,是当中距离最短的,所以此时更新S,S={D,C}。接着我们看与C连接的点,分别有B,E,F,已经在集合S中的不看,dis[C-B]=10,因而dis[B]=dis[C]+10=13,dis[F]=dis[C]+dis[C-F]=9,dis[E]=dis[C]+dis[C-E]=3+5=8>4(初始化时的dis[E]=4)不更新。此时{D(0),C(3),E(4),F(9),G(*),B(13),A(*)}。

第3步:在第2步中,E点的值4最小,更新S={D,C,E},此时看与E点直接连接的点,分别有F,G。dis[F]=dis[E]+dis[E-F]=4+2=6(比原来的值小,得到更新),dis[G]=dis[E]+dis[E-G]=4+8=12(更新)。此时{D(0),C(3),E(4),F(6),G(12),B(13),A(*)}

第4步:在第3步中,F点的值6最小,更新S={D,C,E,F},此时看与F点直接连接的点,分别有B,A,G。dis[B]=dis[F]+dis[F-B]=6+7=13,dis[A]=dis[F]+dis[F-A]=6+16=22,dis[G]=dis[F]+dis[F-G]=6+9=15>12(不更新)。此时{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

第5步:在第4步中,G点的值12最小,更新S={D,C,E,F,G},此时看与G点直接连接的点,只有A。dis[A]=dis[G]+dis[G-A]=12+14=26>22(不更新)。{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

第6步:在第5步中,B点的值13最小,更新S={D,C,E,F,G,B},此时看与B点直接连接的点,只有A。dis[A]=dis[B]+dis[B-A]=13+12=25>22(不更新)。{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

第6步:最后只剩下A值,直接进入集合S={D,C,E,F,G,B,A},此时所有的点都已经遍历结束,得到最终结果{D(0),C(3),E(4),F(6),G(12),B(13),A(22)}.

python实现:

将以上的过程使用python来实现。
首先总结一个Dijkstra算法的核心思想,分成两步走:

  1. 构造一个最短路径数组,每次找到数组中未访问的节点里最小的点
  2. 以上一步的节点为最新节点,更新起始点到所有点的距离

使用python就是实现这两步即可



MAX= float('inf')

matrix = [
    [0,10,MAX,4,MAX,MAX],
    [10,0,8,2,6,MAX],
    [MAX,8,10,15,1,5],
    [4,2,15,0,6,MAX],
    [MAX,6,1,6,0,12],
    [MAX,MAX,5,MAX,12,0]
    ]


def dijkstra(matrix, start_node):
    
    #矩阵一维数组的长度,即节点的个数
    matrix_length = len(matrix)

    #访问过的节点数组
    used_node = [False] * matrix_length

    #最短路径距离数组
    distance = [MAX] * matrix_length

    #初始化,将起始节点的最短路径修改成0
    distance[start_node] = 0
    
    #将访问节点中未访问的个数作为循环值,其实也可以用个点长度代替。
    while used_node.count(False):
        min_value = float('inf')
        min_value_index = 999
        
        #在最短路径节点中找到最小值,已经访问过的不在参与循环。
        #得到最小值下标,每循环一次肯定有一个最小值
        for index in range(matrix_length):
            if not used_node[index] and distance[index] < min_value:
                min_value = distance[index]
                min_value_index = index
        
        #将访问节点数组对应的值修改成True,标志其已经访问过了
        used_node[min_value_index] = True

        #更新distance数组。
        #以B点为例:distance[x] 起始点达到B点的距离,
        #distance[min_value_index] + matrix[min_value_index][index] 是起始点经过某点达到B点的距离,比较两个值,取较小的那个。
        for index in range(matrix_length):
            distance[index] = min(distance[index], distance[min_value_index] + matrix[min_value_index][index])

    return distance




start_node = int(input('请输入起始节点:'))
result = dijkstra(matrix,start_node)
print('起始节点到其他点距离:%s' % result)

数据结构 — 单调队列

顾名思义,单调队列的重点分为 “单调” 和 “队列”

“单调” 指的是元素的的 “规律”——递增(或递减)

“队列” 指的是元素只能从队头和队尾进行操作

  • 单调递增队列:保证队列头元素一定是当前队列的最小值,用于维护区间的最小值。
  • 单调递减队列:保证队列头元素一定是当前队列的最大值,用于维护区间的最大值。

1、对于单调递增队列,对于一个元素如果它大于等于队尾元素, 那么直接把元素进队, 如果这个元素小于队尾元素,则将队尾元素出队列,直到满足这个元素大于等于队尾元素。

2、对于单调递减队列,对于一个元素如果它小于等于队尾元素, 那么直接把元素进队, 如果这个元素大于队尾元素,则将队尾元素出队列,直到满足这个元素小于等于队尾元素。

3、判断队头元素是否在待求解的区间之内,如果不在,就将其删除。(这个很好理解呀,因为单调队列的队头元素就是待求解区间的极值)

4、经过上面两个操作,取出 队列的头元素 ,就是 当前区间的极值

可以发现队列中的元素都是单调递减的(不然也就不叫单调递减队列啦),同时既有入队列的操作、也有出队列的操作。

实现单调队列,主要分为三个部分:

  • 去尾操作队尾元素出队列。当队列有新元素待入队,需要从队尾开始,删除影响队列单调性的元素,维护队列的单调性。(删除一个队尾元素后,就重新判断新的队尾元素)

去尾操作结束后,将该新元素入队列。

  • 删头操作队头元素出队列。判断队头元素是否在待求解的区间之内,如果不在,就将其删除。(这个很好理解呀,因为单调队列的队头元素就是待求解区间的极值)
  • 取解操作 :经过上面两个操作,取出 队列的头元素 ,就是 当前区间的极值

参考代码

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#define maxn 1000100
using namespace std;
int q[maxn], a[maxn];
int n, k;

void getmin() {  // 得到这个队列里的最小值,直接找到最后的就行了
  int head = 0, tail = 0;
  for (int i = 1; i < k; i++) {
    while (head <= tail && a[q[tail]] >= a[i]) tail--;
    q[++tail] = i;
  }
  for (int i = k; i <= n; i++) {
    while (head <= tail && a[q[tail]] >= a[i]) tail--;
    q[++tail] = i;
    while (q[head] <= i - k) head++;
    printf("%d ", a[q[head]]);
  }
}

void getmax() {  // 和上面同理
  int head = 0, tail = 0;
  for (int i = 1; i < k; i++) {
    while (head <= tail && a[q[tail]] <= a[i]) tail--;
    q[++tail] = i;
  }
  for (int i = k; i <= n; i++) {
    while (head <= tail && a[q[tail]] <= a[i]) tail--;
    q[++tail] = i;
    while (q[head] <= i - k) head++;
    printf("%d ", a[q[head]]);
  }
}

int main() {
  scanf("%d%d", &n, &k);
  for (int i = 1; i <= n; i++) scanf("%d", &a[i]);
  getmin();
  printf("\n");
  getmax();
  printf("\n");
  return 0;
}

torch.nn.utils.rnn中的pack_padded_sequence和 pad_packed_sequence and pad_sequence

这两个是互逆的操作。

先来看为什么需要pad和pack操作:

先看一个例子,这个batch中有5个sample

如果不用pack和pad操作会有一个问题,什么问题呢?比如上图,句子“Yes”只有一个单词,但是padding了多余的pad符号,这样会导致LSTM对它的表示通过了非常多无用的字符,这样得到的句子表示就会有误差,更直观的如下图:

那么我们正确的做法应该是怎么样呢?

在上面这个例子,我们想要得到的表示仅仅是LSTM过完单词”Yes”之后的表示,而不是通过了多个无用的“Pad”得到的表示:如下图:

torch.nn.utils.rnn.pack_padded_sequence()

这里的pack,理解成压紧比较好。 将一个 填充过的变长序列 压紧。(填充时候,会有冗余,所以压紧一下)

其中pack的过程为:(注意pack的形式,不是按行压,而是按列压)

pack之后,原来填充的 PAD(一般初始化为0)占位符被删掉了。

输入的形状可以是(T×B×* )。T是最长序列长度,Bbatch size*代表任意维度(可以是0)。如果batch_first=True的话,那么相应的 input size 就是 (B×T×*)

Variable中保存的序列,应该按序列长度的长短排序,长的在前,短的在后。即input[:,0]代表的是最长的序列,input[:, B-1]保存的是最短的序列。

NOTE: 只要是维度大于等于2的input都可以作为这个函数的参数。你可以用它来打包labels,然后用RNN的输出和打包后的labels来计算loss。通过PackedSequence对象的.data属性可以获取 Variable

参数说明:

  • input (Variable) – 变长序列 被填充后的 batch
  • lengths (list[int]) – Variable 中 每个序列的长度。
  • batch_first (bool, optional) – 如果是True,input的形状应该是B*T*size
  • 参数 enforce_sorted ,如果是 True ,则输入应该是按长度降序排序的序列。如果是 False ,会在函数内部进行排序。默认值为 True 。也就是说在输入 pack_padded_sequence 前,我们也可以不对数据进行排序。

返回值:

一个PackedSequence 对象。

torch.nn.utils.rnn.pad_packed_sequence()

填充packed_sequence

上面提到的函数的功能是将一个填充后的变长序列压紧。 这个操作和pack_padded_sequence()是相反的。把压紧的序列再填充回来。填充时会初始化为0。

返回的Varaible的值的size是 T×B×*T 是最长序列的长度,B 是 batch_size,如果 batch_first=True,那么返回值是B×T×*

Batch中的元素将会以它们长度的逆序排列。

参数说明:

  • sequence (PackedSequence) – 将要被填充的 batch
  • batch_first (bool, optional) – 如果为True,返回的数据的格式为 B×T×*

返回值: 一个tuple,包含被填充后的序列,和batch中序列的长度列表

实例 代码:



>>> from torch.nn.utils.rnn import pack_padded_sequence, pad_packed_sequence
>>> seq = torch.tensor([[1,2,0], [3,0,0], [4,5,6]])
#seq的维度是3*3,第一个 维度表示T,is the length of the longest sequence 
#第二维表示B,批次大小,也就是说有3个长度为3的向量,其中每列表示一个序列,序列都是等长的,因为短的序列已经用0补齐了
>>> lens = [2, 1, 3]
>>> packed = pack_padded_sequence(seq, lens, batch_first=True, enforce_sorted=False)
#将补齐的序列压紧成一个序列,将0去掉
>>> packed
PackedSequence(data=tensor([4, 1, 3, 5, 2, 6]), batch_sizes=tensor([3, 2, 1]),
               sorted_indices=tensor([2, 0, 1]), unsorted_indices=tensor([1, 2, 0]))
>>> seq_unpacked, lens_unpacked = pad_packed_sequence(packed, batch_first=True)
>>> seq_unpacked
tensor([[1, 2, 0],
        [3, 0, 0],
        [4, 5, 6]])
>>> lens_unpacked
tensor([2, 1, 3])

pad_sequence

参数

sequences:表示输入样本序列,为 list 类型,list 中的元素为 tensor 类型。 tensor 的 size 为 L * F 。其中,L 为单个序列的长度,F 为序列中每个时间步(time step)特征的个数,根据任务的不同 F 的维度会有所不同。

batch_first:为 True 对应 [batch_size, seq_len, feature];False 对应[seq_len, batch_size, feature],从习惯上来讲一般设置为 True 比较符合我们的认知。

padding_value:填充值,默认值为 0 。

说明

主要用来对样本进行填充,填充值一般为 0 。我们在训练网络时,一般会采用一个一个 mini-batch 的方式,将训练样本数据喂给网络。在 PyTorch 里面数据都是以 tensor 的形式存在,一个 mini-batch 实际上就是一个高维的 tensor ,每个序列数据的长度必须相同才能组成一个 tensor 。为了使网络可以处理 mini-batch 形式的数据,就必须对序列样本进行填充,保证一个 mini-batch 里面的数据长度是相同的。

在 PyTorch 里面一般是使用 DataLoader 进行数据加载,返回 mini-batch 形式的数据,再将此数据喂给网络进行训练。我们一般会自定义一个 collate_fn 函数,完成对数据的填充。

import torch
from torch.utils.data import Dataset, DataLoader
from torch.nn.utils.rnn import pad_sequence,pack_padded_sequence,pack_sequence,pad_packed_sequence

class MyData(Dataset):
    def __init__(self, data):
        self.data = data

    def __len__(self):
        return len(self.data)

    def __getitem__(self, idx):
        return self.data[idx]

def collate_fn(data):
    data.sort(key=lambda x: len(x), reverse=True)
    data = pad_sequence(data, batch_first=True, padding_value=0)
    return data

a = torch.tensor([1,2,3,4])
b = torch.tensor([5,6,7])
c = torch.tensor([7,8])
d = torch.tensor([9])
train_x = [a, b, c, d]

data = MyData(train_x)
data_loader = DataLoader(data, batch_size=2, shuffle=True, collate_fn=collate_fn)
# 采用默认的 collate_fn 会报错
#data_loader = DataLoader(data, batch_size=2, shuffle=True) 
batch_x = iter(data_loader).next()

运行程序,得到 batch_x 的值:

# batch_x
tensor([[1, 2, 3, 4],
        [9, 0, 0, 0]])

从 batch_x 的值可以看出,第二行填充了三个 0 ,使其长度和第一行保持一致。

需要说明的是,对于长度不同的序列,使用默认的 collate_fn 函数,不自定义 collate_fn 函数完成对序列的填充,上面的程序就会报错。

Python3 迭代器与生成器

迭代器

迭代是Python最强大的功能之一,是访问集合元素的一种方式。

迭代器是一个可以记住遍历的位置的对象。

迭代器对象从集合的第一个元素开始访问,直到所有的元素被访问完结束。迭代器只能往前不会后退。

迭代器有两个基本的方法:iter() 用于生成迭代器 和 next() 求迭代器的值

字符串,列表或元组对象都可用于创建迭代器:

>>> list=[1,2,3,4]
>>> it = iter(list)    # 创建迭代器对象
>>> print (next(it))   # 输出迭代器的下一个元素
1
>>> print (next(it))
2
>>>

迭代器对象可以使用常规for语句进行遍历:

list=[1,2,3,4]
it = iter(list)    # 创建迭代器对象
for x in it:
    print (x, end=" ")

执行以上程序,输出结果如下:

1 2 3 4

也可以使用 next() 函数:

import sys         # 引入 sys 模块
 
list=[1,2,3,4]
it = iter(list)    # 创建迭代器对象
 
while True:
    try:
        print (next(it))
    except StopIteration:
        sys.exit()

执行以上程序,输出结果如下:

1
2
3
4

创建一个迭代器

把一个类作为一个迭代器使用需要在类中实现两个方法 __iter__() 与 __next__() 。

如果你已经了解的面向对象编程,就知道类都有一个构造函数,Python 的构造函数为 __init__(), 它会在对象初始化的时候执行。

更多内容查阅:Python3 面向对象

__iter__() 方法返回一个特殊的迭代器对象, 这个迭代器对象实现了 __next__() 方法并通过 StopIteration 异常标识迭代的完成。

__next__() 方法(Python 2 里是 next())会返回下一个迭代器对象。

创建一个返回数字的迭代器,初始值为 1,逐步递增 1:

class MyNumbers:
  def __iter__(self):
    self.a = 1
    return self
 
  def __next__(self):
    x = self.a
    self.a += 1
    return x
 
myclass = MyNumbers()
myiter = iter(myclass)
 
print(next(myiter))
print(next(myiter))
print(next(myiter))
print(next(myiter))
print(next(myiter))

执行输出结果为:

1
2
3
4
5

StopIteration

StopIteration 异常用于标识迭代的完成,防止出现无限循环的情况,在 __next__() 方法中我们可以设置在完成指定循环次数后触发 StopIteration 异常来结束迭代。

在 20 次迭代后停止执行:

class MyNumbers:
  def __iter__(self):
    self.a = 1
    return self
 
  def __next__(self):
    if self.a <= 20:
      x = self.a
      self.a += 1
      return x
    else:
      raise StopIteration
 
myclass = MyNumbers()
myiter = iter(myclass)
 
for x in myiter:
  print(x)

执行输出结果为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

生成器

yield是python的一个关键字,本质上是一个生成器generator。

生成器是一种特殊的函数,它会返回一个迭代器。定义一个生成器函数同定义一个普通函数没有什么区别,特殊之处在于生成器函数内部会包含yield表达式,专门用于生成一个序列。当一个生成器函数被调用时,它会返回一个迭代器。之后就由这个迭代器来控制生成器函数的执行。当生成器函数被调用后,首先会执行到第一个yield表达式处,然后会将生成器函数挂起,将yield生成的表达式的值返回给生成器函数的调用者。当生成器函数被挂起时,它的所有局部状态都会被保存起来,包括当前绑定的局部变量、指令指针、内部栈和异常处理的状态。当调用迭代器的方法时,一般都是调用next()方法,将会恢复生成器函数的执行,并且是从上次被挂起的地方继续执行,直到遇到另外一次yield调用,生成器函数将再次被挂起。

在一个生成器函数中,如果没有 return,则默认执行至函数完毕,如果在执行过程中 return,则直接抛出 StopIteration 终止迭代。

可以中断一个函数的执行,跳转到调用的地方

在 Python 中,使用了 yield 的函数被称为生成器(generator)。

跟普通函数不同的是,生成器是一个返回迭代器的函数,只能用于迭代操作,更简单点理解生成器就是一个迭代器。

在调用生成器运行的过程中,每次遇到 yield 时函数会暂停并保存当前所有的运行信息,返回 yield 的值, 并在下一次执行 next() 方法时从当前位置继续运行。

调用一个生成器函数,返回的是一个迭代器对象。

以下实例使用 yield 实现斐波那契数列:

import sys
 
def fibonacci(n): # 生成器函数 - 斐波那契
    a, b, counter = 0, 1, 0
    while True:
        if (counter > n): 
            return
        yield a
        a, b = b, a + b
        counter += 1
f = fibonacci(10) # f 是一个迭代器,由生成器返回生成
 
while True:
    try:
        print (next(f), end=" ")
    except StopIteration:
        sys.exit()

执行以上程序,输出结果如下:

0 1 1 2 3 5 8 13 21 34 55

数据结构 — 堆

python中直接使用heapq建立 小根堆

http://139.9.1.231/index.php/2022/02/22/python3heapq/

堆简介

堆是一棵树,其每个节点都有一个键值,且每个节点的键值都大于等于/小于等于其父亲的键值。

每个节点的键值都大于等于其父亲键值的堆叫做小根堆,否则叫做大根堆。STL 中的 priority_queue 其实就是一个大根堆。

(小根)堆主要支持的操作有:插入一个数、查询最小值、删除最小值、合并两个堆、减小一个元素的值。

一些功能强大的堆(可并堆)还能(高效地)支持 merge 等操作。

一些功能更强大的堆还支持可持久化,也就是对任意历史版本进行查询或者操作,产生新的版本。

二叉堆

结构

从二叉堆的结构说起,它是一棵二叉树,并且是完全二叉树,每个结点中存有一个元素(或者说,有个权值)。

堆性质:父亲的权值不小于儿子的权值(大根堆)。同样的,我们可以定义小根堆。本文以大根堆为例。

由堆性质,树根存的是最大值(getmax 操作就解决了)。

插入操作

插入操作是指向二叉堆中插入一个元素,要保证插入后也是一棵完全二叉树。

最简单的方法就是,最下一层最右边的叶子之后插入。

如果最下一层已满,就新增一层。

插入之后可能会不满足堆性质?

向上调整:如果这个结点的权值大于它父亲的权值,就交换,重复此过程直到不满足或者到根。

可以证明,插入之后向上调整后,没有其他结点会不满足堆性质。

向上调整的时间复杂度是  的。

二叉堆的插入操作

删除操作

删除操作指删除堆中最大的元素,即删除根结点。

但是如果直接删除,则变成了两个堆,难以处理。

所以不妨考虑插入操作的逆过程,设法将根结点移到最后一个结点,然后直接删掉。

然而实际上不好做,我们通常采用的方法是,把根结点和最后一个结点直接交换。

于是直接删掉(在最后一个结点处的)根结点,但是新的根结点可能不满足堆性质……

向下调整:在该结点的儿子中,找一个最大的,与该结点交换,重复此过程直到底层。

可以证明,删除并向下调整后,没有其他结点不满足堆性质。

减小某个点的权值

很显然,直接修改后,向上调整一次即可

参考代码:


void up(int x) {
  while (x > 1 && h[x] > h[x / 2]) {
    swap(h[x], h[x / 2]);
    x /= 2;
  }
}

void down(int x) {
  while (x * 2 <= n) {
    t = x * 2;
    if (t + 1 <= n && h[t + 1] > h[t]) t++;
    if (h[t] <= h[x]) break;
    std::swap(h[x], h[t]);
    x = t;
  }
}

建堆

方法一:使用 decreasekey(即,向上调整)¶
从根开始,按 BFS 序进行。


void build_heap_1() {
  for (i = 1; i <= n; i++) up(i);
}
方法二:使用向下调整¶
这时换一种思路,从叶子开始,逐个向下调整


void build_heap_2() {
  for (i = n; i >= 1; i--) down(i);
}

并查集

一、概念及其介绍

并查集是一种树型的数据结构,用于处理一些不相交集合的合并及查询问题。

并查集的思想是用一个数组表示了整片森林(parent),树的根节点唯一标识了一个集合,我们只要找到了某个元素的的树根,就能确定它在哪个集合里。

二、适用说明

并查集用在一些有 N 个元素的集合应用问题中,我们通常是在开始时让每个元素构成一个单元素的集合,然后按一定顺序将属于同一组的元素所在的集合合并,其间要反复查找一个元素在哪个集合中。这个过程看似并不复杂,但数据量极大,若用其他的数据结构来描述的话,往往在空间上过大,计算机无法承受,也无法在短时间内计算出结果,所以只能用并查集来处理。

三、并查集的基本数据表示

如上图 0-4 下面都是 05-9 下面都是 1,表示 0、1、2、3、4 这五个元素是相连接的,5、6、7、8、9 这五个元素是相连的。

再如上图 0、2、4、6、8 下面都是 0 这个集合,表示 0、2、4、6、8 这五个元素是相连接的,1、3、5、7、9 下面都是 1 这个集合,表示 0,1、3、5、7、9 这五个元素是相连的。

构造一个类 UnionFind,初始化, 每一个id[i]指向自己, 没有合并的元素:

public UnionFind1(int n) {
        count = n;
        id = new int[n];
        // 初始化, 每一个id[i]指向自己, 没有合并的元素
        for (int i = 0; i < n; i++)
            id[i] = i;
    }

查找

通俗地讲一个故事:几个家族进行宴会,但是家族普遍长寿,所以人数众多。由于长时间的分离以及年龄的增长,这些人逐渐忘掉了自己的亲人,只记得自己的爸爸是谁了,而最长者(称为「祖先」)的父亲已经去世,他只知道自己是祖先。为了确定自己是哪个家族,他们想出了一个办法,只要问自己的爸爸是不是祖先,一层一层的向上问,直到问到祖先。如果要判断两人是否在同一家族,只要看两人的祖先是不是同一人就可以了。

在这样的思想下,并查集的查找算法诞生了。

# Python Version
fa = [0] * MAXN # 记录某个人的爸爸是谁,特别规定,祖先的爸爸是他自己

# 递归
def find(x):
    # 寻找x的祖先
    if fa[x] == x:
        return x # 如果x是祖先则返回
    else:
        return find(fa[x]) # 如果不是则 x 的爸爸问 x 的爷爷

# 非递归
def find(x):
    while x != fa[x]: # 如果 x 不是祖先,就一直往上一辈找
        x = fa[x]
    return x # 如果 x 是祖先则返回

这样的确可以达成目的,但是显然效率实在太低。为什么呢?因为我们使用了太多没用的信息,我的祖先是谁与我父亲是谁没什么关系,这样一层一层找太浪费时间,不如我直接当祖先的儿子,问一次就可以出结果了。甚至祖先是谁都无所谓,只要这个人可以代表我们家族就能得到想要的效果。把在路径上的每个节点都直接连接到根上,这就是路径压缩。

# Python Version
def find(x):
    if x != fa[x]: # x 不是自身的父亲,即 x 不是该集合的代表
        fa[x] = find(fa[x]) # 查找 x 的祖先直到找到代表,于是顺手路径压缩
    return fa[x]

合并

宴会上,一个家族的祖先突然对另一个家族说:我们两个家族交情这么好,不如合成一家好了。另一个家族也欣然接受了。
我们之前说过,并不在意祖先究竟是谁,所以只要其中一个祖先变成另一个祖先的儿子就可以了。

# Python Version
def unionSet(x, y):
    # x 与 y 所在家族合并
    x = find(x)
    y = find(y)
    fa[x] = y # 把 x 的祖先变成 y 的祖先的儿子

AI制药 : 分子对接

任务:

1、安装学习yutobe上面的软件

https://www.youtube.com/watch?v=Sux91FJ3Xe8&t=629s

2、跑通论文代码

https://www.youtube.com/watch?v=Sux91FJ3Xe8&t=629s

Pymol简介

Pymol是一款操作简单,功能强大的分子以及蛋白的可视化软件,由薛定谔公司研发,科研人员可以从官网申请最新教育版本,同时pymo的开源版(https://github.com/schrodinger/pymol-open-source),可以直接从网站上下载,但是版本较老。所以,根据需求选择版本进行下载。 说明: https://cloud.tencent.com/developer/article/1785088

Pymol入门教程:

http://pymol.chenzhaoqiang.com/intro/startManual.html

分子对接教程

1、

https://cloud.tencent.com/developer/inventory/15332

2、https://www.bilibili.com/video/av466685164?from=search&seid=9870338638011316620&spm_id_from=333.337.0.0

vina只负责对接,mgltool负责提供蛋白质分子和配体分子。

先用MGL 生成vina需要的pdbqt文件

MGL tools的作用就是生成pdbqt文件

1、:打开MGL tools,打开受体蛋白的pdb :file-》read molecule

2、蛋白质的pdb(数据库):

关于蛋白质结构的PDB文件,做分子对接,估计大家都知道PDB这个蛋白质数据库啦。这里简单的介绍一下。

蛋白质的三级结构是指整条多肽链的三维空间结构,也就是包括碳骨架和侧链在内的所有原子的空间排列。第一个蛋白质的三维空间结构于 1958 年用 X-射线衍射法(X-ray Crystallography)测定。这种方法目前仍然是获取蛋白质三级结构的主要方法。PDB 数据库中绝大多数蛋白质结构都是用这种方法测定的。另一个测定蛋白质三维空间结构的方法是核磁共振法(Nuclear Magnetic Resonance, NMR)。无法结晶的蛋白质,可以利用核磁共振法在液体环境中进行结构测定。但是核磁共振法只能用于质量小于 70 千道尔顿的分子,大约对应 200 个氨基酸的长度。除此之外,还有一些不太常用的方法也可以测定分子的三维空间结构,比如冷冻电子显微镜技术(Cyro-Electron Microscopy)。无论用什么方法测定的空间结构,都要提交到 PDB 数据库。所以我们获取蛋白质三级结构最直接的办法就是去PDB 搜索(http://www.rcsb.org/)。 从PDB首页的搜索条里,可以通过搜索PDB ID、分子名称、作者姓名等关键词来查找蛋白质三级结构。此外,利用高级搜索工具,可以通过序列相似性搜索获得与输入序列在序列水平上相似的蛋白质的三级结构。搜索方法选 BLAST,输入序列,点击“Result Count”。这里不详细介绍,因为我们做分子对接,通常蛋白名称是已知的。我们重点介绍怎么选择合适的蛋白结构文件。 比如我们搜索PI3K这个蛋白,直接在搜索栏搜索,结果是有很多的。可以看到有393个结构信息。首先我们可以通过左边的栏进行筛选,比如物种信息,我们选择人。当然,结果的显示排序可通过结果上面的选项卡进行选择不同的排序方式。我们筛选合适的蛋白结构,常用Score这个选项.我们选择分辨率较好的在前。这里的0.9Å,Å是光波长度和分子直径的常用计量单位,值越小,分辨率越高,结构越准确。页面往下拉,可以看见这个值越来越大,我们优先选择值小的。我们可以从页面里面看见一下基本信息,比如方法,物种以及被解析的时间等。这里5GJI这个结构获取的方法就是X-RAY。我们点击这个蛋白,进入后可以看见详细的信息。然后我们还要看这个蛋白的描述是不是我们想要的蛋白,从这里面感觉看起来比较费劲。这里我们借助uniprot这个数据库来选择是比较方便的。这里简单介绍一下这个数据库,可能有的同学是第一次知道。翻了多年前的笔记,粘贴在下面。 UniProt 数据库有三个层次。

第一层叫 UniParc,收录了所有 UniProt 数据库子库中的蛋白质序列,量大,粗糙。

第二层是 UniRef,他归纳了 UniProt 几个主要数据库并且是将重复序列去除后的数据库。

第三层是 UniProtKB,他有详细注释并与其他数据库有链接,分为 UniProtKB 下的 Swiss-Prot和 UniProtKB 下的 TrEMBL 数据库。

关系稍有点复杂,但实际上我们最常用的就是 UniProtKB下的 Swiss-Prot 数据库。

从 UniProt 数据库查看一条蛋白质序列(http://www.uniprot.org/)。在UniProt数据库的首页上有一个关于 UniProtKB 数据库的统计表。可以看到,TrEMBL 数据库里存储的序列数量远远大于 Swiss-Prot 中的。统计表里清楚的写着:TrEMBL 是自动注释的,没有经过检查,而 Swiss-Prot 是人工注释的,并且经过检查。

然后点击下载文件就可以直接下载PDB格式的蛋白结构文件。下载的PDB文件可以用pymol或者VMD观察结构。能够实现蛋白质三维结构可视化的软件非常多。比专业级的PyMOL(https://pymol.org/2/)。这个软件已经被世界上著名的生物医药软件公司“薛定谔公司(Schrödinger)”收购。这种专业级的可视化软件不仅能够做出非常漂亮的图片,它还有强大的插件支持各种各样的蛋白质结构分析,这款软件需要购买,如果你发表的文章里提到某些内容是使用PyMOL制作的,而文章中所有作者和作者单位都没有PyMOL的购买记录的话,你可能会面临薛定谔公司的追责。 如果要对接的蛋白没有结构,我们又要对接,那就只能是自己通过软件预测了。蛋白质结构预测的方法有从头计算法,同源建模法,穿线法和综合法。常用的是同源建模法,SWISS-MODEL(www.swissmodel.expasy.org)就是一款用同源建模法预测蛋白质三级结构的全自动软件,这里不详细介绍了,预测的模型还要涉及模型好坏的评价,后续有时间,再介绍蛋白质三级结构的预测。

接下来我们打开AutoDockTools(ADT),打开我们前面保存的文件1E8Y_PYMOL.pdb

删除水分子和其他配体,常规操作不用解释

然后计算电荷和添加原子类型

Edit–Charges–Compute Gasteiger

Edit–Atoms–Assign AD4 type

就可以导出成pbdqt格式的文件了

然后右键吧蛋白删除掉,导入配体小分子,随便从ZINC下了一个

ZINC(http://zinc.docking.org/) 还有一个数据库能下载mol2格式的文件。ZINC),这里就不介绍了,你要是能从上面的数据库下载到你配体小分子的mol2格式文件,就直接用,如果不能,那就是去PubChem数据库(https://pubchem.ncbi.nlm.nih.gov/)下载sdf文件,然后进行转换,这也是我这里要介绍的。

Ligand-input-open

子对接教程

1 分子对接的工作环境
1.1 对接基本需求
受体的 pdb 文件、知道配体的结构
1.2 对接软件
使用 AutoDock 进行对接,该软件由 AutoDock 和 MGLTools 两
部分组成,AutoDock 为主程序,仅提供了命令行接口,而 MGLTools
中的 AutoDockTools(ADT)可以看作是 AutoDock 的图形用户界面。
配体使用 ChemDraw 和 Chem3D 绘制,该软件为收费软件,可
使用其它免费的结构式绘制工具进行结构绘制,使用 Open Babel 转
3D 结构和文件格式的转换。
使用 PyMol 开源版本对受体蛋白的 pdb 文件进行处理。

2 准备受体、配体的 pdbqt 文件


AutoDock 只接受 pdbqt 格 式 的 文 件 , 所 以 需 要 通 过
AutoDockTools 将 pdb 或者其他格式的文件转化为 pdbqt 文件。


2.1 准备受体的 pdbqt 文件
2.1.1 使用 PyMol 软件进行处理


pdb 下载的蛋白需要先使用 PyMol 软件删除多余的离子、水分
子,同源多聚体蛋白还可以只保留一条链。直接从 pdb 下载蛋白可以
使用 fetch [protein name]命令。导入后 pdb 文件后,点击右下角的 S
以显示结构信息、

2.1.2 关于 pdb 文件的格式

该格式省略了一切氢原子,包括游离的水分子都只保留了一个 O,所以后续需要加氢操作。该文件的格式较复杂,最好使用成熟的软件进行编辑,而不上自行编辑。该文件中的氨基酸残基的原子全部记录在 ATOM 行中,后面的每一项分别为原子序号、原子名称(第一个字母为原子的元素符号,第二个字母为远近标识符 A、B、G、D、E、Z、H 分别对应有机化合物命名系统中的 α、β、γ、δ、ε、ζ、η)、残基名称、链编号、残基序号、原子坐标等

而离子、水分子以及结合在蛋白中的配体(如抑制剂等)等非蛋白质的部分记录在 HETNAM(非标准残基的名称)中:

同样的依次是原子编号、原子名称、基团名称(如水这里起名为HOH)、链编号(但是它本身不依附于哪条链)、原子编号、坐标等。

2.1.3 使用 PyMol 命令进行选择和删除
选择也可以使用命令选择,select 命令可以用于选择链如:
select chain A
indicate 命令的格式则为 indicate [element_type] [name],如根据残基
名称 HOH 选择所有的水分子:
indicate resname HOH
需要注意的是,选择水分子不要使用:
indicate name O
否则其会选中所有氨基酸残基中的氧原子。根据残基名称选择配体
Mr-Greyfun
(示例中配体名称为 STI):
indicate resname STI
根据原子名称选择氯原子:
indicate name CL
可以使用 remove 命令进行删除操作,命令格式与之类似:
remove resname HOH
remove name CL
remove chain B

2.1.1 保存
在 PyMol 中完成上述处理后使用 Export Molecule-save 保存成pdb 格式的文件。
2.1.2 加氢操作
pdb 文件会省略氢原子,所以需要进行加氢,可以在 PyMol 中使用 h_add 命 令 加氢 (删除为 remove hydrogen ), 但最 好 使 用AutoDockTools(ADT)进行加氢。处理好后的文件用 ADT 打开后,点击 Edit-Hydrogens-Add-OK。加氢。
2.1.3 转化为 pdbqt 格式的文件
点击 Grid-Macromolecule-Choose,选择好后点 select molecule,
它会自动加电荷等,点 OK 即可,完成后,会弹出保存窗口,直接保
存 pdbqt 格式即可

2.1.1 保存
在 PyMol 中完成上述处理后使用 Export Molecule-save 保存成
pdb 格式的文件。
2.1.2 加氢操作
pdb 文件会省略氢原子,所以需要进行加氢,可以在 PyMol 中使
用 h_add 命 令 加氢 (删除为 remove hydrogen ), 但最 好 使 用
AutoDockTools(ADT)进行加氢。
处理好后的文件用 ADT 打开后,点击 Edit-Hydrogens-Add-OK。
加氢。
2.1.3 转化为 pdbqt 格式的文件
点击 Grid-Macromolecule-Choose,选择好后点 select molecule,
它会自动加电荷等,点 OK 即可,完成后,会弹出保存窗口,直接保
存 pdbqt 格式即可

设置好后点 Done,然后点击 Ligand-Output-Save as PDBQT 即可
保存。

2.2.1 附:用 PyMol 导出 pdb 文件中的配体的方法
在界面上选中小分子,点击(sele)的 A-copy to object-new。红框内
的区域可以点击后将其对应的对象隐藏起来:
然后点击 Export Molecule-save 保存成 pdb 格式的文件。注意,这个
内置的配体也是没有加氢的,需要用 ADT 加氢。

3 对接
3.1 对接操作步骤
3.1.1 导入受体和配体
在 ADT 里面,点击 Grid-Macromolecule-Open 打开受体蛋白
点击 Grid-Set-Map-Types-Open Ligand 打开受体的 pdbqt 文件
3.1.2 定义对接盒子
点击 Grid-Grid Box 定义对接盒子
3.1.3 生成 config.txt 文件
点击 Docking-Output-Vina config 生成 config.txt 文件
3.1.4 启动对接

视频中的步骤:

1、导入protein蛋白质

2、删除水分子 如果你的对接区域有水分子,会影响对接结果

pdq该格式省略了一切氢原子,包括游离的水分子都只保留了一个 O,所以后续需要加氢操作。

3、edit ->hydrogens->add->polar only 此时结构中发亮的就是氢键(加氢)

4、加电荷 edit -> charges->add kollman charges

5 报存 grid -> macromolecule ->choose->select

准备配体文件:

以sdf结尾的文件直接拖进 AUTO软件中会报错,需要转换成pdb文件

可以使用pymol可视化工具转换(注意:配体 英文 ligend)

1、将该文件拖动到pymol中打开,file ->molecule

配体文件:

1、打开auto dock,将配体文件导入:

2、 点击 Ligand(配体)->input->choose (这一步就是生成了配体文件)

3、点击 Ligand(配体)->output->save as pdbqt

接下来就可以进行分子对接:

1、将两个(受体和配体)pdpqt导入

2、重新选择蛋白质分子作为受体

点击NO

接下来设置对接盒子:

grid -> grid box

设置盒子位置(spacing设置为1)

然后 grid box 弹出设置中选择 file->output grid dimensions file(保存盒子设置)

保存:

新建config.txt:用于启动vina

receptor 蛋白质名(生成的蛋白质文件名)

ligand 配体名

center 和size在上一步的grid中有

receptor:指定受体分子的路径

ligand:配体分子的路径

center_x,center_y,center_z:搜索空间中心的坐标

size_x,size_y,size_z:指定搜索空间的大小。这里设置的大小基本就是把整个受体分子都包含了,属于blind docking。如何更准确确定结合口袋的位置,我们稍后再说。

energy_range:默认4,与最优结合模型相差的最大能量值,单位是kcal/mol。比方说,最优模型的能量计算出来是-8.5kcal/mol,那么vina也就最多计算到-4.5kcal/mol的模型就终止了,也就意味着这个值决定了生成模型的最大个数。

exhaustiveness:用来控制对接的细致程度,默认值是8. 大致与时间成正比。

num_modes:最多生成多少个模型。实际生成的模型数由num_modes和energy_range共同决定。

energy_range
 maximum energy difference between the best binding 
 mode and the worst one displayed (kcal/mol)

最后启动vina 进行分子对接:

在config目录中打开cmd->输入vina

其中:这条命令就是利用config.txt文件进行分子对接(cmd 必须在config文件目录下打开)

执行 命令:

D:\Autodock\pdbqt>”D:\Vina\vina.exe” –receptor protein.pdbqt –ligand ligend.pdbqt –config config.txt –log log.txt –out output.pdbqt

“D:\Vina\vina.exe” –receptor selected_prediction_ready.pdbqt –ligand Conformer3D_CID_65536.pdbqt –config config.txt –log log.txt –out output.pdbqt

如果上述命令报错,一般是生成的config文件有问题:

正确的config内容如下:

执行完毕可以看到生成一个log文件

vina 中的affinity是亲合力结果排名

最后:

查看生成的模型:

将protein和output输出导入pymol

点击 all -> s查看表面 结构

点击左右箭头查看不同的结构:

5.我们挑选第一个模型,看看结构方式是怎样的,见下图。很显然与真实的结合方式相差甚远,可以说是完全错误。

-1.jpg

为什么会出现这种情况,很大程度上是因为search space太大,可能需要设置更大的exhaustiveness。

如果我们大致知道binding pocket在什么位置,那准确性应该会高不少,如何大致确定binding pocket的位置呢?我们接着试验。

可以通过实验的方式,比如某个点突变对结合或者活性影响非常大,那么大概率这个残基是结合口袋的一部分。

可以通过软件预测,比如蛋白与配体(底物)结合位点预测:https://zhanglab.ccmb.med.umich.edu/COACH/

再比如Discovery studio软件(专业版的),可以很方便的根据受体分子的表面形状来预测结合口袋位置。

20190507175152.png

口袋的坐标为:

34.3356,14.9412,26.9615

我们修改下对接参数,新的参数如下:

receptor = r.pdbqt

ligand = nap.pdbqt

center_x = 34.3356

center_y = 14.9412

center_z = 26.9615

size_x = 30.0

size_y = 30.0

size_z = 30.0

energy_range = 4

exhaustiveness = 10

num_modes = 10

最优结果与真实模型的RMSD为1.795埃,可以说非常精准了

比对结果如下:

Souce: 纽普生物    2019-05-07

PyMOL 相关操作:

1、导入蛋白质

先在NCBI子数据库structure检索所需要的蛋白结构,https://www.ncbi.nlm.nih.gov/structure/?term=
记录下该蛋白的PDB ID,然后打开PyMOL,命令行输入:

fetch 5ocn#foxn1
fetch 3uf0#
fetch 1si4#血红蛋白

去除水分子

remove solvent

分离得到蛋白

remove organic

分离配体:

../_images/splitlig2020-03-12_202822.573956.png

pymol教程:

http://pymol.chenzhaoqiang.com/intro/startManual.html

分子对接进阶教程:选定对接位点区域:

如果我们想要将配体和蛋白质受体的某几个位点部位对接:

1、标记这些点 select->select from string ,选择对应的molecule和链,以及residue(寻找的位点),点击add

2、此时,会出现currenr selection ,选择该列中 c,点击绿点,会将对应的分子三D化。

3、接下来设置盒子,spacing 设置为1(相当于比例尺),xyz一般设置为20-22

这样就可以了。