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是最长序列长度，B是batch 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 下面都是 0，5-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

这样就可以了。

Python中collections.defaultdict()使用

1、以一个例子开始：统计一个列表里各单词重复次数

words = ['hello', 'world', 'nice', 'world']
counter = dict()
for kw in words:
    counter[kw] += 1

这样写肯定会报错的，因为各词的个数都没有初始值，引发KeyError

2、改进一下：加入if判断

words = ['hello', 'world', 'nice', 'world']
counter = dict()
for kw in words:
    if kw in words:
        counter[kw] += 1
    else:
        counter[kw] = 0

3、再改进：使用setdefault()方法设置默认值

words = ['hello', 'world', 'nice', 'world']
counter = dict()
for kw in words:
    counter.setdefault(kw, 0)
    counter[kw] += 1

setdefault()，需提供两个参数，第一个参数是键值，第二个参数是默认值，每次调用都有一个返回值，如果字典中不存在该键则返回默认值，如果存在该键则返回该值，利用返回值可再次修改代码。

4、接着改进
一种特殊类型的字典本身就保存了默认值defaultdict()，defaultdict类的初始化函数接受一个类型作为参数，当所访问的键不存在的时候，可以实例化一个值作为默认值。

from collections import defaultdict
dd = defaultdict(list)
defaultdict(<type 'list'>, {})
#给它赋值，同时也是读取
dd['hh']
defaultdict(<type 'list'>, {'hh': []})
dd['hh'].append('haha')
defaultdict(<type 'list'>, {'hh': ['haha']})

该类除了接受类型名称作为初始化函数的参数之外，还可以使用任何不带参数的可调用函数，到时该函数的返回结果作为默认值，这样使得默认值的取值更加灵活。

>>> from collections import defaultdict
>>> def zero():
...     return 0
...
>>> dd = defaultdict(zero)
>>> dd
defaultdict(<function zero at 0xb7ed2684>, {})
>>> dd['foo']
0
>>> dd
defaultdict(<function zero at 0xb7ed2684>, {'foo': 0})

最终代码：

from collections import defaultdict
words = ['hello', 'world', 'nice', 'world']
#使用lambda来定义简单的函数：函数返回值就是默认key的value值
counter = defaultdict(lambda: 0) 
for kw in words:
    counter[kw] += 1

python 字典树、前缀匹配

在python中可以用字典来实现字典树这一结构。

class Trie:

    def __init__(self):
        """
        Initialize your data structure here.
        """
        self.root = {}
        self.isEnd = '#'

    def insert(self, word: str) -> None:
        """
        Inserts a word into the trie.
        """
        node = self.root
        for char in word:
            if char not in node:
                node[char] = {}  #如果没有这一字母则新建一个字典
            node = node[char]
        node[self.isEnd] = '#'   #加上结束符

    def search(self, word: str) -> bool:
        """
        Returns if the word is in the trie.
        """
        node = self.root
        for char in word:
            if char not in node:
                return False
            node = node[char]
        return self.isEnd in node  #如果没有# 说明只是前缀


    def startsWith(self, prefix: str) -> bool:
        """
        Returns if there is any word in the trie that starts with the given prefix.
        """
        node = self.root
        for char in prefix:
            if char not in node:
                return False
            node = node[char]
        return True



# Your Trie object will be instantiated and called as such:
# obj = Trie()
# obj.insert(word)
# param_2 = obj.search(word)
# param_3 = obj.startsWith(prefix)

2.采用defaultdict创建trie



from collections import defaultdict
from functools import reduce
TrieNode = lambda: defaultdict(TrieNode)
class Trie:
    def __init__(self):
        self.trie = TrieNode()
    def insert(self, word):
        reduce(dict.__getitem__, word, self.trie)['end'] = True
    def search(self, word):
        return reduce(lambda d,k: d[k] if k in d else TrieNode(), word, self.trie).get('end', False)
    def startsWith(self, word):
        return bool(reduce(lambda d,k: d[k] if k in d else TrieNode(), word, self.trie).keys())

collections，解释是数据类型容器模块。这里面有一个collections.defaultdict()经常被用到

这里的defaultdict(function_factory)构建的是一个类似dictionary的对象，其中keys的值，自行确定赋值，但是values的类型，是function_factory的类实例，而且具有默认值。比如default(int)则创建一个类似dictionary对象，里面任何的values都是int的实例，而且就算是一个不存在的key, d[key] 也有一个默认值，这个默认值是int()的默认值0.

