如何利用MATLAB中PCA得分矩阵定位RNAseq异常样本索引并过滤数据

先明确矩阵对应关系 #

你的normCounts是24427个基因 × 136个样本，转置后传入pca()，因此输出的score矩阵是136行 × 136列：

• 每一行对应1个原始样本（顺序和normCounts的列完全一致）
• 每一列对应一个主成分（按方差解释率从高到低排序）

步骤1：定位异常样本索引 #

方法1：用Hotelling T²统计量（最直接）

pca()输出的tsquared是每个样本的T²值，衡量样本到所有主成分中心的距离，超过统计阈值的即为异常：

% 设置显著性水平，常用0.05
alpha = 0.05;
% 计算T²统计量的阈值（基于F分布）
num_pcs = size(score, 2);
num_samples = size(score, 1);
threshold = finv(1 - alpha, num_pcs, num_samples - num_pcs);
% 找出异常样本的索引
outlier_idx = find(tsquared > threshold);

方法2：结合前几个主成分的得分（直观易验证）

通常前2-3个主成分解释了大部分方差，先可视化确认异常点：

% 绘制PC1 vs PC2散点图，标记样本索引
scatter(score(:,1), score(:,2));
text(score(:,1), score(:,2), string(1:136), 'VerticalAlignment', 'bottom');

手动观察离群点后，用标准差法量化筛选（可扩展到多个主成分）：

% 计算PC1得分的Z分数，用3σ原则筛选异常
pc1_z = (score(:,1) - mean(score(:,1))) / std(score(:,1));
outlier_idx = find(abs(pc1_z) > 3);

% 若结合PC1+PC2+PC3，计算马氏距离筛选
selected_scores = score(:,1:3);
md = mahal(selected_scores, selected_scores);
md_threshold = finv(1 - alpha, 3, 136-3);
outlier_idx = find(md > md_threshold);

方法3：聚类识别离群点

用DBSCAN聚类标记远离主簇的样本：

% eps和MinPts需根据你的数据调整
[idx, ~] = dbscan(score(:,1:2), 2, 3);
outlier_idx = find(idx == 0); % idx=0代表噪声点（异常样本）

步骤2：过滤原始归一化计数矩阵 #

得到outlier_idx后，剔除原始矩阵中对应列的样本（normCounts的列对应样本）：

% 保留非异常样本的计数矩阵
filtered_normCounts = normCounts(:, ~ismember(1:size(normCounts,2), outlier_idx));
% 验证过滤后的样本数
size(filtered_normCounts) % 应为24427 × (136 - 异常样本数)

内容的提问来源于stack exchange，提问作者GrimPillBilly