METHODS · PROGRAMMING

R / Python & Machine Learning

医学科研编程与算法工具箱

从数据对象、统计单位和研究问题出发，理解常用包如何被调用，以及模型输出应该怎样验证。

10 个常用包10 个经典算法No browser execution

R / PYTHON · METHODS TOOLBOX

代码不是黑箱：先确认对象，再读输出

这里不提供浏览器端执行环境。每段代码都标注输入对象、统计单位、关键参数、预期输出和不能证明的结论；复制到本地后，再用自己的数据和版本复现。

调用包确认对象验证泛化

01 · COMMON PACKAGES

医学科研常用包

按研究问题选择包，不按“哪个包最强”排序。

tidyverse

数据整理 · 如何把原始临床记录整理成可审计的分析表？

输入对象: 患者/样本/细胞长表
最小调用: library(tidyverse)
关键参数: group_by() 的统计单位summarise() 的分母join() 的键
如何读输出: 输出表每行代表一个明确统计单位，并保留缺失和重复记录。
统计单位: 患者、donor 或样本，取决于研究问题
常见误用: 把长表行数当成独立样本数。
不能单独证明: 数据整理完成不代表测量无偏或研究设计可识别。

R · REPRODUCE

R：按 donor 汇总临床表

把细胞级或就诊级记录还原到独立 donor/患者统计单位，避免把观测行误当作生物学重复。

下载 .R

INPUTclinical_cells.csvdonor_idconditionvalue

ENVIRONMENTR 4.5 · dplyr 1.1 · readr 2.1

1library(readr)
2library(dplyr)
3 
4records <- read_csv('clinical_cells.csv')
5donor_summary <- records %>%
6  group_by(donor_id, condition) %>%
7  summarise(
8    n_cells = n(),
9    median_value = median(value, na.rm = TRUE),
10    .groups = 'drop'
11  )

SOURCESEnhancing Reproducibility through Rigor and Transparency

ggplot2

科研绘图 · 如何把统计对象和视觉编码对应起来？

输入对象: tidy data frame
最小调用: library(ggplot2) ggplot(df, aes(x, y, colour = group)) + geom_point()
关键参数: x/y 的单位colour 的编码facet 的分母
如何读输出: 图例、坐标、点代表的统计单位和分母必须能回到数据表。
统计单位: 图中每个点/线/柱对应的观测单位
常见误用: 只调整颜色和主题，却不说明每个点代表什么。
不能单独证明: 图形美观或分组分离不等于统计显著或因果。

此包的完整可执行课例尚未收录。上方仅为 API 调用示意，需自行准备输入对象；不以其他包的代码代替。

DESeq2

Bulk RNA-seq · 如何在 raw counts 和设计矩阵上比较表达？

输入对象: gene × sample raw count matrix
最小调用: library(DESeq2)
关键参数: designsize factordispersionpadj
如何读输出: 同时读取 log2FC、CI/SE、P 值和 FDR，并回到样本级重复。
统计单位: 样本，而不是基因矩阵中的列数以外的伪重复
常见误用: 把 TPM 或细胞级行直接送进 DESeq2。
不能单独证明: 差异表达不能单独证明调控机制。

R · REPRODUCE

R：用 DESeq2 做 count-level 差异分析

在样本级 count 矩阵和设计矩阵上估计离散度与效应，不把 TPM 或细胞行直接当作独立重复。

下载 .R

INPUTgene × sample raw countssample metadatadesign formula

ENVIRONMENTR 4.5 · Bioconductor DESeq2

1library(DESeq2)
2dds <- DESeqDataSetFromMatrix(countData = counts,
3  colData = metadata, design = ~ batch + condition)
4dds <- DESeq(dds)
5res <- results(dds, contrast = c('condition', 'PF', 'Control'))

SOURCESDESeq2 Bioconductor package manual Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

Seurat

scRNA-seq · 如何组织细胞级表达、QC、降维和聚类？

