2021-07-28 论文复现-ACS Catal. 2018, 8, 3, 2402–2412(6)

-为啥么要死磕 hierarchical clustering 呢？直接用 Jalview 不香吗？
-希望做到可视化与算法互为表里，不仅仅是软件呈现一个结果，还希望可以去利用结果设计一些算法。
-那还不如去看看(1) Jalview里面Tree的设置和导出 (2) scikit-learn如何做hierarchical clustering。
-那也先把biopython的cluster部分看完吧。要不太没有信心了。
-好吧。

Biopython

15.4 hierarchical clustering

dendrogram 树状图
可视化方法，文中提到了 Java TreeView，可以尝试一下。

https://sourceforge.net/projects/jtreeview/

原文是从 tree 的构造、方法来讲，最后才讲了如何从数据来构建 tree。在此，我们从数据出发，正着来走教程。

from Bio.Cluster import treecluster
tree = treecluster(data) # 从数据，如基因表达的 array 数据出发.

'''
可以定义的参数：datamaskweighttransposemethod : {'s','m','c','a'}
'''# 或者直接使用distance数据进行计算的方法
tree = treecluster(distancematrix=distance)
'''
可以定义的参数：data: 2D Numerical Python array. 只需要矩阵左下角的数据：distance = array([[0.0, 1.1, 2.3], [1.1, 0.0, 4.5], [2.3, 4.5, 0.0]])distance = array([1.1, 2.3, 4.5])distance = [array([]), array([1.1]), array([2.3, 4.5])]method : {'s','m','a'}
tree 将是 Bio.Cluster.Tree 对象。
'''

测试了 distancematrix 和 method 输入。

import numpy as np
'== 2D-array初始化 =='
distance1 = np.array([[0.0, 1.1, 2.3], [1.1, 0.0, 4.5], [2.3, 4.5, 0.0]])
# array([[0. , 1.1, 2.3],
#        [1.1, 0. , 4.5],
#        [2.3, 4.5, 0. ]])
tree1 = treecluster(data=None,distancematrix=distance1) # 默认 method = max
print(tree1)
# (1, 0): 1.1
# (2, -1): 4.5
'tree的节点，正数表示叶子节点(原始数据点)，负数表示subnode'
distance1
# array([[0. , 1.1, 2.3],
#        [4.5, 0. , 4.5],
#        [4.5, 4.5, 0. ]])
'distance在计算后被改变，需要注意''== 1D-array初始化 =='
distance4 = np.array([1.1, 2.3, 4.5])
# array([[0. , 1.1, 2.3],
#        [1.1, 0. , 4.5],
#        [2.3, 4.5, 0. ]])
tree4 = treecluster(data=None,distancematrix=distance4)
print(tree4)
# (1, 0): 1.1
# (2, -1): 4.5
distance4
# array([4.5, 4.5, 4.5])'==下三角 2D-array初始化 =='
distance5 = np.array([[0.0, 0.0, 0.0], [1.1, 0.0, 0.0], [2.3, 4.5, 0.0]])
# array([[0. , 1.1, 2.3],
#        [1.1, 0. , 4.5],
#        [2.3, 4.5, 0. ]])
tree5 = treecluster(data=None,distancematrix=distance5)
print(tree5)
# (1, 0): 1.1
# (2, -1): 4.5
distance5
# array([[0. , 0. , 0. ],
#        [4.5, 0. , 0. ],
#        [4.5, 4.5, 0. ]])'== method = average =='
distance2 = np.array([[0.0, 1.1, 2.3], [1.1, 0.0, 4.5], [2.3, 4.5, 0.0]])
tree2 = treecluster(data=None,method='a',distancematrix=distance2)
print(tree2)
# (1, 0): 1.1
# (2, -1): 3.4
distance2
# array([[0. , 1.1, 2.3],
#        [3.4, 0. , 4.5],
#        [3.4, 4.5, 0. ]])'== method = min =='
distance3 = np.array([[0.0, 1.1, 2.3], [1.1, 0.0, 4.5], [2.3, 4.5, 0.0]])
tree3 = treecluster(data=None,method='s',distancematrix=distance3)
print(tree3)
# (1, 0): 1.1
# (-1, 2): 2.3
distance3
# array([[0. , 1.1, 2.3],
#        [1.1, 0. , 4.5],
#        [2.3, 4.5, 0. ]])

Bio.Cluster.Tree 对象有一些适用的方法：

1. tree.cut(nclusters=i)将叶子节点分为i个类别，返回的numpy.ndarray表明每个节点归属的类别。

tree = tree1
print(tree) # Tree是一个 imutable 对象，可以通过一些方式进行调整。
# (1, 0): 1.1
# (2, -1): 4.5
tree[0] # 返回 Node 对象
# (0, 1): 1.1
tree[1]
# (2, -1): 4.5
indices1 = tree.cut(nclusters=1)
indices1
# array([0, 0, 0], dtype=int32)
indices2 = tree.cut(nclusters=2)
indices2
# array([1, 1, 0], dtype=int32)
indices3 = tree.cut(nclusters=3)
indices3
# array([2, 1, 0], dtype=int32)
indices4 = tree.cut(nclusters=4)
# 报错：ValueError: more clusters requested than items available

2. tree.sort()
   一个树有2^(n-1)种组织方式，通过sort可以固定到某一种。

3. tree.scale()
   某些可视化软件，如 Java TreeView 需要树的距离在0-1之间，通过scale函数，使得结果满足可视化软件要求。

15.7 Handling Cluster/TreeView-type files

