质量控制与异常样本剔除：别让几颗老鼠屎坏了一锅粥

数据拿到手，第一件事是什么？跑模型？别急。我见过太多人，上来就把数据塞进流水线，结果跑出来的结果一塌糊涂。说白了，问题往往不在算法，而在数据本身。

这一章我们聊聊怎么给数据「体检」。质量控制，英文叫QC，核心就三件事：

• 找出离群的样本——那些明显不对劲的数据点
• 评估样本整体质量——用PCA看看有没有明显的结构问题
• 检查样本间一致性——相关系数矩阵能帮你发现配错或污染的样本

嗯，我们一个个来看。

4.1 为什么要先做质量控制？

你想想看，一个实验里难免有几个样本出了状况：加样时污染了、测序深度不够、或者干脆就是标签贴错了。如果不做质控，这些「问题样本」会带偏整个分析结果。

我个人习惯把质控分成两步走：

1. 技术质控——检查测序质量、比对率、基因检出数等基础指标
2. 统计质控——用PCA、离群点检测、相关系数矩阵找出统计上的异常

这一章我们重点讲第二步。技术质控通常是上游流程做的，但统计质控往往是数据科学家自己的责任。

4.2 PCA：一眼看穿数据的大致结构

PCA，主成分分析，说白了就是把高维数据投影到二维平面。你一眼就能看出样本有没有明显的集群结构、有没有明显偏离的个体。

4.2.1 什么时候用PCA做质控？

• 批次效应检测——如果样本按批次聚集，说明批次效应存在
• 离群样本检测——远离95%置信椭圆的样本值得怀疑
• 分组验证——处理组和对照组是否自然分离？

嗯，这里有个坑。我遇到过这样的情况：PCA一看，两个处理组完美分开。当时我还挺高兴，结果后来发现是批次效应搞的鬼——同一个批次的样本全跑一边去了。所以看PCA图时，记得用不同颜色标记批次和分组，两者重叠才是好信号。

4.2.2 实操代码片段

用Python做PCA质控，我通常这样写：

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

# 输入：log2归一化后的表达矩阵，行=基因，列=样本
expr = pd.read_csv("expression_log2norm.csv", index_col=0)

# 转置：PCA按样本做
expr_t = expr.T

# 标准化
scaler = StandardScaler()
expr_scaled = scaler.fit_transform(expr_t)

# PCA计算
pca = PCA(n_components=10)
pca_result = pca.fit_transform(expr_scaled)

# 整合为DataFrame
pca_df = pd.DataFrame(pca_result[:, :2], 
                       columns=['PC1', 'PC2'],
                       index=expr_t.index)

# 计算每个样本到中心的距离（马氏距离近似）
from scipy.spatial.distance import mahalanobis
cov = np.cov(expr_scaled.T)
inv_cov = np.linalg.inv(cov)
mean_vec = np.mean(expr_scaled, axis=0)

distances = []
for i in range(expr_scaled.shape[0]):
    d = mahalanobis(expr_scaled[i, :], mean_vec, inv_cov)
    distances.append(d)

# 卡方检验阈值：自由度=10，p=0.01
threshold = 23.209  # chi2.ppf(0.99, df=10)
outliers_pca = np.array(distances) > threshold
print(f"PCA检测到的离群样本数：{np.sum(outliers_pca)}")

为什么不直接画图看？因为样本多了以后，肉眼判断主观性太强。用马氏距离配合卡方阈值，可以量化地给出判断标准。

4.3 离群点检测：多重方法交叉验证

单靠PCA是不够的。我记得有一次，PCA图看着挺正常，但LOF（局部离群因子）检测揪出了好几个样本。后来查实验记录，发现这些样本的RNA提取时间确实和别的不一样。

常用的离群点检测方法有三种，各有优劣：

4.3.1 IQR法，我最常用的方法

为什么呢？因为它不需要假设数据分布。你想想看，基因表达数据有几个是正态分布的？

# IQR离群点检测 - 按样本的基因表达均值/中位数等统计量
import numpy as np

def detect_outliers_iqr(data, factor=1.5):
    """
    data: 一维数组，例如每个样本的表达量中位数
    """
    q1 = np.percentile(data, 25)
    q3 = np.percentile(data, 75)
    iqr = q3 - q1
    lower_bound = q1 - factor * iqr
    upper_bound = q3 + factor * iqr
    outliers = (data < lower_bound) | (data > upper_bound)
    return outliers

# 示例：每个样本的基因数(检测到的基因数)
gene_counts = expr.apply(lambda col: np.sum(col > 0))  # 每列非零基因数
outliers_iqr = detect_outliers_iqr(gene_counts.values)
print(f"IQR检测到的离群样本数：{np.sum(outliers_iqr)}")

4.3.2 LOF的实战调参经验

LOF有一个关键参数：k，也就是邻居数。我踩过的坑是：k设得太小（比如k=5），会误判正常样本为离群；设得太大（k=100），又会漏掉真正的离群点。

我个人习惯：样本数小于100时，k取min(20, n_samples-1)；样本数大于100时，k取20~30之间。

