PCA和PLS-DA
我们使用SRBCT数据集来说明PCA和sPLS-DA。
安装并加载mixOmics包
BiocManager::install('mixOmics')
library(mixOmics)
示例数据集
示例数据是mixOmics包自带的经过标准化处理过的可以直接使用的数据，来自小圆蓝细胞肿瘤(SRBCT)。数据集包括63个样本的2308个基因的表达水平。样本分为四类:Burkitt淋巴瘤(BL) 8例，尤文氏肉瘤(EWS) 23例，成神经细胞瘤(NB) 12例，横纹肌肉瘤(RMS) 20例。
srbct数据集包含以下内容:
$gene:一个63行2308列的数据框。63个样本中2308个基因的表达水平。
$class：包含每个个体的类肿瘤的分类(共4个类)。
$gene.name：一个包含2308行和2列的数据框，包含关于基因的进一步信息。
data(srbct)
X = srbct$gene #the gene expression data
dim(X)
X[1:6,1:10]
> dim(X)
[1] 63 2308
> X[1:6,1:10]
g1 g2 g3 g4 g5 g6 g7 g8 g9 g10
EWS.T1 3.2025 0.0681 1.0460 0.1243 0.4941 3.1207 3.7106 1.8416 1.2607 2.9001
EWS.T2 1.6547 0.0710 1.0409 0.0520 0.2045 2.1609 2.4452 1.1473 0.7371 1.9989
EWS.T3 3.2779 0.1160 0.8926 0.1014 0.2818 1.9773 3.2590 1.4106 0.9548 2.0775
EWS.T4 1.0060 0.1906 0.4302 0.1035 0.2984 1.6804 5.8901 0.2958 0.7381 1.6610
EWS.T6 2.7098 0.2367 0.3693 0.2190 0.3711 1.7800 3.2376 0.6769 0.8546 0.6808
EWS.T7 2.0588 0.0823 0.9021 0.1288 0.3961 1.7199 2.1729 1.5609 0.6581 1.5038
summary(srbct$class)
> summary(srbct$class)
EWS BL NB RMS
23 8 12 20
主成分分析(PCA)
对基因表达数据进行初步的PCA分析，可以首次探索数据变异的主要来源。PCA是一种无监督分析，没有提供关于肿瘤类别的信息。为了理解所解释的变化量，我们将主成分数目(ncomp = 10)设置为一个相当大的数字。在PCA中，centering(中心化)使所有基因具有相同的零平均值，利于关注样本之间的差异。Scaling(缩放)的目的是在分析中给所有基因相似的权重，因为高方差的基因会被认为在PCA中有影响，但不一定有生物学相关性。
pca.srbct = pca(X, ncomp = 10, center = TRUE, scale = TRUE)
# 输出每个主成分可解释方差(变异)
plot(pca.srbct)
上面的柱状图显示了两个成分足以解释来自数据的大部分或信息)。
在下面的样本图中，样本用前两个主成分表示，并根据肿瘤类型着色。这里我们观察到变异的主要来源可能不能用肿瘤类型来解释。注意，由于PCA是无监督的，为了可视化目的，我们只考虑PCA之后的样本类型信息。
plotIndiv(pca.srbct, group = srbct$class, ind.names = FALSE,
legend = TRUE, title = 'PCA on SRBCT')
PLS-DA 分析
对于判别分析，我们设置因子Y来表示每个样本的类别隶属度。在PLS-DA过程中，将Y因子转化为一个虚拟矩阵。
Y = srbct$class
summary(Y) #结果分类
> summary(Y)
EWS BL NB RMS
23 8 12 20
PLS-DA模型采用10个成分来评估最终模型所需的性能和成分数量(见下文)。
样品图:将样本投影到前两个成分所形成的子空间中。
srbct.plsda <- plsda(X, Y, ncomp = 10)
plotIndiv(srbct.plsda , comp = 1:2,
group = srbct$class, ind.names = FALSE,
ellipse = TRUE, legend = TRUE, title = 'PLSDA on SRBCT')
从样本图中可以看出，与无监督的PCA样本图相比，四种肿瘤类型明显分离。绘制每个类的置信椭圆以突出区分的强度(置信水平默认设置为95%，参数ellipse.level)。
在覆盖样本图之前，可以通过计算背景面来可视化预测区域。
# with background
background = background.predict(srbct.plsda, comp.predicted=2, dist = "max.dist")
#optional: xlim = c(-40,40), ylim = c(-30,30))
plotIndiv(srbct.plsda, comp = 1:2,
group = srbct$class, ind.names = FALSE, title = "Maximum distance",
legend = TRUE, background = background)
性能
PLS-DA模型的分类性能通过重复10次的5折交叉验证来评估。重复次数对于确保分类错误率的良好估计是必要的(因为cv -fold是以随机方式确定的)。从性能结果中我们可以决定选择最终的PLS模型的成分数量。
set.seed(2543) # 设置种子，便于重复
perf.plsda.srbct <- perf(srbct.plsda, validation = "Mfold", folds = 5,
progressBar = FALSE, auc = TRUE, nrepeat = 10)
# perf.plsda.srbct$error.rate 错误率
plot(perf.plsda.srbct, col = color.mixo(5:7), sd = TRUE, legend.position = "horizontal")
从性能图中，可以看出总体错误率和平衡错误率(BER)相似，从1个成分急剧下降到3个成分。6个成分后错误率趋于稳定。ncomp = 6时，BER和最大距离足以实现良好的性能(0.06错误率)。
auc.plsda = auroc(srbct.plsda, roc.comp = 6)
