PCoA分析

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')