Alok Saldanha’s http://jtreeview.sourceforge.net/Java TreeView program, which can display hierarchical as well as k-means clustering results.

实际上包装了一个record，整合了数据输入，数据处理，数据输出过程。但对于distancematrix直接输入输出的方式不太适用。

from Bio import Cluster
with open("cyano.txt") as handle:record = Cluster.read(handle)genetree = record.treecluster(method="s")
genetree.scale()
exptree = record.treecluster(dist="u", transpose=1)
record.save("cyano_result", genetree, exptree) # 保存的结果用于 Java Treeview 的输入。

fasta > USEARCH blast6out > Biopython

## 数据载入
# 载入 fasta数据，获取annotation
from Bio import SeqIO
import numpy as np
import pandas as pdfilename = 'your_self_filepath'
records = list(SeqIO.parse(filename, "fasta"))   
n_sequence = len(records) 
records_dict = {records[i].description:i for i in range(len(records))} # 构造字典{annotation->i}file_blast6out = 'your_self_filepath'
blast6out = pd.read_table(file_blast6out, header = None) # 载入 blast6out数据，构造下三角矩阵
list_i = [ records_dict[des] for des in blast6out[0]] # 可以用pd的apply或者map函数，更紧凑一些
list_j = [ records_dict[des] for des in blast6out[1]]
list_s = [ 1 - score/100.0 for score in blast6out[2]]distance_matrix = np.zeros((n_sequence, n_sequence))
for i,j,s in zip(list_i,list_j,list_s):distance_matrix[i,j] = sdistance_matrix[j,i] = s## 数据处理 - hierarchical clustering
from Bio.Cluster import treecluster
tree = treecluster(data=None,distancematrix=distance_matrix.copy()) # 希望不改变原数据，因此传入复制
indices = tree.cut(nclusters=4) # 4是随意的一个数字
indices## 数据可视化(略)
# 能直接调包是最好，但是如果要自己写轮子那就。。。
# 所以还是先看 JalView 吧。## 数据分析
# 1. 例如，根据某种准则确定cut。
# 2. 例如，确定某一个序列在某一种cut中位于哪一个cluster中，并且得到与之相邻的cluster。
# 3. 例如，确定某几条已知序列，筛选cut，然后通过heatmap来可视化不同cut下的分类情况，然后根据某种准则确定最好的cut情况。
# 既然如此，为什么要把所有blast得到的结果均做聚类呢？
# 4. 根据依赖于相似性的某种准则输出若干条备选序列。

JalView 中的 Tree

树的构建

File -> Input Alignment -> From File, 选择 100.txt。这是一个fasta文件。
在打开的窗口里 Web Service -> Alignment -> Clustal -> With default。
任务结束后出现新的窗口，关闭原窗口和任务提交窗口，外层 File -> Save Project as … ，保存为 test-100.jvp
Calculate -> Calculate Tree of PCA … > Tree
聚类算法可以选择 Neighbor Joining 或者 Average Distance (UPGMA).
度量方法可以选择 BLOSUM62, PAM250, PID (% Identity), Sequence Feature Similarity.

BLOSUM62和PAM250是从进化的角度看待序列，相关概念任何一门生物信息学课程都会涉及。
在选择上，有这么个建议：

PAM	BLOSUM
To compare closely related sequences, PAM matrices with lower numbers are created.	To compare closely related sequences, BLOSUM matrices with higher numbers are created.
To compare distantly related proteins, PAM matrices with high numbers are created.	To compare distantly related proteins, BLOSUM matrices with low numbers are created.

https://en.wikipedia.org/wiki/BLOSUM
也就是，更相关的序列用更高的BLOSUM，更不相关的序列用更高的PAM。
在这里，更相关的序列用BLOSUM62(一个中等的)，更不相关的序列用PAM250。

结果可以导出为Newick 文件。得到的是一行文件，用括号匹配和数字表示了树的分割以及distance。

Tree Window 是可以交互的，可以对树进行分割，得到的各个 Cluster 被划分为不同的颜色。通过 Calculate -> Sort -> By Tree Order -> … 可以将同一个 Cluster 的序列排布到一起，结合 MSA 的颜色看，相当直观。

Average Distance Using PID

计算也可以针对部分行和部分列进行。

下图针对pentad进行树的构建。> Neighbour Joining Using BLOSUM62, 41,111,112,135,275
可以确定这100条序列中三条其它的序列：

WP_110315832.1	NDWEH
WP_003910922.1	NDWDH
MBV8838173.1	NDWEY
AEO24280.1	NDWE-

可以从树的页面恢复计算该树时的数据来源。 File -> Input Data
可以选择与树(颜色)相链接的 View 的页面。 View -> Associated Nodes With -> …

Tree Based Conservation Analysis

在Tree Window中分组后，在Alignment Window中通过
Annotations -> Autocalculation -> Group Conservation, Group Consensus 可以得到对于序列保守性的分析。

Redundancy Removal

可以移除PID大于某一阈值的序列。注意对所有的view是同时生效的。
Edit -> Remove Redundancy …

根据若干位点的残基进行分组 Subdividing the Alignment According to Specific Mutations

若是为了考察某几个点的突变，相比于建树，更简单的方法是选中这几列，然后用Select -> Make groups for selection， 然后 Calculate -> Sort -> By Group。然后再用 Annotations -> Autocalculated Annotation -> Group Consensus 来考察。

其它

做生信分析的话为什么不去搞R语言呢…大无语。

scipy dendrogram

https://docs.scipy.org/doc/scipy/reference/generated/scipy.cluster.hierarchy.dendrogram.html