输入对象: gene × cell count matrix + metadata
最小调用: library(Seurat) obj <- CreateSeuratObject(counts)
关键参数: assay/layernfeaturesdimsresolution
如何读输出: 每个分析层都要记录输入矩阵、参数和统计单位。
统计单位: 细胞用于表示，donor 用于疾病推断
常见误用: 只保存 UMAP 图片而不保存原始 counts、metadata 和参数。
不能单独证明: Seurat 对象完成不等于注释或机制成立。

R · REPRODUCE

R：计算 scRNA QC 指标

在逐样本分布中检查 nFeature、nCount 和线粒体比例，而不是直接套用一个全局 cutoff。

下载 .R

INPUTSeurat objectRNA assaysample_id

ENVIRONMENTR 4.5 · Seurat 5

1library(Seurat)
2obj[['percent.mt']] <- PercentageFeatureSet(obj, pattern = '^MT-')
3qc <- FetchData(obj, vars = c('nFeature_RNA', 'nCount_RNA', 'percent.mt', 'sample_id'))
4qc_by_sample <- aggregate(. ~ sample_id, data = qc, FUN = median)

SOURCESComprehensive Integration of Single-Cell Data Current best practices in single-cell RNA-seq analysis: a tutorial

DropletUtils

scRNA-seq · 如何从 raw droplets 识别细胞关联 barcode？

输入对象: raw barcode × gene count matrix
最小调用: library(DropletUtils) emptyDrops(counts(raw))
关键参数: ambient poollowerFDRraw versus filtered
如何读输出: 读取 barcode-level FDR 和背景模型，不把 knee 当作单独证据。
统计单位: barcode/droplet
常见误用: 从 filtered matrix 反推 emptyDrops。
不能单独证明: cell call 不能单独证明细胞健康或没有 doublet。

R · REPRODUCE

R：用 EmptyDrops 检查 barcode 是否偏离 ambient background

在 raw droplet 矩阵中计算背景模型和 FDR；不能从 filtered matrix 反推空滴概率。

下载 .R

INPUTraw_feature_bc_matrixraw barcode × gene counts

ENVIRONMENTR 4.5 · Bioconductor DropletUtils

1library(DropletUtils)
2raw <- read10xCounts('raw_feature_bc_matrix')
3res <- emptyDrops(counts(raw))
4called <- !is.na(res$FDR) & res$FDR <= 0.001
5summary(called)

SOURCESDropletUtils Bioconductor package manual Current best practices in single-cell RNA-seq analysis: a tutorial

Python

pandas

数据整理 · 如何在患者级整理 metadata、缺失和重复？

输入对象: DataFrame
最小调用: import pandas as pd
关键参数: merge keysgroupby unitmissingnessduplicate rows
如何读输出: 输出每行统计单位、缺失机制和 join 后行数变化。
统计单位: 患者、样本或 donor
常见误用: 不检查 merge 后行数膨胀。
不能单独证明: 表格整洁不代表变量定义正确。

Python · REPRODUCE

Python：患者级 Logistic 回归与校准

演示患者级二分类预测，区分模型排序能力、概率校准和临床效用。

下载 .py

INPUTpatient × feature tablebinary outcomepatient-level split

ENVIRONMENTPython 3.12 · scikit-learn 1.6

1from sklearn.model_selection import GroupKFold, cross_val_predict
2from sklearn.pipeline import make_pipeline
3from sklearn.preprocessing import StandardScaler
4from sklearn.linear_model import LogisticRegression
5from sklearn.metrics import roc_auc_score, brier_score_loss
6 
7model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000))
8cv = GroupKFold(n_splits=5)
9prob = cross_val_predict(model, X, y, groups=patient_id, cv=cv, method='predict_proba')[:, 1]
10print(roc_auc_score(y, prob), brier_score_loss(y, prob))

SOURCESTRIPOD+AI reporting guideline scikit-learn User Guide

Python

scanpy/anndata

scRNA-seq · 如何保留单细胞数据的矩阵层和邻居图？

