单细胞 RNA-seq 质量控制：工作原理

P(X = k) = (λ^k * e^(-λ)) / k!

import numpy as np

# 示例数据（包含异常值）
data = np.array([1000, 1200, 1100, 1300, 1050, 50000])  # 50000 是异常值

# 标准差
mean = np.mean(data)
std = np.std(data)
upper_threshold_std = mean + 3 * std

print(f"标准差方法：")
print(f"  均值：{mean:.0f}")
print(f"  标准差：{std:.0f}")
print(f"  上限：{upper_threshold_std:.0f}")
print(f"  被标记为异常：{np.sum(data > upper_threshold_std)} 个")

# MAD
median = np.median(data)
mad = np.median(np.abs(data - median))
upper_threshold_mad = median + 3 * mad

print(f"\nMAD 方法：")
print(f"  中位数：{median:.0f}")
print(f"  MAD: {mad:.0f}")
print(f"  上限：{upper_threshold_mad:.0f}")
print(f"  被标记为异常：{np.sum(data > upper_threshold_mad)} 个")

标准差方法:
  均值: 10775  ← 被异常值拉高
  标准差: 19517
  上限: 69326
  被标记为异常: 0 个  ← 未能检测到异常值！

MAD 方法:
  中位数: 1150  ← 不受异常值影响
  MAD: 75
  上限: 1375
  被标记为异常: 1 个  ← 成功检测！

MAD = median(|xi - median(x)|)

σ ≈ 1.4826 × MAD

# 检查过滤后的分布
adata_filtered = adata[filter_mask]

# 计算偏度和峰度
from scipy.stats import skew, kurtosis
skewness = skew(adata_filtered.obs['n_genes_by_counts'])
kurt = kurtosis(adata_filtered.obs['n_genes_by_counts'])

print(f"偏度：{skewness:.2f}")  # 期望接近 0
print(f"峰度：{kurt:.2f}")      # 期望接近 0（正态分布）

# 检查高 mt% 细胞是否表达细胞类型标志物
high_mt = adata.obs['pct_counts_mt'] > 20

# 心肌细胞标志物
if 'MYH6' in adata.var_names:
    cardiomyocyte = adata[:, 'MYH6'].X > 0
    overlap = (high_mt & cardiomyocyte).sum()
    print(f"高 mt% 且表达 MYH6 的细胞：{overlap}")
    # 如果重叠高，说明这些是正常心肌细胞，不应过滤

# 计算残差
from scipy.stats import linregress

slope, intercept, r_value, p_value, std_err = linregress(
    adata.obs['n_counts'],
    adata.obs['n_genes_by_counts']
)

# 预测值
predicted_genes = intercept + slope * adata.obs['n_counts']
residuals = adata.obs['n_genes_by_counts'] - predicted_genes

# 标记异常（残差 > 3 MAD）
mad_residuals = np.median(np.abs(residuals - np.median(residuals)))
outliers = np.abs(residuals) > 3 * mad_residuals

print(f"异常细胞数：{outliers.sum()}")

# 核糖体比例 vs 细胞周期
sc.tl.score_genes_cell_cycle(adata, ...)

# 可视化
import seaborn as sns
sns.boxplot(
    data=adata.obs,
    x='phase',
    y='pct_counts_ribo'
)
# 期望：S 期和 G2M 期核糖体比例高

理论双细胞率 ≈ (加载细胞数 / 回收细胞数) × 1000

例如：
  - 目标回收: 5000 细胞
  - 实际加载: 6000 细胞
  - 双细胞率: (6000 / 5000) × 0.8% ≈ 1.0%

import scrublet as scr

# 1. 创建 Scrublet 对象
scrub = scr.Scrublet(counts_matrix)

# 2. 检测双细胞
doublet_scores, predicted_doublets = scrub.scrub_doublets()

# 3. 双细胞得分分布
# - 低得分：单细胞
# - 高得分：可能的双细胞
# - 双峰分布：单细胞和双细胞分离良好

