MATLAB和Python及R聚类和亚群识别

MATLAB片段

在 MATLAB 中进行聚类和亚群识别可以使用一些内置函数和工具箱，如 Statistics and Machine Learning Toolbox 和 Bioinformatics Toolbox。以下是如何使用 MATLAB 进行聚类和亚群识别的详细步骤，尤其是用于检测和分析选择性扰动效应。

1. 数据导入和预处理

首先，需要导入数据并进行预处理。假设数据是一个表达矩阵，行是基因，列是样本。

% 导入表达数据
expressionData = readmatrix('expression_data.csv'); % 行表示基因，列表示样本% 可选：导入样本的元数据（如扰动标签）
metadata = readtable('metadata.csv'); % 假设包含亚群和扰动标签

2. 数据标准化和降维

使用 zscore 进行标准化处理，并使用 PCA 或 t-SNE 进行降维，以便进行聚类。

% 对数据进行标准化
standardizedData = zscore(expressionData, 0, 2);% 进行 PCA 降维
[coeff, score, ~] = pca(standardizedData');% 使用前两个主成分进行可视化
figure;
scatter(score(:,1), score(:,2), 20, 'filled');
title('PCA Visualization of Samples');
xlabel('PC1');
ylabel('PC2');

3. 聚类分析

可以使用 k-means 或 层次聚类 等方法对样本进行聚类。

% 使用 k-means 聚类
numClusters = 3; % 假设我们预期有 3 个亚群
[idx, C] = kmeans(score(:, 1:10), numClusters); % 使用前 10 个主成分% 可视化聚类结果
figure;
gscatter(score(:,1), score(:,2), idx);
title('k-means Clustering of Samples');
xlabel('PC1');
ylabel('PC2');

4. 层次聚类

层次聚类可以用于更详细的分析，尤其适用于树状图。

% 计算样本之间的欧几里得距离
distances = pdist(standardizedData');% 进行层次聚类
linkageTree = linkage(distances, 'ward');% 绘制树状图
figure;
dendrogram(linkageTree, 0);
title('Hierarchical Clustering Dendrogram');

5. 亚群标注和分析

根据聚类结果，将样本分配到不同的亚群，并使用元数据进行标注和扰动效应分析。

% 将聚类结果添加到元数据中
metadata.Cluster = idx;% 使用散点图查看扰动效应在亚群中的分布
figure;
scatter(score(:,1), score(:,2), 20, metadata.Perturbation, 'filled');
title('Perturbation Distribution Across Clusters');
xlabel('PC1');
ylabel('PC2');

6. 差异表达分析（可选）

可以使用 t 检验或其他方法来分析不同亚群之间的差异表达情况。

% 初始化结果存储
pValues = zeros(size(expressionData, 1), numClusters);% 对每个亚群进行差异表达分析
for cluster = 1:numClustersclusterSamples = standardizedData(:, metadata.Cluster == cluster);otherSamples = standardizedData(:, metadata.Cluster ~= cluster);for i = 1:size(expressionData, 1)[~, pValues(i, cluster)] = ttest2(clusterSamples(i, :), otherSamples(i, :));end
end% 应用多重比较校正（Benjamini-Hochberg 法）
adjustedPValues = mafdr(pValues(:), 'BHFDR', true);

7. 结果可视化

绘制热图以显示不同亚群中显著差异表达的基因。

% 选择显著差异表达的基因
significantGenes = find(adjustedPValues < 0.05);% 绘制热图
figure;
heatmap(standardizedData(significantGenes, :));
title('Heatmap of Differentially Expressed Genes');
xlabel('Samples');
ylabel('Significant Genes');

8. 下游分析和解释

根据聚类和差异表达分析的结果，解释亚群的生物学意义，特别是扰动效应在这些亚群中的表现。可以使用外部工具（如 Enrichr 或 DAVID）进行基因富集分析。

这个流程提供了如何使用 MATLAB 进行聚类和亚群识别，以及如何检测和分析选择性扰动效应的完整指南。你可以根据需要调整聚类参数和分析方法，以便更好地适应具体的数据和研究问题。

Python片段

在 Python 中，进行聚类和亚群识别以检测和分析选择性扰动效应可以使用诸如 scikit-learn、pandas、numpy、seaborn 和 matplotlib 等库。以下是一个详细的流程和代码示例，展示如何使用这些工具进行此类分析。

步骤概述：

导入数据和库。
数据预处理（如标准化）。
降维（使用 PCA 或 t-SNE）。
聚类分析（使用 KMeans 或层次聚类）。
结果可视化。
亚群标注和差异分析。
绘制热图和分析扰动效应。

代码示例：

# 导入库
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.cluster import KMeans
from sklearn.manifold import TSNE
from scipy.cluster.hierarchy import linkage, dendrogram, fcluster# 读取数据
expression_data = pd.read_csv('expression_data.csv', index_col=0)  # 行是基因，列是样本
metadata = pd.read_csv('metadata.csv', index_col=0)  # 包含样本的元数据，如亚群和扰动标签# 数据标准化
scaler = StandardScaler()
scaled_data = scaler.fit_transform(expression_data.T)  # 转置以按样本标准化# PCA 降维
pca = PCA(n_components=2)
pca_result = pca.fit_transform(scaled_data)
plt.figure(figsize=(10, 7))
plt.scatter(pca_result[:, 0], pca_result[:, 1], c='blue', alpha=0.6)
plt.title('PCA Visualization')
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.show()# t-SNE 降维可选项
tsne = TSNE(n_components=2, perplexity=30, random_state=42)
tsne_result = tsne.fit_transform(scaled_data)
plt.figure(figsize=(10, 7))
sns.scatterplot(x=tsne_result[:, 0], y=tsne_result[:, 1], hue=metadata['Group'], palette='viridis')
plt.title('t-SNE Visualization')
plt.show()

3. 聚类分析

使用 KMeans 聚类分析样本。

# KMeans 聚类
kmeans = KMeans(n_clusters=3, random_state=42)
kmeans_labels = kmeans.fit_predict(pca_result)# 将聚类结果添加到元数据中
metadata['Cluster'] = kmeans_labels# 可视化聚类结果
plt.figure(figsize=(10, 7))
sns.scatterplot(x=pca_result[:, 0], y=pca_result[:, 1], hue=metadata['Cluster'], palette='Set2')
plt.title('KMeans Clustering on PCA Result')
plt.show()

4. 层次聚类

# 计算欧几里得距离并进行层次聚类
linked = linkage(scaled_data, method='ward')
plt.figure(figsize=(12, 7))
dendrogram(linked, labels=metadata.index, leaf_rotation=90, leaf_font_size=8)
plt.title('Hierarchical Clustering Dendrogram')
plt.show()# 提取聚类标签
hier_labels = fcluster(linked, t=3, criterion='maxclust')  # 3 个聚类
metadata['HierCluster'] = hier_labels

5. 差异表达分析

使用 t 检验来分析不同亚群的差异表达。

from scipy.stats import ttest_ind# 比较两个聚类组中的基因表达
cluster1 = expression_data.loc[:, metadata['Cluster'] == 0]
cluster2 = expression_data.loc[:, metadata['Cluster'] == 1]p_values = []
for gene in expression_data.index:_, p = ttest_ind(cluster1.loc[gene], cluster2.loc[gene], equal_var=False)p_values.append(p)# 调整 P 值（多重比较校正）
adjusted_p_values = np.array(p_values) * len(p_values)  # Bonferroni 校正
significant_genes = expression_data.index[adjusted_p_values < 0.05]print("显著差异表达的基因数量:", len(significant_genes))

6. 热图可视化显著基因

# 使用显著差异基因绘制热图
sns.clustermap(expression_data.loc[significant_genes], metric='euclidean', cmap='vlag', col_cluster=False)
plt.title('Heatmap of Differentially Expressed Genes')
plt.show()

7. 分析扰动效应

可视化扰动效应在不同聚类中的分布。

plt.figure(figsize=(10, 7))
sns.scatterplot(x=pca_result[:, 0], y=pca_result[:, 1], hue=metadata['Perturbation'], palette='coolwarm')
plt.title('Perturbation Effects Across Clusters')
plt.show()