输入对象: AnnData
最小调用: import scanpy as sc adata = sc.read_10x_mtx(path)
关键参数: obs/varlayersn_top_genesn_pcsn_neighbors
如何读输出: 区分 raw、normalized、PCA 和 neighbors 对象，不能只看 UMAP。
统计单位: 细胞用于表示，donor 用于比较
常见误用: 覆盖原始 counts 或把 UMAP 当作统计输入。
不能单独证明: AnnData 对象不自动证明批次已校正。

Python · REPRODUCE

Python：Scanpy 的最小 scRNA 表示流程

从已完成 cell calling 和样本级 QC 的计数对象继续分析；显式复制 counts 与 lognorm 层，再计算 HVG、PCA 和邻居图。不会自动进行批次校正或组间检验。

下载 .py

INPUTqc_counts.h5ad：X 是尚未归一化的非负整数计数，已完成 cell calling/QCobs.sample_id 已按唯一细胞 ID 对齐且无缺失var_names 为人类基因符号；其他物种或 Ensembl ID 必须改用对应线粒体基因注释

ENVIRONMENTPython + Scanpy/AnnData；需支持 pca(mask_var=...)，运行时记录实际包版本；未提供环境锁文件。

1import numpy as np
2import scanpy as sc
3from scipy import sparse
4 
5adata = sc.read_h5ad('qc_counts.h5ad')
6assert adata.obs_names.is_unique and adata.var_names.is_unique
7assert 'sample_id' in adata.obs and adata.obs['sample_id'].notna().all()
8values = adata.X.data if sparse.issparse(adata.X) else np.asarray(adata.X)
9assert np.isfinite(values).all() and (values >= 0).all()
10assert np.allclose(values, np.round(values)), 'X must contain counts, not log-expression'
11assert (np.asarray(adata.X.sum(axis=1)).ravel() > 0).all()
12adata.layers['counts'] = adata.X.copy()
13adata.var['mt'] = adata.var_names.str.startswith('MT-')
14assert adata.var['mt'].any(), 'Check species and gene identifier annotation'
15sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, inplace=True)
16sc.pp.normalize_total(adata, target_sum=1e4)
17sc.pp.log1p(adata)
18adata.layers['lognorm'] = adata.X.copy()
19sc.pp.highly_variable_genes(adata, flavor='seurat', n_top_genes=min(2000, adata.n_vars))
20n_pcs = min(30, adata.n_obs - 1, int(adata.var['highly_variable'].sum()) - 1)
21assert n_pcs >= 2, 'Too few cells or variable genes for this example'
22sc.pp.pca(adata, n_comps=n_pcs, mask_var='highly_variable', random_state=0)
23sc.pp.neighbors(adata, n_neighbors=min(15, adata.n_obs - 1), n_pcs=min(20, n_pcs), random_state=0)

关键行说明

L5这里从事先准备的 QC 后计数对象起步，不是从 raw droplets 起步；保留所有待分析基因和 sample_id。整数性检查不能代替检查数据来源。
L12copy() 建立独立的 counts 层；后续原位归一化 X 不会覆盖这份原始计数。
L14MT- 仅适用于本例约定的人类基因符号；找不到匹配基因时停止，不能把错误注释解释为线粒体比例为零。
L18lognorm 明确保存总量归一化再 log1p 的结果；它不是原始 counts，也不是差异分析的独立生物学重复。
L19seurat 风格 HVG 使用对数归一化输入；2000 是演示上限，不是每个数据集的最优基因数。
L22PCA 明确只使用 HVG，但不删除对象中其他基因；counts 和 lognorm 仍保留同一完整基因轴。
L23邻居图用 PCA 表示。维数和邻居数的 min 只防止小对象越界，不能代替参数敏感性分析。

EXPECTED OUTPUT

layers['counts'] 保存输入计数，layers['lognorm'] 与 X 保存对数归一化矩阵；var['highly_variable']、obsm['X_pca']、obsp['distances']/['connectivities'] 分别记录特征、坐标和图。本段未计算 UMAP、聚类或差异表达。

BEFORE YOU TRUST IT