# 4. 自动阈值选择
# Scrublet 自动在双峰之间的谷处设定阈值

# 查看双细胞得分分布
scrub.plot_doublet_score_histogram()

# 如果阈值不合理，手动调整
threshold = 0.25  # 根据直方图选择
predicted_doublets = doublet_scores > threshold

# 检查数据分布
for metric in ['n_genes_by_counts', 'n_counts', 'pct_counts_mt']:
    median = np.median(adata.obs[metric])
    mad = np.median(np.abs(adata.obs[metric] - median))

    print(f"\n{metric}:")
    print(f"  中位数：{median:.1f}")
    print(f"  MAD: {mad:.1f}")
    print(f"  3 MAD 阈值：[{median-3*mad:.1f}, {median+3*mad:.1f}]")

# 尝试不同 n_mads
for n_mads in [2, 3, 5]:
    filter_mask = (
        mad_filter(adata.obs['n_genes_by_counts'], n_mads, 'lower') &
        mad_filter(adata.obs['pct_counts_mt'], n_mads, 'upper')
    )
    retention = filter_mask.sum() / len(filter_mask) * 100
    print(f"n_mads={n_mads}: 保留率 {retention:.1f}%")

# 检查已知细胞群是否保留
marker_genes = {
    'T_cell': 'CD3D',
    'B_cell': 'MS4A1',
    'Monocyte': 'LYZ'
}

for cell_type, marker in marker_genes.items():
    if marker in adata.var_names:
        before = (adata[:, marker].X > 0).sum()
        after = (adata_qc[:, marker].X > 0).sum()
        retention = after / before * 100
        print(f"{cell_type} 保留率：{retention:.1f}%")

组织类型？
├─ 免疫细胞（PBMC）
│  └─ 标准 MAD (3)
│     └─ mt% 阈值: 中位数 + 3MAD (≈ 10-12%)
│
├─ 高能量组织（心脏、肌肉、肾）
│  └─ 提高 mt% 阈值 (5 MAD)
│     └─ 验证: 检查组织特异标志物
│
├─ 肿瘤
│  └─ 极宽松 (5 MAD)
│     └─ 不过滤低 UMI（可能是 CTC）
│
├─ 胚胎/干细胞
│  └─ 极宽松 (5 MAD)
│     └─ 允许低 n_genes
│
└─ 非模式生物
   └─ 验证线粒体基因前缀

# 模拟：严格过滤 vs 宽松过滤
strict_filter = mad_filter(adata.obs['pct_counts_mt'], 2, 'upper')  # 严格
permissive_filter = mad_filter(adata.obs['pct_counts_mt'], 5, 'upper')  # 宽松

# 分别聚类
adata_strict = adata[strict_filter].copy()
adata_permissive = adata[permissive_filter].copy()

sc.pp.pca(adata_strict)
sc.pp.neighbors(adata_strict)
sc.tl.leiden(adata_strict, resolution=0.5)

sc.pp.pca(adata_permissive)
sc.pp.neighbors(adata_permissive)
sc.tl.leiden(adata_permissive, resolution=0.5)

print(f"严格过滤聚类数：{adata_strict.obs['leiden'].nunique()}")
print(f"宽松过滤聚类数：{adata_permissive.obs['leiden'].nunique()}")
# 预期：宽松过滤保留更多稀有簇

# 比较：过滤前 vs 过滤后的 DE 结果
# （假设 sc.tl.rank_genes_groups 已运行）

# 高 mt% 细胞的伪基因
high_mt_genes = adata[adata.obs['pct_counts_mt'] > 20].var_names
print(f"高 mt% 细胞富集的基因：{high_mt_genes[:10]}")

# 这些"差异表达"可能是技术伪影

1. 初始 QC（宽松）
   ↓
2. 初步聚类和细胞类型注释
   ↓
3. 检查可疑簇
   - 高 mt%？
   - 低 n_genes？
   - 表达应激反应基因？
   ↓
4. 调整 QC 参数（如果需要）
   ↓
5. 重新分析