python 二分查找 ：bisect模块

bisect模块的源码：

"""Bisection algorithms."""

def insort_right(a, x, lo=0, hi=None):

    if lo < 0:
        raise ValueError('lo must be non-negative')
    if hi is None:
        hi = len(a)
    while lo < hi:
        mid = (lo+hi)//2
        if x < a[mid]: hi = mid
        else: lo = mid+1
    a.insert(lo, x)

insort = insort_right   # backward compatibility

def bisect_right(a, x, lo=0, hi=None):

    if lo < 0:
        raise ValueError('lo must be non-negative')
    if hi is None:
        hi = len(a)
    while lo < hi:
        mid = (lo+hi)//2
        if x < a[mid]: hi = mid
        else: lo = mid+1
    return lo

bisect = bisect_right   # backward compatibility

def insort_left(a, x, lo=0, hi=None):

    if lo < 0:
        raise ValueError('lo must be non-negative')
    if hi is None:
        hi = len(a)
    while lo < hi:
        mid = (lo+hi)//2
        if a[mid] < x: lo = mid+1
        else: hi = mid
    a.insert(lo, x)


def bisect_left(a, x, lo=0, hi=None):

    if lo < 0:
        raise ValueError('lo must be non-negative')
    if hi is None:
        hi = len(a)
    while lo < hi:
        mid = (lo+hi)//2
        if a[mid] < x: lo = mid+1
        else: hi = mid
    return lo

# Overwrite above definitions with a fast C implementation
try:
    from _bisect import *
except ImportError:
    pass

方法介绍

bisect模块采用经典的二分算法查找元素。模块提供下面几个方法：

bisect.bisect_left(a, x, lo=0, hi=len(a))：返回待插入值应该放置的位置下标，如果 a 中有和x相同的元素，就在返回该元素左边的下标，即在左边插入

定位x在序列a中的插入点，并保持原来的有序状态不变。参数lo和hi用于指定查找区间。如果x已经存在于a中，那么插入点在已存在元素的左边。函数的返回值是列表的整数下标。

bisect.bisect_right(a, x, lo=0, hi=len(a))

和上面的区别是插入点在右边。函数的返回值依然是一个列表下标整数。

bisect.bisect(a, x, lo=0, hi=len(a))

等同于bisect_right()。

注意，前面这三个方法都是获取插入位置，也就是列表的某个下标的，并不实际插入。但下面三个方法实际插入元素，没有返回值。

bisect.insort_left(a, x, lo=0, hi=len(a))

将x插入有序序列a内，并保持有序状态。相当于a.insert(bisect.bisect_left(a, x, lo, hi), x)。碰到相同的元素则插入到元素的左边。这个操作通常是 O(log n) 的时间复杂度。

bisect.insort_right(a, x, lo=0, hi=len(a))

同上，只不过碰到相同的元素则插入到元素的右边。

bisect.insort(a, x, lo=0, hi=len(a))

等同于insort_right()。

单调栈

单调栈

单调栈实际上就是栈，只是限制要比普通的栈更严格而已了。要求是每次入栈的元素必须要有序（如果新元素入栈不符合要求，则将之前的元素出栈，直到符合要求再入栈），使之形成单调递增/单调递减的一个栈。

  • 单调递增栈：只有比栈顶小的才能入栈，否则就把栈顶出栈后，再入栈。出栈时可能会有一些计算。适用于求解第一个大于该位置元素的数。
  • 单调递减栈：与单调递增栈相反。适用于求解第一个小于该位置元素的数。

如何判断

单调递增/递减栈是根据出栈后的顺序来决定的。例如，栈内顺序[1, 2, 6]，出栈后顺序[6, 2, 1]，这就是单调递减栈。

适用场景

单调栈一般用于解决 第一个大于 xxx 或者 第一个小于 xxx 这种题目。

int[] ans = new int[T.Length];
Stack<int> stack = new Stack<int>();

for(int i = 0; i < T.Length; i++) {
    while(stack.Count > 0 && T[i] > T[stack.Peek()]) {
        int curI = stack.Pop();
        ans[curI] = i - curI;
    }
    stack.Push(i);
}