确认 counts 与输入计数逐元素相同，且层形状和 obs/var 轴一致
按 sample_id 检查 QC 与 PCA，不能把成功创建 AnnData 当作质控已通过
比较候选 PC/邻居数对邻域和注释的影响；需要保存时再显式调用 write_h5ad

COMMON ERROR

只调用 normalize_total/log1p 就声称自动保存 counts 或 normalized layer；AnnData 层必须显式赋值。

SOURCESScanpy official API documentation SCANPY: large-scale single-cell gene expression data analysis Scanpy API: normalize_total Scanpy API: pca and mask_var AnnData API: AnnData.layers

Python

scikit-learn

机器学习 · 如何训练和验证经典预测模型？

输入对象: patient × feature matrix + outcome
最小调用: from sklearn.model_selection import GroupKFold
关键参数: group splitpipelinecross-validationmetrics
如何读输出: 分别报告 discrimination、calibration 和外部验证，不把 AUC 当成临床效用。
统计单位: 患者或其他独立研究单位
常见误用: 在全体数据上先选特征再交叉验证。
不能单独证明: 预测性能不能单独证明治疗机制或临床获益。

Python · REPRODUCE

Python：患者级 Logistic 回归与校准

演示患者级二分类预测，区分模型排序能力、概率校准和临床效用。

下载 .py

INPUTpatient × feature tablebinary outcomepatient-level split

ENVIRONMENTPython 3.12 · scikit-learn 1.6

1from sklearn.model_selection import GroupKFold, cross_val_predict
2from sklearn.pipeline import make_pipeline
3from sklearn.preprocessing import StandardScaler
4from sklearn.linear_model import LogisticRegression
5from sklearn.metrics import roc_auc_score, brier_score_loss
6 
7model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000))
8cv = GroupKFold(n_splits=5)
9prob = cross_val_predict(model, X, y, groups=patient_id, cv=cv, method='predict_proba')[:, 1]
10print(roc_auc_score(y, prob), brier_score_loss(y, prob))

SOURCESTRIPOD+AI reporting guideline scikit-learn User Guide

Python

statsmodels

统计推断 · 如何在 Python 中得到系数、CI 和模型诊断？

输入对象: 分析表 + 预先定义的设计矩阵
最小调用: import statsmodels.api as sm
关键参数: formula/designlinkrobust SEdiagnostics
如何读输出: 同时报告效应、CI、模型假设和统计单位。
统计单位: 由研究设计定义的独立单位
常见误用: 只取系数表中的 P 值。
不能单独证明: 模型拟合不等于因果识别。

此包的完整可执行课例尚未收录。上方仅为 API 调用示意，需自行准备输入对象；不以其他包的代码代替。

Python

lifelines

生存分析 · 如何把删失和时间零点纳入生存模型？

输入对象: patient × time/event/covariates
最小调用: from lifelines import CoxPHFitter
关键参数: duration_colevent_coltime zeroproportional hazards
如何读输出: 读取 HR、CI、事件数和假设诊断。
统计单位: 患者
常见误用: 只报告 HR 不说明删失和随访起点。
不能单独证明: Cox 关联不自动等于治疗因果效应。

Python · REPRODUCE

Python：生存数据的 Cox 模型

把事件时间、删失指示和协变量放在患者级 Cox 模型中，并检查比例风险假设。

下载 .py

INPUTpatient tabletimeeventcovariates

ENVIRONMENTPython 3.12 · lifelines

1from lifelines import CoxPHFitter
2 
3cph = CoxPHFitter()
4cph.fit(df, duration_col='time', event_col='event')
5cph.print_summary()
6cph.check_assumptions(df, p_value_threshold=0.05)

SOURCESTRIPOD+AI reporting guideline survival R package documentation Regression Models and Life-Tables

02 · CLASSICAL MACHINE LEARNING

先问医学问题，再选算法

算法输出的是可验证的统计对象，不是自动生成的机制结论。

supervised

线性回归

连续结局与协变量的平均差异是多少？

统计单位患者/样本

结局：连续数值

模型前提