from sklearn.neighbors import LocalOutlierFactor

# LOF检测
clf = LocalOutlierFactor(n_neighbors=20, contamination=0.05)
y_pred = clf.fit_predict(expr_scaled)  # -1 = 离群, 1 = 正常

outliers_lof = y_pred == -1
print(f"LOF检测到的离群样本数：{np.sum(outliers_lof)}")

4.4 样本相关性矩阵：发现配对错误或污染的利器

这一点容易被忽略。我曾经在一个项目中，发现两个样本的相关系数异常低。一开始以为是批次效应，后来查记录，发现是标注配错了——一个是肿瘤组织，一个是正常组织，标签颠倒了。

样本相关性矩阵能告诉你什么？

• 组内相关性高（>0.9）：同一处理组内样本应该高度相似
• 组间相关性低（<0.8）：不同组间应该有一定差异
• 全部样本相关性接近1：可能是有严重的批次效应，或者数据没有生物学差异
• 某个样本与所有样本的相关性都特别低（<0.7）：疑似污染或技术失败

4.4.1 用热图直观判断

import seaborn as sns
import matplotlib.pyplot as plt

# 计算样本间Pearson相关系数
corr_matrix = expr.corr(method='pearson')

# 热图可视化
plt.figure(figsize=(12, 10))
sns.heatmap(corr_matrix, 
            cmap='RdBu_r', 
            vmin=0.5, vmax=1.0,
            xticklabels=True, yticklabels=True,
            square=True)
plt.title("样本相关性矩阵热图")
plt.tight_layout()
plt.savefig("sample_correlation_heatmap.png", dpi=150)
plt.close()

# 找出异常样本：与组内其他样本的平均相关性低于阈值的
group_labels = pd.read_csv("metadata.csv")["group"]  # 分组信息
threshold_inner = 0.85  # 组内最低平均相关性

anomaly_samples = []
for group in group_labels.unique():
    group_samples = group_labels[group_labels == group].index
    if len(group_samples) < 2:
        continue
    sub_corr = corr_matrix.loc[group_samples, group_samples]
    mean_corr = sub_corr.mean(axis=1)
    bad = mean_corr[mean_corr < threshold_inner]
    anomaly_samples.extend(bad.index.tolist())

4.5 综合决策：如何确定最终剔除名单？

我通常的做法是：如果一个样本被至少两个方法标记为异常，就列入剔除名单。为什么不相信单一方法？你想想看，Z-score可能对厚尾分布不敏感，PCA可能被少量离群点拉偏，LOF的参数选择又主观。三重保险，心里踏实。

# 综合三种方法的标记结果
decision_matrix = pd.DataFrame({
    'sample': expr.columns,
    'pca_outlier': outliers_pca,
    'iqr_outlier': outliers_iqr,
    'lof_outlier': outliers_lof
})

# 计算每个样本被标记为异常的次数
decision_matrix['n_votes'] = decision_matrix.iloc[:, 1:].sum(axis=1)

# 剔除：至少2票异常
samples_to_remove = decision_matrix[decision_matrix['n_votes'] >= 2]['sample'].tolist()

print(f"建议剔除的样本数：{len(samples_to_remove)}")
print("剔除样本列表：", samples_to_remove)

# 保存清洗后的表达矩阵
expr_clean = expr.drop(columns=samples_to_remove)
expr_clean.to_csv("expression_clean.csv")

4.6 质控报告：记录是专业性的体现

做了这么多检测，结果要记录下来。我个人习惯生成一份简洁的质控报告，包含以下内容：

• 样本总数：原始样本数、剔除后样本数
• 各方法检测到的离群数：PCA、IQR、LOF、相关性矩阵各标记了多少
• 最终剔除名单及被标记的原因
• 质控前后的数据概况对比：基因数、中位数表达量等
• 保留的PCA和相关性热图

# 简单的质控报告生成
report = f"""
========================================
   质量控制报告 - {pd.Timestamp.now().strftime('%Y-%m-%d')}
========================================

原始样本数: {expr.shape[1]}
清洗后样本数: {expr_clean.shape[1]}
剔除样本数: {len(samples_to_remove)}

离群点检测结果:
- PCA (马氏距离): {np.sum(outliers_pca)} 个
- IQR (因子=1.5): {np.sum(outliers_iqr)} 个
- LOF (k=20): {np.sum(outliers_lof)} 个

剔除样本列表:
{samples_to_remove}

被标记原因:
{pd.DataFrame(decision_matrix[decision_matrix['n_votes']>=2])}

质控前后数据概况:
- 原始数据: 基因数={expr.shape[0]}, 样本数={expr.shape[1]}
- 清洗后: 基因数={expr_clean.shape[0]}, 样本数={expr_clean.shape[1]}
"""

with open("qc_report.txt", "w") as f:
    f.write(report)

print("质控报告已保存至 qc_report.txt")

好，这一套流程走下来，你的数据应该干净多了。质控这件事，说白了就是一个「宁可错杀一千，不可放过一个」的活——但也不能杀太多，否则统计功效会下降。嗯，多练几次就会有手感，知道你的数据在什么情况下算是「干净」了。