{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.11.11","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":99575,"databundleVersionId":11943919,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"先安装依赖：`anndata` 和 `scanpy` .","metadata":{}},{"cell_type":"code","source":"!pip install anndata scanpy --quiet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T08:32:47.125257Z","iopub.execute_input":"2025-05-06T08:32:47.125558Z","iopub.status.idle":"2025-05-06T08:32:55.168205Z","shell.execute_reply.started":"2025-05-06T08:32:47.125535Z","shell.execute_reply":"2025-05-06T08:32:55.166938Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# scDNAm Cluster Competition ","metadata":{}},{"cell_type":"markdown","source":"## 读取数据","metadata":{}},{"cell_type":"code","source":"import anndata as ad\nimport numpy as np\nimport warnings\nwarnings.filterwarnings('ignore')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:08:42.742699Z","iopub.execute_input":"2025-05-06T09:08:42.743446Z","iopub.status.idle":"2025-05-06T09:08:42.750206Z","shell.execute_reply.started":"2025-05-06T09:08:42.743411Z","shell.execute_reply":"2025-05-06T09:08:42.748943Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def read_data(file: str) -> tuple[np.ndarray]:\n    # 读取数据并返回生物数据和批次\n    data = ad.read_h5ad(file)\n    batch = data.obs[\"batch\"]\n    print(f\"阅读数据完毕，发现共有 {np.unique(batch)} 个批次！\")\n    \n    \n    \n    d = {\"batch1\": 0, \"batch2\": 1}\n    batch = np.array(list(map(d.get, batch)))\n    print(\"批次数据处理完毕！\")\n\n    return data.X.astype(np.float16), batch\n\ndata, batch = read_data(\"/kaggle/input/data-mining-hw-2/final_dataset.h5ad/final_dataset.h5ad\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T08:46:49.060884Z","iopub.execute_input":"2025-05-06T08:46:49.061999Z","iopub.status.idle":"2025-05-06T08:48:55.768199Z","shell.execute_reply.started":"2025-05-06T08:46:49.061931Z","shell.execute_reply":"2025-05-06T08:48:55.767109Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"batch","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T08:51:08.453376Z","iopub.execute_input":"2025-05-06T08:51:08.453715Z","iopub.status.idle":"2025-05-06T08:51:08.4627Z","shell.execute_reply.started":"2025-05-06T08:51:08.453681Z","shell.execute_reply":"2025-05-06T08:51:08.461585Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 数据预处理","metadata":{}},{"cell_type":"markdown","source":"### 查看缺失比例和缺失数量","metadata":{}},{"cell_type":"code","source":"def cal_nan_num(X: np.ndarray):\n    row_nan_counts = np.isnan(X).sum(axis=1)\n    \n    # 统计每一列的缺失值数量\n    col_nan_counts = np.isnan(X).sum(axis=0)\n    \n    # 对行缺失值数量进行排序并输出排行\n    row_ranking = np.argsort(row_nan_counts)[::-1]\n    print(\"行缺失值数量排行（索引）:\", row_ranking)\n    print(\"行缺失值数量:\", row_nan_counts[row_ranking])\n    \n    # 对列缺失值数量进行排序并输出排行\n    col_ranking = np.argsort(col_nan_counts)[::-1]\n    print(\"列缺失值数量排行（索引）:\", col_ranking)\n    print(\"列缺失值数量:\", col_nan_counts[col_ranking])\n\n    no_nan_col_count = np.sum(col_nan_counts == 0)\n    print(\"没有缺失值的列的数量:\", no_nan_col_count)\n\ncal_nan_num(data)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:02:16.823699Z","iopub.execute_input":"2025-05-06T09:02:16.824888Z","iopub.status.idle":"2025-05-06T09:02:38.387882Z","shell.execute_reply.started":"2025-05-06T09:02:16.824851Z","shell.execute_reply":"2025-05-06T09:02:38.386851Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"缺失值数量非常庞大，可以看到我们每一行都有上万的缺失值，某些行甚至缺失 200000 列数据，如何填补缺失值，如何挑选列对数据分析非常关键.","metadata":{}},{"cell_type":"code","source":"# 先丢弃缺失值比例过高的列\nratio = np.isnan(data).sum(axis=0) / data.shape[0]\nkeep = ratio <= 0.95     # 如果缺失值超过 95% 就没有必要保留了.\ndata = data[:, keep]\ndata.shape","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:13:15.76704Z","iopub.execute_input":"2025-05-06T09:13:15.767384Z","iopub.status.idle":"2025-05-06T09:13:30.856288Z","shell.execute_reply.started":"2025-05-06T09:13:15.767361Z","shell.execute_reply":"2025-05-06T09:13:30.855394Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"然后我们在这里使用均值填充，这里使用截断均值.","metadata":{}},{"cell_type":"code","source":"# 进度条\nfrom rich.progress import Progress","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:21:11.203069Z","iopub.execute_input":"2025-05-06T09:21:11.203472Z","iopub.status.idle":"2025-05-06T09:21:11.511404Z","shell.execute_reply.started":"2025-05-06T09:21:11.203445Z","shell.execute_reply":"2025-05-06T09:21:11.51042Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def impute_missing_with_truncated_mean(data, block_size = 5000):\n    num_cols = data.shape[1]\n    \n    with Progress() as progress:\n        task = progress.add_task(\"[green]Imputing columns...\", total=num_cols)\n        for start in range(0, num_cols, block_size):\n            end = min(start + block_size, num_cols)\n            block = data[:, start:end]\n            for col_idx in range(block.shape[1]):\n                col = block[:, col_idx]\n                non_nan_col = col[~np.isnan(col)]\n                if non_nan_col.size > 0:\n                    lower_bound = np.percentile(non_nan_col, 5)\n                    upper_bound = np.percentile(non_nan_col, 95)\n                    truncated_col = non_nan_col[(non_nan_col >= lower_bound) & (non_nan_col <= upper_bound)]\n                    mean_val = np.mean(truncated_col)\n                    block[np.isnan(col), col_idx] = mean_val\n            data[:, start:end] = block\n            progress.update(task, advance=end - start)\n    return data","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:21:50.275696Z","iopub.execute_input":"2025-05-06T09:21:50.276734Z","iopub.status.idle":"2025-05-06T09:21:50.284944Z","shell.execute_reply.started":"2025-05-06T09:21:50.276689Z","shell.execute_reply":"2025-05-06T09:21:50.283905Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data = impute_missing_with_truncated_mean(data)\ndata","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:22:08.877423Z","iopub.execute_input":"2025-05-06T09:22:08.877745Z","iopub.status.idle":"2025-05-06T09:25:21.955368Z","shell.execute_reply.started":"2025-05-06T09:22:08.87772Z","shell.execute_reply":"2025-05-06T09:25:21.954039Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 选取特征","metadata":{}},{"cell_type":"markdown","source":"特征的选取非常重要，并且还要考虑到批次效应，我们这里考虑分批次筛选.","metadata":{}},{"cell_type":"code","source":"from scipy.stats import median_abs_deviation","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:31:51.495699Z","iopub.execute_input":"2025-05-06T09:31:51.496041Z","iopub.status.idle":"2025-05-06T09:31:51.915291Z","shell.execute_reply.started":"2025-05-06T09:31:51.49602Z","shell.execute_reply":"2025-05-06T09:31:51.914127Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def select_highly_variable_features(data, batch_labels=None, method='variance', top_n=1000):\n    if batch_labels is None:\n        # 无批次效应，直接计算整个数据集的统计量\n        if method == 'variance':\n            stats = np.var(data, axis=0)\n        elif method == 'mad':\n            stats = median_abs_deviation(data, axis=0)\n        else:\n            raise ValueError(\"方法必须是 'variance' 或 'mad'。\")\n    else:\n        unique_batches = np.unique(batch_labels)\n        stats_per_batch = []\n        for batch in unique_batches:\n            batch_data = data[batch_labels == batch]\n            if method == 'variance':\n                batch_stats = np.var(batch_data, axis=0)\n            elif method == 'mad':\n                batch_stats = median_abs_deviation(batch_data, axis=0)\n            else:\n                raise ValueError(\"方法必须是 'variance' 或 'mad'。\")\n            stats_per_batch.append(batch_stats)\n        stats = np.mean(stats_per_batch, axis=0)\n\n    # 选取 top_n 个特征\n    top_indices = np.argsort(stats)[-top_n:]\n    return top_indices","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:31:18.130596Z","iopub.execute_input":"2025-05-06T09:31:18.130934Z","iopub.status.idle":"2025-05-06T09:31:18.138624Z","shell.execute_reply.started":"2025-05-06T09:31:18.13091Z","shell.execute_reply":"2025-05-06T09:31:18.137343Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"selected = select_highly_variable_features(data, batch_labels=batch,\n                                      method = 'mad', top_n = 30000)\ndata = data[:, selected]\ndata","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:31:54.593346Z","iopub.execute_input":"2025-05-06T09:31:54.594186Z","iopub.status.idle":"2025-05-06T09:34:28.431316Z","shell.execute_reply.started":"2025-05-06T09:31:54.59416Z","shell.execute_reply":"2025-05-06T09:34:28.430032Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 消除批次效应","metadata":{}},{"cell_type":"markdown","source":"批次效应是本次作业最难处理的地方，我们这里使用已经有的批次处理方法.","metadata":{}},{"cell_type":"code","source":"!pip install combat --quiet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:35:12.639053Z","iopub.execute_input":"2025-05-06T09:35:12.639471Z","iopub.status.idle":"2025-05-06T09:35:20.031415Z","shell.execute_reply.started":"2025-05-06T09:35:12.639446Z","shell.execute_reply":"2025-05-06T09:35:20.030116Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"from combat.pycombat import pycombat\nimport pandas as pd","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:36:24.304436Z","iopub.execute_input":"2025-05-06T09:36:24.305794Z","iopub.status.idle":"2025-05-06T09:36:24.311143Z","shell.execute_reply.started":"2025-05-06T09:36:24.305764Z","shell.execute_reply":"2025-05-06T09:36:24.309893Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data_df = pd.DataFrame(data)\ndata_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:36:36.604215Z","iopub.execute_input":"2025-05-06T09:36:36.60454Z","iopub.status.idle":"2025-05-06T09:36:36.652375Z","shell.execute_reply.started":"2025-05-06T09:36:36.604519Z","shell.execute_reply":"2025-05-06T09:36:36.651386Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"data = pycombat(data_df.T,batch).values.T","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:37:07.349893Z","iopub.execute_input":"2025-05-06T09:37:07.350242Z","iopub.status.idle":"2025-05-06T09:37:16.103815Z","shell.execute_reply.started":"2025-05-06T09:37:07.350221Z","shell.execute_reply":"2025-05-06T09:37:16.102923Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## PCA 降维","metadata":{}},{"cell_type":"code","source":"from sklearn.decomposition import PCA","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:37:27.651696Z","iopub.execute_input":"2025-05-06T09:37:27.652068Z","iopub.status.idle":"2025-05-06T09:37:28.103025Z","shell.execute_reply.started":"2025-05-06T09:37:27.652044Z","shell.execute_reply":"2025-05-06T09:37:28.101933Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pca = PCA(n_components = 40)\ndata = pca.fit_transform(data)\ndata","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:38:00.958635Z","iopub.execute_input":"2025-05-06T09:38:00.959026Z","iopub.status.idle":"2025-05-06T09:38:04.884666Z","shell.execute_reply.started":"2025-05-06T09:38:00.958993Z","shell.execute_reply":"2025-05-06T09:38:04.883498Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"最后绘图查看我们的数据.","metadata":{}},{"cell_type":"code","source":"!pip install umap-learn --quiet","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:39:59.368684Z","iopub.execute_input":"2025-05-06T09:39:59.369122Z","iopub.status.idle":"2025-05-06T09:40:03.590604Z","shell.execute_reply.started":"2025-05-06T09:39:59.369096Z","shell.execute_reply":"2025-05-06T09:40:03.58942Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nimport umap","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:40:03.592704Z","iopub.execute_input":"2025-05-06T09:40:03.593053Z","iopub.status.idle":"2025-05-06T09:40:46.331929Z","shell.execute_reply.started":"2025-05-06T09:40:03.59302Z","shell.execute_reply":"2025-05-06T09:40:46.3311Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def plot_umap_clustering(data, labels, title=\"UMAP Clustering Plot\"):\n    # 设置绘图风格为 ggplot\n    plt.style.use('ggplot')\n\n    # 使用 UMAP 进行降维\n    reducer = umap.UMAP(random_state=42)\n    embedding = reducer.fit_transform(data)\n\n    # 绘制聚类图\n    plt.figure(figsize=(10, 8))\n    scatter = plt.scatter(embedding[:, 0], embedding[:, 1], c=labels, cmap='Set1')\n    plt.colorbar(scatter)\n    plt.title(title)\n    plt.xlabel('UMAP 1')\n    plt.ylabel('UMAP 2')\n    plt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:46:39.870781Z","iopub.execute_input":"2025-05-06T09:46:39.87117Z","iopub.status.idle":"2025-05-06T09:46:39.87874Z","shell.execute_reply.started":"2025-05-06T09:46:39.871144Z","shell.execute_reply":"2025-05-06T09:46:39.877544Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"plot_umap_clustering(data, batch)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:46:41.958121Z","iopub.execute_input":"2025-05-06T09:46:41.95847Z","iopub.status.idle":"2025-05-06T09:46:53.366688Z","shell.execute_reply.started":"2025-05-06T09:46:41.958445Z","shell.execute_reply":"2025-05-06T09:46:53.365314Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## 高斯混合聚类","metadata":{}},{"cell_type":"code","source":"from sklearn.mixture import GaussianMixture","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"gm = GaussianMixture(n_components=10, random_state=42).fit(data)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"labels = gm.predict(data)","metadata":{"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"pd.DataFrame({'ID': range(5052),\n                      'TARGET': labels}).to_csv(f'submission.csv',index=False)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-05-06T09:52:50.983208Z","iopub.execute_input":"2025-05-06T09:52:50.983606Z","iopub.status.idle":"2025-05-06T09:52:51.004523Z","shell.execute_reply.started":"2025-05-06T09:52:50.983578Z","shell.execute_reply":"2025-05-06T09:52:51.003443Z"}},"outputs":[],"execution_count":null}]}