线性关系或已建模的非线性
残差和方差诊断
独立单位定义正确

验证与指标

残差图
CI
外部或时间验证
RMSE
MAE
R²
95% CI

数据泄漏风险

在全体样本上先做缺失填补或特征选择。

系数关联不能单独证明因果。

该算法的可执行课例尚未收录；当前提供方法说明与验证要求，不展示其他算法的代码作为替代。

supervised

Logistic 回归

患者是否发生预定义二分类结局？

统计单位患者

结局：0/1 结局

模型前提

结局定义和时间窗固定
连续变量函数形式可诊断
事件数足以支撑参数

验证与指标

ROC-AUC
校准曲线
Brier score
外部验证
OR
95% CI
AUC
Brier

数据泄漏风险

把结局发生后的变量或同一患者重复记录放入预测特征。

高 AUC 不等于临床获益。

Python · REPRODUCE

Python：患者级 Logistic 回归与校准

演示患者级二分类预测，区分模型排序能力、概率校准和临床效用。

下载 .py

INPUTpatient × feature tablebinary outcomepatient-level split

ENVIRONMENTPython 3.12 · scikit-learn 1.6

1from sklearn.model_selection import GroupKFold, cross_val_predict
2from sklearn.pipeline import make_pipeline
3from sklearn.preprocessing import StandardScaler
4from sklearn.linear_model import LogisticRegression
5from sklearn.metrics import roc_auc_score, brier_score_loss
6 
7model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000))
8cv = GroupKFold(n_splits=5)
9prob = cross_val_predict(model, X, y, groups=patient_id, cv=cv, method='predict_proba')[:, 1]
10print(roc_auc_score(y, prob), brier_score_loss(y, prob))

SOURCESTRIPOD+AI reporting guideline scikit-learn User Guide

supervised

LASSO / Elastic Net

高维特征中哪些变量能形成可验证的预测模型？

统计单位患者 × 特征矩阵

结局：二分类、连续或生存结局

模型前提

特征处理只在训练折完成
lambda 通过交叉验证选择
保留选择不确定性

验证与指标

嵌套交叉验证
外部队列
校准
CV loss
AUC/RMSE
校准
系数路径

数据泄漏风险

用全体样本预筛特征或标准化。

被选中不等于生物学驱动因子。

R · REPRODUCE

R：用 glmnet 做正则化预测

用两个预先指定的数值变量演示二分类 LASSO：训练集内选择 lambda，验证集只用于预测。高维数据沿用相同输入约束，不能在全体样本上先筛特征。

下载 .R

INPUTtrain/test：已按患者划分，每人一行，patient_id 不重叠age、baseline_score：预测时已知的数值变量，无缺失train$outcome：0/1；本例 5 折要求两类各至少 5 人，仅为代码运行条件

ENVIRONMENTR + glmnet；运行时记录 R.version.string 与 packageVersion('glmnet')，此字段不是环境锁文件。

1library(glmnet)
2predictors <- c('age', 'baseline_score')
3stopifnot(all(c('patient_id', predictors, 'outcome') %in% names(train)),
4          all(c('patient_id', predictors) %in% names(test)))
5stopifnot(!anyNA(train$patient_id), !anyNA(test$patient_id))
6stopifnot(!anyDuplicated(train$patient_id), !anyDuplicated(test$patient_id))
7stopifnot(length(intersect(train$patient_id, test$patient_id)) == 0L)
8stopifnot(all(vapply(train[predictors], is.numeric, logical(1))),
9          all(vapply(test[predictors], is.numeric, logical(1))))
10x <- as.matrix(train[, predictors, drop = FALSE])
11x_test <- as.matrix(test[, predictors, drop = FALSE])
12y <- train$outcome
13stopifnot(all(is.finite(x)), all(is.finite(x_test)),
14          !anyNA(y), all(y %in% c(0, 1)), length(unique(y)) == 2L)
15stopifnot(min(table(y)) >= 5L, identical(colnames(x), colnames(x_test)))
16set.seed(2026)
17foldid <- integer(length(y))
18for (value in c(0, 1)) {
19  rows <- which(y == value)
20  foldid[rows] <- sample(rep(seq_len(5L), length.out = length(rows)))
21}
22fit <- cv.glmnet(x, y, family = 'binomial', alpha = 1,
23                 foldid = foldid, type.measure = 'deviance')
24prob <- drop(predict(fit, newx = x_test, s = 'lambda.1se', type = 'response'))

