生物角度分析

1. 差异化表达分析
2. 细胞共表达
3. 聚类稳定性

import scanpy as sc
import pandas as pd

# 1. 加载数据（无行名/列名）
data_matrix = pd.read_csv("data.csv", header=None)
labels = pd.read_csv("pred_labels.csv", header=None, names=["cluster"])

# 2. 生成默认行名（假设细胞ID为cell1, cell2,...）
data_matrix.index = [f"cell{i+1}" for i in range(data_matrix.shape[0])]

# 3. 生成默认列名（假设基因为gene1, gene2,...）
data_matrix.columns = [f"gene{i+1}" for i in range(data_matrix.shape[1])]

# 4. 关联标签到数据矩阵
adata = sc.AnnData(X=data_matrix.values, obs=labels)
#5. 将 'cluster' 列转换为 category 类型
adata.obs['cluster'] = adata.obs['cluster'].astype('category')


# # 差异表达分析
# 数据预处理
sc.pp.filter_cells(adata, min_genes=200)    # 过滤低质量细胞
sc.pp.filter_genes(adata, min_cells=3)     # 过滤低表达基因
sc.pp.normalize_total(adata, target_sum=1e4)  # 归一化
sc.pp.log1p(adata)                         # 对数转换

# 差异表达分析
sc.tl.rank_genes_groups(adata, groupby="cluster", method="wilcoxon")

# 提取结果
de_genes = sc.get.rank_genes_groups_df(adata, group=None)
de_genes.to_csv("DE_genes.csv")


# 火山图
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 筛选显著差异基因（示例阈值：padj < 0.05 & log2FC > 1）
significant = de_genes[de_genes["pvals_adj"] < 0.05]
significant = significant[abs(significant["logfoldchanges"]) > 1]

fig=plt.figure(figsize=(10, 6))
sns.scatterplot(
    x=de_genes["logfoldchanges"],
    y=-np.log10(de_genes["pvals_adj"]),
    hue=de_genes.index.isin(significant.index),
    palette={True: "red", False: "gray"},
    alpha=0.6
)
plt.xlabel("log2(Fold Change)")
plt.ylabel("-log10(Adjusted p-value)")
plt.title("DE Genes Identified by  Clustering")
plt.legend(title="Significant", bbox_to_anchor=(1, 1))
# 保存图像为tif格式，dpi设置为600
fig.savefig('01.tif', format='tif', dpi=600)
plt.show()

#热图
# 提取前20个DE基因（按padj排序）
top_genes = de_genes.sort_values("padj").index[:20]
subset = adata[:, top_genes].to_df()

# 按簇分组并标准化表达量
subset["cluster"] = adata.obs["cluster"]
subset = subset.groupby("cluster").mean().T  # 计算簇平均表达
subset = (subset - subset.mean(axis=1)) / subset.std(axis=1)  # Z-score标准化

plt.figure(figsize=(12, 8))
sns.heatmap(subset, cmap="viridis", annot=False, fmt=".2f")
plt.title("Top DE Genes Across Clusters")
plt.xlabel("Clusters")
plt.ylabel("Genes")
plt.show()


#细胞共表达
corr_matrix = np.corrcoef(adata.X, rowvar=False)
high_var_genes = adata.var_names[np.argsort(adata.X.var(axis=0))[-50:]]
subset_corr = corr_matrix[np.ix_(adata.var_names.isin(high_var_genes), 
                                adata.var_names.isin(high_var_genes))]
fig=plt.figure(figsize=(12, 8))
plt.imshow(subset_corr, cmap='coolwarm', interpolation='none')
plt.colorbar()
plt.title("Gene Co-expression Network")
fig.savefig('01.tif', format='tif', dpi=600)
plt.show()

# 聚类稳定性
from sklearn.metrics import silhouette_score
silhouette_avg = silhouette_score(adata.X, adata.obs['cluster'])
print(f"Silhouette Score: {silhouette_avg:.3f}")
sc.tl.tsne(adata)
# 绘制 t-SNE 图
sc.pl.tsne(adata, color='cluster', palette='tab10', title="Clustering Stability")
# 保存图形为 DPI 为 600 的 TIFF 文件
plt.gcf().savefig('03.tif', format='tif', dpi=600)
# 显示图形
plt.show()