PCoA全称为Principal Coordinates Analysis，是一种常用的多元分析方法。在生物信息学中，PCoA通常用于探索样本之间的相似性和差异性，比如基于OTU表格或者不同生物数据分析后得到的距离矩阵，通过PCoA方法将高维的样本数据映射成二维或三维空间中进行可视化。本文将对PCoA分析方法进行详细介绍，包括距离矩阵计算、PCoA主成分计算、坐标生成和可视化等方面。
一、距离矩阵的计算
距离矩阵是PCoA分析的关键输入，通常依据不同的数据类型和目的选择合适的距离度量方法。其中常用的距离度量方法包括曼哈顿距离、欧几里得距离、皮尔逊相关系数、Jaccard相似性系数等。下面以OTU表格为例演示如何基于欧式距离计算距离矩阵：
from scipy.spatial.distance import pdist, squareform
import pandas as pd
otu_table = pd.read_csv('otu_table.csv')
otu_table = otu_table.set_index('OTUs')
# 计算欧式距离
distance_matrix = pdist(otu_table.values, metric='euclidean')
# 将距离向量转化为距离矩阵
distance_matrix = squareform(distance_matrix)
# 将距离矩阵保存为csv文件
pd.DataFrame(distance_matrix, index=otu_table.index, columns=otu_table.index).to_csv('distance_matrix.csv')
二、PCoA主成分计算
距离矩阵一般为对称矩阵，包含n个样本之间的距离值。在进行PCoA分析前，需要对这个距离矩阵进行特征值分解（eigendecomposition）和奇异值分解（singular value decomposition，SVD），计算出其前k个主成分（principal components），其中k是分析者根据数据的特征和目的自行拟定的参数。下面以scipy库为例，展示如何进行PCoA主成分计算：
from scipy.spatial import distance_matrix
from scipy.linalg import eigh
# 读取距离矩阵
distance_matrix = distance_matrix.read_csv('distance_matrix.csv', index_col=0)
# 计算内积矩阵
double_centering_matrix = -0.5 * (distance_matrix ** 2 - \
(distance_matrix ** 2).mean(axis=0)[:, np.newaxis] - \
(distance_matrix ** 2).mean(axis=1)[:, np.newaxis] + \
(distance_matrix ** 2).mean())
# 特征值分解
eigenvalues, eigenvectors = eigh(double_centering_matrix)
# 获取前k个主成分
k = 3
principal_components = eigenvectors[:, -k:]
三、坐标生成与可视化
在进行了PCoA主成分计算后，我们可以将k个主成分对应的坐标生成为一个k维矩阵，进行可视化展示。其中，生成的坐标可以用于不同的生物数据的可视化，比如16S rRNA数据、代谢组学数据和基因表达谱数据等。下面以matplotlib库为例，展示如何进行坐标生成和可视化：
import matplotlib.pyplot as plt
# 读取距离矩阵
distance_matrix = pd.read_csv('distance_matrix.csv', index_col=0)
# 计算内积矩阵
double_centering_matrix = -0.5 * (distance_matrix ** 2 - \
(distance_matrix ** 2).mean(axis=0)[:, np.newaxis] - \
(distance_matrix ** 2).mean(axis=1)[:, np.newaxis] + \
(distance_matrix ** 2).mean())
# 特征值分解
eigenvalues, eigenvectors = eigh(double_centering_matrix)
# 获取前k个主成分
k = 3
principal_components = eigenvectors[:, -k:]
# 绘制二维或三维图
if k == 2:
plt.scatter(principal_components[:, 0], principal_components[:, 1])
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.show()
elif k == 3:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(principal_components[:, 0], principal_components[:, 1], principal_components[:, 2])
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_zlabel('PC3')
plt.show()
四、聚类分析与差异分析
除了PCoA主成分分析外，还可以结合聚类分析和差异分析等方法，更好地探索样本之间的相似性和差异性。其中，聚类分析可以利用距离矩阵将样本聚成几个类，不同类别间的样本具有明显差异性；差异分析则可以找出不同类别间显著差异的OTUs或生物差异标记，进行生物学解释。下面以sklearn库为例，展示如何进行层次聚类分析和差异分析：
# 读取距离矩阵
distance_matrix = pd.read_csv('distance_matrix.csv', index_col=0)
# 层次聚类分析
from sklearn.cluster import AgglomerativeClustering
cluster = AgglomerativeClustering(n_clusters=3, affinity='euclidean', linkage='ward')
cluster_labels = cluster.fit_predict(distance_matrix)
# 差异分析
from sklearn.feature_selection import f_classif
import numpy as np
# 读取OTU表格和分类标签
otu_table = pd.read_csv('otu_table.csv').set_index('OTUs')
metadata = pd.read_csv('metadata.csv').set_index('SampleID')
# 过滤低频OTUs
otu_count = (otu_table > 0).sum(axis=0)
otu_table = otu_table.loc[:, otu_count >= 5]
# 计算差异分析p值和FDR校正p值
pval, _, _, _ = f_classif(otu_table.T, metadata['Group'])
pval = pd.Series(pval, index=otu_table.columns)
pval_corrected = pd.Series(statsmodels.stats.multitest.fdrcorrection(pval)[1], index=pval.index)
# 可视化聚类结果和差异分析结果
import seaborn as sns
sns.clustermap(otu_table, row_cluster=False, col_cluster=False)
sns.catplot(x='Group', y='pval_corrected', data=pd.concat([pval_corrected, metadata['Group']], axis=1).melt(id_vars='Group'), kind='box')