关键行说明

L2明确列出预测变量；outcome 和 patient_id 不进入训练或验证特征矩阵。示例不处理分类编码或缺失填补，这些步骤需在训练折内拟合。
L3先确认必需字段存在；R 对缺列的 $ 访问会返回 NULL，仅检查缺失值或重复值不能发现缺列。
L7患者不能跨训练/验证；这里要求每人一行。重复测量数据需要患者分组交叉验证，不应直接使用本例分折。
L11验证矩阵使用与训练完全相同的列及顺序；即使 test 中保存了 outcome，也不会将它作为预测变量。
L15每类至少 5 人只是本例构造 5 折的最低运行条件，不是模型开发样本量建议。
L20只用训练集的结局分层分折；5 折和种子均是演示设置，应按研究设计预先指定。
L22alpha=1 为 LASSO。glmnet 的标准化在各次拟合内完成；外部预筛、填补和调参仍需放入相应训练折。
L24lambda.1se 按训练交叉验证损失选择较强收缩方案；输出为 outcome=1 的预测概率，还未验证校准或临床效用。

EXPECTED OUTPUT

fit 含训练交叉验证偏差与 lambda；prob 按 test 行序保存 outcome=1 的概率。代码不计算外部验证指标，不能把训练 CV 损失称为外部验证结果。

BEFORE YOU TRUST IT

核对 patient_id 唯一且训练/验证无交集
确认预测变量测量时间早于预测时点，特征矩阵不含 outcome
核对每折两类人数；另用保留的验证结局评估区分度、校准与不确定性

COMMON ERROR

用 model.matrix(~ . - 1, test) 把 test$outcome 或患者 ID 也加入预测矩阵；列数相同也不保证列含义和顺序相同。

SOURCESTRIPOD+AI reporting guideline glmnet CRAN documentation glmnet reference manual: cv.glmnet and predict.cv.glmnet

supervised

Random Forest

非线性和交互关系能否提高患者级预测？

统计单位患者/样本

结局：分类或回归结局

模型前提

训练与验证按独立单位切分
特征测量时间点一致
缺失处理不泄漏

验证与指标

重复交叉验证
外部验证
校准和决策曲线
AUC/PR-AUC
RMSE
Brier
calibration

数据泄漏风险

同一患者多个时间点随机分到不同折。

变量重要性不是机制证据。

该算法的可执行课例尚未收录；当前提供方法说明与验证要求，不展示其他算法的代码作为替代。

supervised

Gradient Boosting

逐步拟合残差能否改善预测误差？

统计单位患者/样本

结局：分类或回归结局

模型前提

学习率、树深和迭代数需在训练数据内选择
验证集独立
输出概率需校准

验证与指标

嵌套 CV
外部验证
calibration
AUC/PR-AUC
Brier
MAE/RMSE
calibration slope

数据泄漏风险

用验证集反复挑选 boosting 轮数。

性能提升不能证明病因或治疗有效。

该算法的可执行课例尚未收录；当前提供方法说明与验证要求，不展示其他算法的代码作为替代。

supervised

Support Vector Machine

高维特征是否存在可泛化的分类边界？

统计单位患者/样本

结局：分类结局

模型前提

特征尺度在训练折内处理
核函数与 C/γ 通过验证选择
类别不平衡显式处理

验证与指标

分层或分组 CV
外部验证
概率校准
AUC
PR-AUC
sensitivity/specificity
calibration

数据泄漏风险

全体数据标准化后再切分。

分类边界不等于临床诊断机制。

该算法的可执行课例尚未收录；当前提供方法说明与验证要求，不展示其他算法的代码作为替代。

