异常检测算法¶

import numpy as np
import pandas as pd
from sklearn.ensemble import IsolationForest
from sklearn.neighbors import LocalOutlierFactor
from sklearn.preprocessing import StandardScaler

# ===== 1. 简单 Z-score 异常检测（单细胞 QC 例子）=====
def z_score_filter(values, threshold=3.0):
    """
    基于 Z-score 检测异常值
    threshold: 超过几个标准差算异常（通常3）
    返回: 是否为正常样本的布尔数组
    """
    mean = np.mean(values)   # 均值
    std = np.std(values)     # 标准差
    z_scores = np.abs((values - mean) / std)  # 每个值的 Z-score
    return z_scores <= threshold  # True=正常，False=异常

# 模拟 scRNA-seq QC 指标
np.random.seed(42)
n_cells = 500
n_genes_per_cell = np.random.normal(2000, 500, n_cells)  # 每个细胞检测到的基因数
n_genes_per_cell = np.clip(n_genes_per_cell, 0, None)    # 非负

# 插入几个异常细胞（空液滴或双联体）
n_genes_per_cell[490:495] = [20, 15, 10, 8000, 9000]    # 低基因=空液滴，高基因=双联体

is_normal = z_score_filter(n_genes_per_cell, threshold=3)
print(f'正常细胞: {is_normal.sum()} / {len(is_normal)}')
print(f'异常细胞索引: {np.where(~is_normal)[0]}')

# ===== 2. Isolation Forest 高维异常检测 =====
# 构造多维 QC 矩阵（实际用 scanpy 的 obs 数据）
qc_features = np.column_stack([
    n_genes_per_cell,                                      # 检测到的基因数
    np.random.normal(0.1, 0.05, n_cells).clip(0, 1),      # 线粒体基因比例
    np.random.normal(5000, 1000, n_cells).clip(0, None)   # 总 UMI 计数
])
# 手动设置异常值
qc_features[490:495, 1] = [0.9, 0.85, 0.8, 0.02, 0.01]  # 高线粒体=死细胞

# 归一化
scaler = StandardScaler()
qc_scaled = scaler.fit_transform(qc_features)

# 训练 Isolation Forest
iso_forest = IsolationForest(
    n_estimators=100,     # 树的数量
    contamination=0.05,   # 预期异常比例（5%）
    random_state=42
)
iso_pred = iso_forest.fit_predict(qc_scaled)  # -1=异常，1=正常
anomaly_idx = np.where(iso_pred == -1)[0]
print(f'\nIsolation Forest 检测到 {len(anomaly_idx)} 个异常细胞')
print(f'异常细胞索引（前10个）: {anomaly_idx[:10]}')

# ===== 3. 自编码器异常检测 =====
import torch
import torch.nn as nn

class AnomalyAutoencoder(nn.Module):
    """
    自编码器：正常数据重建误差小，异常数据重建误差大
    用于高维基因表达矩阵的异常样本检测
    """
    def __init__(self, input_dim, hidden_dim=32, bottleneck_dim=8):
        super().__init__()
        # 编码器：压缩到低维
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),    # 输入 -> 隐藏
            nn.ReLU(),
            nn.Linear(hidden_dim, bottleneck_dim)  # 隐藏 -> 瓶颈（最低维）
        )
        # 解码器：从低维重建
        self.decoder = nn.Sequential(
            nn.Linear(bottleneck_dim, hidden_dim),  # 瓶颈 -> 隐藏
            nn.ReLU(),
            nn.Linear(hidden_dim, input_dim)         # 隐藏 -> 重建输出
        )

    def forward(self, x):
        z = self.encoder(x)   # 压缩
        x_recon = self.decoder(z)  # 重建
        return x_recon

# 模拟基因表达数据（50个基因）
gene_dim = 50
normal_data = torch.randn(400, gene_dim)    # 正常样本
anomaly_data = torch.randn(10, gene_dim) * 5  # 异常样本（尺度不同）
all_data = torch.cat([normal_data, anomaly_data], dim=0)

# 只用正常数据训练（无监督）
ae_model = AnomalyAutoencoder(input_dim=gene_dim)
optimizer = torch.optim.Adam(ae_model.parameters(), lr=1e-3)
loss_fn = nn.MSELoss()

for epoch in range(100):
    ae_model.train()
    recon = ae_model(normal_data)        # 重建正常数据
    loss = loss_fn(recon, normal_data)   # 重建误差
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()

# 计算所有样本的重建误差（异常样本重建误差应更大）
ae_model.eval()
with torch.no_grad():
    recon_all = ae_model(all_data)
    recon_errors = ((recon_all - all_data) ** 2).mean(dim=1)  # 每个样本的重建误差

# 设定阈值（训练集的 95 百分位）
with torch.no_grad():
    normal_errors = ae_model(normal_data)
    normal_errors = ((normal_errors - normal_data) ** 2).mean(dim=1)
    threshold = torch.quantile(normal_errors, 0.95)  # 第 95 百分位数为阈值

predicted_anomalies = (recon_errors > threshold).sum().item()
print(f'\n自编码器检测到 {predicted_anomalies} 个异常样本（真实异常: 10）')