{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# scDNAm Cluster Competition\n\n本 Notebook 解决 **scDNAm 聚类竞赛任务（数据挖掘作业二）**：对单细胞 DNA 甲基化数据进行无监督聚类，并在过程中消除批次效应。\n\n**核心思路**（基于 scanpy 单细胞分析流程）：\n1. 使用均值填充缺失值；\n2. `normalize_total` + `log1p` 进行文库归一化与对数转换；\n3. `scale` 标准化后做 PCA 降维；\n4. 使用 **Harmony** 校正批次效应，并基于校正后的表示构建 KNN 邻居图；\n5. UMAP 降维可视化；\n6. 通过 **KMeans + 轮廓系数** 评估合适的聚类数，并用 **Leiden** 算法进行最终聚类；\n7. 输出 Leiden 聚类结果用于提交。","metadata":{}},{"cell_type":"markdown","source":"## 环境准备\n\n运行核心流程需要 `harmonypy`（Harmony 批次校正）与 `leidenalg`（Leiden 聚类）。","metadata":{}},{"cell_type":"code","source":"# 安装依赖：Harmony（批次校正）与 leidenalg（Leiden 聚类）\n!pip install harmonypy leidenalg --quiet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:40:45.87917Z","iopub.execute_input":"2026-08-19T03:40:45.87939Z","iopub.status.idle":"2026-08-19T03:44:08.729117Z","shell.execute_reply.started":"2026-08-19T03:40:45.879365Z","shell.execute_reply":"2026-08-19T03:44:08.728041Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 第三方库导入","metadata":{}},{"cell_type":"code","source":"# 单细胞数据分析\nimport scanpy as sc\n\n# 科学计算与表格\nimport numpy as np\nimport pandas as pd\nimport scipy.sparse as sp\n\n# 机器学习：聚类与评估\nfrom sklearn.cluster import KMeans\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import (adjusted_rand_score, adjusted_mutual_info_score,\n                             normalized_mutual_info_score, homogeneity_score,\n                             silhouette_score)\n\n# 绘图\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n# 警告过滤\nimport warnings\nwarnings.filterwarnings('ignore')\n\n# 固定随机种子，保证结果可复现\nnp.random.seed(42)\n\nprint(\"[INFO] 模块导入完成！\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.731553Z","iopub.execute_input":"2026-08-19T03:44:08.73192Z","iopub.status.idle":"2026-08-19T03:44:08.747825Z","shell.execute_reply.started":"2026-08-19T03:44:08.731877Z","shell.execute_reply":"2026-08-19T03:44:08.745805Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 读取数据\n\n读取 `.h5ad` 格式的甲基化数据，并将数据类型降为 `float32` 以节省内存。数据附带 `batch` 列，记录样本所属批次。","metadata":{}},{"cell_type":"code","source":"# 读取数据\nadata = sc.read_h5ad(r\"D:\\下载\\data-mining-hw-2\\final_dataset.h5ad\\final_dataset.h5ad\")\n# 降低精度以节省内存\nadata.X = adata.X.astype(np.float32)\n\nprint(f\"[INFO] 数据形状: 样本数 {adata.shape[0]}, 特征数 {adata.shape[1]}\")\nprint(f\"[INFO] 批次信息: {np.unique(adata.obs['batch'])}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.748638Z","iopub.status.idle":"2026-08-19T03:44:08.749025Z","shell.execute_reply.started":"2026-08-19T03:44:08.748847Z","shell.execute_reply":"2026-08-19T03:44:08.748873Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 缺失值处理\n\n原始数据存在大量缺失值，这里先查看整体缺失比例，再用**列均值**进行填充。填充后构建新的 `AnnData` 对象用于后续预处理流程。","metadata":{}},{"cell_type":"code","source":"# 查看整体缺失值比例（兼容稀疏矩阵）\nX_arr = adata.X.toarray() if sp.issparse(adata.X) else adata.X\nnan_ratio = np.isnan(X_arr).mean()\nprint(f\"[INFO] 整体缺失值比例: {nan_ratio:.4f}\")\n\n# 使用均值填充缺失值\nimputer = SimpleImputer(strategy='mean')\nX_imputed = imputer.fit_transform(adata.X)\n\n# 构建新的 AnnData 对象，保留原始 obs（含 batch 信息）\nadata_pp = sc.AnnData(X=X_imputed, obs=adata.obs)\nprint(f\"[INFO] 缺失值填充完成，数据形状: {adata_pp.shape}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.750547Z","iopub.status.idle":"2026-08-19T03:44:08.751198Z","shell.execute_reply.started":"2026-08-19T03:44:08.750969Z","shell.execute_reply":"2026-08-19T03:44:08.750998Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 归一化与对数转换\n\n- `normalize_total`：文库归一化，使每个样本的总计数一致（target_sum=1e4），消除测序深度差异；\n- `log1p`：对数转换，使数据分布更接近正态，压缩长尾。","metadata":{}},{"cell_type":"code","source":"# 文库归一化\nsc.pp.normalize_total(adata_pp, target_sum=1e4)\n# 对数转换\nsc.pp.log1p(adata_pp)\nprint(\"[INFO] 归一化与对数转换完成\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.752861Z","iopub.status.idle":"2026-08-19T03:44:08.753214Z","shell.execute_reply.started":"2026-08-19T03:44:08.753073Z","shell.execute_reply":"2026-08-19T03:44:08.753094Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 标准化与 PCA 降维\n\n先对特征进行标准化（clip 到 `max_value=10` 避免极端值主导），再做 PCA 降维至 30 个主成分。","metadata":{}},{"cell_type":"code","source":"# 标准化，max_value 截断极端值\nsc.pp.scale(adata_pp, max_value=10)\n# PCA 降维\nsc.tl.pca(adata_pp, n_comps=30)\nprint(\"[INFO] PCA 降维完成，主成分数: 30\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.754735Z","iopub.status.idle":"2026-08-19T03:44:08.755085Z","shell.execute_reply.started":"2026-08-19T03:44:08.75495Z","shell.execute_reply":"2026-08-19T03:44:08.754966Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 查看各主成分的方差解释比例\nsc.pl.pca_variance_ratio(adata_pp, n_pcs=30, log=True)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.756425Z","iopub.status.idle":"2026-08-19T03:44:08.756722Z","shell.execute_reply.started":"2026-08-19T03:44:08.756592Z","shell.execute_reply":"2026-08-19T03:44:08.756608Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 批次效应校正（Harmony）\n\n批次效应是聚类任务中的关键难点。使用 **Harmony** 在 PCA 空间上校正批次效应，得到 `X_pca_harmony` 表示，并基于该表示构建 KNN 邻居图供后续聚类使用。","metadata":{}},{"cell_type":"code","source":"# 使用 Harmony 校正批次效应，结果存于 obsm['X_pca_harmony']\nsc.external.pp.harmony_integrate(adata_pp, 'batch')\n\n# 基于校正后的表示构建 KNN 邻居图\nsc.pp.neighbors(adata_pp, use_rep='X_pca_harmony')\nprint(\"[INFO] Harmony 批次校正与邻居图构建完成\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.759211Z","iopub.status.idle":"2026-08-19T03:44:08.75975Z","shell.execute_reply.started":"2026-08-19T03:44:08.759574Z","shell.execute_reply":"2026-08-19T03:44:08.759601Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## UMAP 降维可视化\n\n在邻居图基础上计算 UMAP 嵌入，用于二维可视化，直观检查批次效应是否已被消除以及潜在聚类结构。","metadata":{}},{"cell_type":"code","source":"# 计算 UMAP 嵌入\nsc.tl.umap(adata_pp)\n# 可视化 UMAP，按 batch 着色检查批次混合情况\nsc.pl.umap(adata_pp, color='batch')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.761607Z","iopub.status.idle":"2026-08-19T03:44:08.762041Z","shell.execute_reply.started":"2026-08-19T03:44:08.761839Z","shell.execute_reply":"2026-08-19T03:44:08.761863Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# 查看 PCA 空间中的样本分布\nsc.pl.pca(adata_pp)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.764118Z","iopub.status.idle":"2026-08-19T03:44:08.764551Z","shell.execute_reply.started":"2026-08-19T03:44:08.764339Z","shell.execute_reply":"2026-08-19T03:44:08.764365Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 聚类数选择（KMeans + 轮廓系数）\n\n在 Harmony 校正后的表示上，用 KMeans 尝试不同聚类数（2~20），并以**轮廓系数（Silhouette Score）**作为评估指标，选取最优聚类数 `best_k`。\n\n> 该步骤仅用于评估合适的聚类数，最终提交结果采用 Leiden 算法。","metadata":{}},{"cell_type":"code","source":"# 在 Harmony 校正后的表示上评估不同聚类数\nX = adata_pp.obsm['X_pca_harmony']\nmin_clusters = 2\nmax_clusters = 20\nscores = []\nbest_k = min_clusters\nbest_score = -np.inf\n\nfor k in range(min_clusters, max_clusters + 1):\n    kmeans = KMeans(n_clusters=k, random_state=42)\n    cluster_labels = kmeans.fit_predict(X)\n    if len(np.unique(cluster_labels)) > 1:  # 确保有多个聚类\n        score = silhouette_score(X, cluster_labels)\n        scores.append(score)\n        if score > best_score:\n            best_score = score\n            best_k = k\n\nprint(f\"最佳聚类数: {best_k}，Silhouette 分数: {best_score:.4f}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.765576Z","iopub.status.idle":"2026-08-19T03:44:08.765963Z","shell.execute_reply.started":"2026-08-19T03:44:08.765736Z","shell.execute_reply":"2026-08-19T03:44:08.765758Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## Leiden 聚类\n\nLeiden 算法是基于模块度的社区检测方法，能够稳健地发现图结构中的社区，是单细胞聚类的主流方法。这里在邻居图上以 `resolution=1.0` 进行聚类，结果存入 `obs['leiden']`。","metadata":{}},{"cell_type":"code","source":"# Leiden 社区检测\nsc.tl.leiden(adata_pp, resolution=1.0, key_added='leiden')\n# UMAP 上按 Leiden 聚类着色\nsc.pl.umap(adata_pp, color=[\"leiden\"])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.767593Z","iopub.status.idle":"2026-08-19T03:44:08.767926Z","shell.execute_reply.started":"2026-08-19T03:44:08.767739Z","shell.execute_reply":"2026-08-19T03:44:08.767755Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## KMeans 聚类对比\n\n使用上一步得到的最优 `best_k` 进行 KMeans 聚类，作为与 Leiden 结果的对比参考。","metadata":{}},{"cell_type":"code","source":"# 使用最佳 k 的 KMeans 聚类结果作为对比\nkmeans = KMeans(n_clusters=best_k, random_state=42)\nadata_pp.obs['kmeans'] = kmeans.fit_predict(adata_pp.obsm['X_pca_harmony']).astype(str)\nsc.pl.pca(adata_pp, color='kmeans')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.769289Z","iopub.status.idle":"2026-08-19T03:44:08.76969Z","shell.execute_reply.started":"2026-08-19T03:44:08.769496Z","shell.execute_reply":"2026-08-19T03:44:08.769517Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 结果输出\n\n将 **Leiden 聚类结果**整理为 `ID, TARGET` 格式并保存为 `result3.csv` 用于提交。","metadata":{}},{"cell_type":"code","source":"# 使用 Leiden 聚类结果作为最终提交\nsubmission = pd.DataFrame({\n    'ID': range(5052),\n    'TARGET': adata_pp.obs['leiden'].values\n})\nsubmission.to_csv(\"./result3.csv\", index=False)\n\nprint(f\"[INFO] 结果已保存至 result3.csv，共 {len(submission)} 条\")\nsubmission.head()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-08-19T03:44:08.771541Z","iopub.status.idle":"2026-08-19T03:44:08.77186Z","shell.execute_reply.started":"2026-08-19T03:44:08.771689Z","shell.execute_reply":"2026-08-19T03:44:08.771704Z"}},"outputs":[],"execution_count":null}]}