unsupervised

PCA

数据中最大的线性变化轴由什么驱动？

统计单位细胞/样本 × 特征矩阵

结局：低维表示和 loading

模型前提

输入层和缩放方式明确
高方差技术因素已检查
PC 数量通过下游稳定性解释

验证与指标

loading 检查
协变量相关
邻居稳定性
explained variance
loading
neighbor overlap

数据泄漏风险

用结局标签挑选 PCA 轴并称为无监督发现。

PC 方向不是因果方向。

Python · REPRODUCE

Python：Scanpy 的最小 scRNA 表示流程

从已完成 cell calling 和样本级 QC 的计数对象继续分析；显式复制 counts 与 lognorm 层，再计算 HVG、PCA 和邻居图。不会自动进行批次校正或组间检验。

下载 .py

ENVIRONMENTPython + Scanpy/AnnData；需支持 pca(mask_var=...)，运行时记录实际包版本；未提供环境锁文件。

1import numpy as np
2import scanpy as sc
3from scipy import sparse
4 
5adata = sc.read_h5ad('qc_counts.h5ad')
6assert adata.obs_names.is_unique and adata.var_names.is_unique
7assert 'sample_id' in adata.obs and adata.obs['sample_id'].notna().all()
8values = adata.X.data if sparse.issparse(adata.X) else np.asarray(adata.X)
9assert np.isfinite(values).all() and (values >= 0).all()
10assert np.allclose(values, np.round(values)), 'X must contain counts, not log-expression'
11assert (np.asarray(adata.X.sum(axis=1)).ravel() > 0).all()
12adata.layers['counts'] = adata.X.copy()
13adata.var['mt'] = adata.var_names.str.startswith('MT-')
14assert adata.var['mt'].any(), 'Check species and gene identifier annotation'
15sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, inplace=True)
16sc.pp.normalize_total(adata, target_sum=1e4)
17sc.pp.log1p(adata)
18adata.layers['lognorm'] = adata.X.copy()
19sc.pp.highly_variable_genes(adata, flavor='seurat', n_top_genes=min(2000, adata.n_vars))
20n_pcs = min(30, adata.n_obs - 1, int(adata.var['highly_variable'].sum()) - 1)
21assert n_pcs >= 2, 'Too few cells or variable genes for this example'
22sc.pp.pca(adata, n_comps=n_pcs, mask_var='highly_variable', random_state=0)
23sc.pp.neighbors(adata, n_neighbors=min(15, adata.n_obs - 1), n_pcs=min(20, n_pcs), random_state=0)

关键行说明

L5这里从事先准备的 QC 后计数对象起步，不是从 raw droplets 起步；保留所有待分析基因和 sample_id。整数性检查不能代替检查数据来源。
L12copy() 建立独立的 counts 层；后续原位归一化 X 不会覆盖这份原始计数。
L14MT- 仅适用于本例约定的人类基因符号；找不到匹配基因时停止，不能把错误注释解释为线粒体比例为零。
L18lognorm 明确保存总量归一化再 log1p 的结果；它不是原始 counts，也不是差异分析的独立生物学重复。
L19seurat 风格 HVG 使用对数归一化输入；2000 是演示上限，不是每个数据集的最优基因数。
L22PCA 明确只使用 HVG，但不删除对象中其他基因；counts 和 lognorm 仍保留同一完整基因轴。
L23邻居图用 PCA 表示。维数和邻居数的 min 只防止小对象越界，不能代替参数敏感性分析。

EXPECTED OUTPUT

BEFORE YOU TRUST IT

确认 counts 与输入计数逐元素相同，且层形状和 obs/var 轴一致
按 sample_id 检查 QC 与 PCA，不能把成功创建 AnnData 当作质控已通过
比较候选 PC/邻居数对邻域和注释的影响；需要保存时再显式调用 write_h5ad

COMMON ERROR

只调用 normalize_total/log1p 就声称自动保存 counts 或 normalized layer；AnnData 层必须显式赋值。

SOURCESScanpy official API documentation SCANPY: large-scale single-cell gene expression data analysis scikit-learn User Guide Scanpy API: normalize_total Scanpy API: pca and mask_var AnnData API: AnnData.layers

unsupervised

k-means

在给定 k 下，哪些观测靠近同一质心？

统计单位患者/细胞 × 数值特征

结局：cluster labels

模型前提

距离度量与尺度有科学解释
k 是预设或敏感性分析对象
标签需跨 seed/重采样检查

验证与指标

ARI/Jaccard
轮廓系数
donor 覆盖
marker/外部标签
inertia
silhouette
ARI
Jaccard

数据泄漏风险

先看结局再决定 k 或挑最漂亮的聚类。

聚类标签不是天然细胞类型或病因。

Python · REPRODUCE

Python：聚类与稳定性检查

演示聚类标签如何依赖特征空间、随机种子和参数，并用 ARI 比较重复结果。

下载 .py

INPUTpatient/cell × feature matrix预处理后的特征seed list

ENVIRONMENTPython 3.12 · scikit-learn 1.6

1from sklearn.cluster import KMeans
2from sklearn.metrics import adjusted_rand_score
3 
4labels_a = KMeans(n_clusters=4, n_init=20, random_state=1).fit_predict(X)
5labels_b = KMeans(n_clusters=4, n_init=20, random_state=2).fit_predict(X)
6ari = adjusted_rand_score(labels_a, labels_b)
7print(f'ARI={ari:.3f}')

SOURCESscikit-learn User Guide

survival

Cox 回归

协变量与事件发生率的相对关系如何？

统计单位患者

结局：生存时间 + 事件/删失

模型前提

time zero 明确
比例风险可检查
事件数和删失机制可描述

验证与指标

Schoenfeld/比例风险诊断
HR CI
外部验证
HR
95% CI
concordance
calibration

数据泄漏风险

用随访后信息预测基线结局。

观察性 HR 不能自动证明治疗因果。

Python · REPRODUCE

Python：生存数据的 Cox 模型

把事件时间、删失指示和协变量放在患者级 Cox 模型中，并检查比例风险假设。

下载 .py

INPUTpatient tabletimeeventcovariates

ENVIRONMENTPython 3.12 · lifelines

1from lifelines import CoxPHFitter
2 
3cph = CoxPHFitter()
4cph.fit(df, duration_col='time', event_col='event')
5cph.print_summary()
6cph.check_assumptions(df, p_value_threshold=0.05)

SOURCESTRIPOD+AI reporting guideline survival R package documentation Regression Models and Life-Tables

evaluation

模型评估与校准

模型概率是否能在新患者中保持可解释？

统计单位独立患者级预测

结局：预测概率与实际结局

模型前提

验证数据未参与模型选择
时间/机构外推范围明确
事件数可支持区间

验证与指标

外部验证
bootstrap CI
calibration-in-the-large
decision curve
AUC
PR-AUC
Brier
calibration slope/intercept

数据泄漏风险

用测试集挑模型、阈值或报告指标。

单一 AUC 不能证明临床获益、因果或可实施性。

Python · REPRODUCE

Python：患者级 Logistic 回归与校准

演示患者级二分类预测，区分模型排序能力、概率校准和临床效用。

下载 .py

INPUTpatient × feature tablebinary outcomepatient-level split

ENVIRONMENTPython 3.12 · scikit-learn 1.6

1from sklearn.model_selection import GroupKFold, cross_val_predict
2from sklearn.pipeline import make_pipeline
3from sklearn.preprocessing import StandardScaler
4from sklearn.linear_model import LogisticRegression
5from sklearn.metrics import roc_auc_score, brier_score_loss
6 
7model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000))
8cv = GroupKFold(n_splits=5)
9prob = cross_val_predict(model, X, y, groups=patient_id, cv=cv, method='predict_proba')[:, 1]
10print(roc_auc_score(y, prob), brier_score_loss(y, prob))

SOURCESTRIPOD+AI reporting guideline scikit-learn User Guide