数据降维
在进行数据挖掘或者机器学习时，我们面临的数据往往是高维数据。
相较于低维数据，高维数据为我们提供了更多的信息和细节，也更好的描述了样本；但同时，很多高效且准确的分析方法也将无法使用。
处理高维数据和高维数据可视化是数据科学家们必不可少的技能。解决这个问题的方法便是降低数据的维度。在数据降维时，要使用尽量少的维度来表达较多原数据的特性和结构。
降维是通过减少数据中的指标（或变量）以化简数据的过程。这里的减少指标，并不是随意加减，而是用复杂的数理知识，得到几个“综合指标”来代表整个数据。
PCA原理
主成分分析（Principal Component Analysis, PCA）是一种线性降维算法，也是一种常用的数据预处理（Pre-Processing）方法。
它的目标是是用方差（Variance）来衡量数据的差异性，并将差异性较大的高维数据投影到低维空间中进行表示。
绝大多数情况下，我们希望获得两个主成分因子：分别是从数据差异性最大和次大的方向提取出来的，称为PC1(Principal Component 1) 和 PC2（Principal Component 2）。
PCA原理简述
PCA具体是如何实现的？PC1和PC2又是如何计算的呢？
举个简单的例子，Scores.xlsx 包含了约70名学生的全科考试成绩。其中每名学生是一个独立的样本，每门学科的成绩都是一个数据维度（共有13门成绩）。我们的目的是通过分析学生的考试成绩来判断学生的类别（理科、文科生，和体育、艺术特长生）。
下表为随机选取的六名学生的语文和数学考试成绩：
Step 1 散点图
制作为散点图：
上图中每个点代表了一个学生，X轴代表语文成绩，Y轴代表数学成绩。
Step 2 计算均值
然后分别取所有样本的X平均值和Y平均值，并将这两个值变为X、Y坐标，在图中画出这个点（用五角星表示）:
Step 3 中心化
按照图中箭头所示方向，将整个坐标系平移，使原点与五角星重叠。这样就获得了一个新的平面直角坐标系：
尽管此时坐标系和每个点的值都发生了变化，点与点之间的相对位置仍保持一致。
这个过程称为“中心化”。“中心化”处理的原因是，这些数字后续会参与统计运算，比如求样本方差。
Step 4 找到最优拟合线
找到这些点的最优拟合线（Line of Best Fit），也就找到了PC1，再通过原点做PC1的垂线，就找到了PC2：
处理三维数组时便会产生第三个因子（PC3），以此类推，数据的维度越大，因子的数量也就越多。当维度大于等于4的时候，我们是无法想象出图像的，但PC4确实存在；假设有x个维度，便可以做x-1条垂线，就能得到PCx。接下来要做的便是选取最能代表数据差异性的两个因子，作为PC1和PC2。
Step 5 投影(Projection)
按照下图所示，将点A投影到PC1上（六角星的位置），并计算其与原点之间的距离称为d1：
其余的五个点也做同样操作，得出d2至d5，再求这六个距离的平方和，称为PC1的特征值（Eigenvalue）。
Step 6 计算差异值
然后将PC1的特征值除以总样本数量减一（n-1），就计算出了PC1的差异值（Variation）：
Step 7 碎石图(scree plot)
以此类推，并选择差异值最大的两个因子作为PC1 和 PC2。假设在某个三维数组中，获得了PC1、PC2和PC3的差异值分别为18，7，5。通过计算（18+7）/ (18+7+5) ≈ 83.3% 得到结论：PC1 和 PC2 代表了这个三维数组83.3%的差异性。在本次分析的13个因子中，PC1和PC2描述了整组数据约81%的差异性：
Step 8 PCA图(PCA plot)
最后，再通过选中的PC1和PC2将样本映射回本身所在的坐标，就可以得到降维后的图像（PCA Plot）。
上图是在经过PCA降维处理后，得到的浮点图；图中每个点都代表了一名学生，其中文科生为绿色，理科生为红色，体育生和艺术生分别为黄色和蓝色。坐标系的X轴为PC1，Y轴为PC2。从图中我们可以明显地看出同类别学生的聚类（Clustering）趋势，也证明了PCA在降维的同时，尽可能地保留了原数据的特性。
PCA的理论推导
上述原理只是PCA原理的简单化说明，实际上，主成分不是原来的指标中的任何一个，而是由所有原有指标数据线性组合而来。
PCA是将一组变量通过正交变换转变成另一组变量的分析方法，以此来实现数据降维的目的。转换后得到的这一组变量，即是我们所说的主成分。也就是把多指标的数据用少数几个综合指标（主成分）替代，还原数据最本质特征。然后根据这些主成分对样本进行聚类，代表样本的点在坐标轴上距离越远，说明样本差异越大。同时，还可以看两个或多个分组之间能否分得开，能否找到差别，也是一种质控的手段。
这些主成分会依据方差的大小进行排序，称作主成分（PC）1、主成分2、……主成分r。而每个主成分的方差在这一组变量中的总方差中所占的比例，即是该主成分的可解释差异贡献度。一般我们要求前K个主成分的累计可解释差异贡献度大于80%，k小于4，因为最多作出三维图像，即仅考察贡献度前2或者前3的主成分，经过可视化后，即得到了二维或三维PCA散点图。在二维图像中，把横纵坐标百分比加起来，这就是它们可解释差异贡献度，这个数字越大模型的效果越好，三维图像是X+Y+Z，这个数字要求大于80%。
以最简单的模型为例：如何把二维数据压缩至一维。
数据压缩的目标是：让新数据的方差尽可能地大。这个标准能使得新数据尽可能地不丢失原有数据的信息，因为方差越大，数据间的差异越大。
如下图所示：有六个点，每个点有两个特征，分别对应x轴和y轴。需要把他们压缩成一维的数据，即每个点只有一个特征。因此要寻一条直线，让所有点投影到该直线上。该直线上的刻度即为新数据的值。
首先进行中心化处理。中心化的好处在于，我们寻求的直线必定经过原点。如下图所示，只需从所有经过原点的直线中，找一条直线，使得各个数据的方差最大。
如下图所示，设P是样本的一个数据。注意到这样一个性质：由于OP的距离是固定的，因此过原点作任何一条直线，记Q是P在该直线上的投影点，都有 OQ2 + PQ2 = 定值。
显然知道 OQ （带正负号的长度）是数据 P(xp, yp) 压缩后的值。
那么如下图所示：对于 A1 到 A6 六个点，分别作直线的垂线，交直线于 B1 到 B6 六个点。直线可以看做数轴，坐标系的原点就是数轴的原点。记 B1到 B6 六个点的刻度分别为 b1,…,b6。那么 b1,…,b6 就是降维后的数据。
也就说明 b 的均值为0。那么 b 的方差为：
我们的目标是让压缩后的数据方差最大，即让
最大。记点
到直线的距离为
，那么有 
= 定值。于是
= 定值。使
最大，那么等价于
最小。
这个结果很美妙，在这里和一元线性回归作比较。如下图所示，一共有九个点。我们要进行一元线性回归，所采取的策略是让平方损失最小，通俗地来说，就是让九条绿线的平方和最小。
如果是采用PCA，如下图所示，我们采取的策略是让九条红线的平方和最小。
事实上，一元线性回归也能采取“使点到直线距离的平方和最小”的策略。两种策略各有各的优势与缺陷。只不过“平方误差最小原则”已经被广泛地接受。
类似地，对于n维数据，投影到更小的k维空间也能同样地进行类比。
为什么PCA可以代替所有数据？
主成分为什么拽到可以代替所有数据？换句话说，为什么较多的变量经过数据变换之后变成了较少的几个变量呢？这不是会有大量的信息丢失掉吗？
这是因为部分变量其实是相互关联的，这就会造成数据冗余。而降维就可以帮助我们去除这些指标中重叠、多余的信息，把数据最本质和关键的信息提取出来。在我们最开始得到的一组变量中，变量之间并不是完全相互独立的。例如一个位点发生了变异，那么与之连锁的几个位点也大概率会发生变异；或者一个基因的表达量发生了变化，同一条通路中的其他基因的表达量也大概率会发生变化，即变量之间是存在相关性的。极端一点，假设两个位点完全连锁，那么我们去掉其中一个突变的所有信息，并不会影响总的信息含量。主成分分析也即是基于这样一种思想开展的，将变量之间根据相关性进行分解、合并和降维，类似于从n维空间到r维空间的投影。
PCA优缺点
优点：
1、以方差衡量信息的无监督学习，不受样本标签限制。
2、由于协方差矩阵对称，因此k个特征向量之间两两正交，也就是各主成分之间正交，正交就肯定线性不相关，可消除原始数据成分间的相互影响
3. 可减少指标选择的工作量
4.用少数指标代替多数指标，利用PCA降维是最常用的算法
5. 计算方法简单，易于在计算机上实现。
缺点：
1、主成分解释其含义往往具有一定的模糊性，不如原始样本完整
2、贡献率小的主成分往往可能含有对样本差异的重要信息，也就是可能对于区分样本的类别（标签）更有用
3、特征值矩阵的正交向量空间是否唯一有待讨论
4、无监督学习
PCA结果解读
在分析过程中，PCA可以让我们非常直观地看出各个样本之间的相似性。例如在一张PCA散点图中，数个样本的点聚在一起，那么就说明这几个样本之间的相似性非常高；反之，如果几个样本的点非常分散，则说明这几个样本之间的相似性比较低。
例1：
下图：几个组的样本对应的散点在组内呈现相互聚集的情况，说明组内的重复性比较好，样本数据非常相似，而组间则有较好的区分度。有的时候为了说明组内样本的相似程度，还会用一个椭圆将同一组的样本对应的散点全部囊括起来。
例2：
并不是所有的PCA结果都是这种明显的分群结果会比较好，还是要根据分析的目的而定。例如在GWAS分析当中，这种“天女散花”一般的PCA散点图，正说明了样本之间不具备明显的亚群分化，适宜进行后续的GWAS分析。
例3：
下图：通过主成分分析方法(PCA)分析9种食物的蛋白质消耗量(变量)与25个欧洲国家(样本)之间的关系
样本点连线距离长 = 样本之间差异性大
样本点连线距离短 = 样本之间差异性小
由图可得，大部分欧洲国家蛋白摄入习惯是：吃鸡蛋、红肉（猪牛羊等畜肉）、白肉（禽、鱼肉及水产品），喝牛奶。
1、各样本点连线的距离：体现各国家蛋白摄入习惯的相似性。
2、主成分与原变量之间的关系：箭头对应的原始变量在投影到水平和垂直方向上后的值，可以分别体现该变量与PC1和PC2的相关性（正负相关性及其大小）（例如，Eggs对PC1具有较大的贡献，而Nuts则与PC1之间呈较大的负相关性）。
3、样本点和箭头之间的距离：反映样本与原始变量的关系。（对于图中用蓝色粗箭头所指的样本点而言，该国的蛋白质来源主要为Fruits and Vegetables）。
例4：
通过主成分分析方法(PCA)分析9种食物的蛋白质消耗量(变量)与25个欧洲国家(样本)之间的关系（3D图）
例3的三维图（取三个主成分的PA图），读图方法不变，只不过是分析样本点在三维空间的距离，而非平面。
例5：
下图：以肠道菌群种类为变量分析肠道菌群相似性
样本间有一定程度的聚类，CC组和HC组的肠道菌群存在明显分离，重叠部分不大（即两个种群相似度不大）。由图可得，结肠癌后肠道菌群发生了一定程度变化。
例6：
下图：以多种蛋白质的表达量为指标区分样本
30例健康组样本点（Healthy）聚集在右上角区域，与75例HCC组和LC组的样本点明显区分。而且HCC组的7个样本点和LC组的2个样本点位于健康者区域，没有健康者被归类为HCC或LC组。
PCA图绘制，R语言代码
数据来源
从TCGA下载结肠癌+直肠癌的表达数据并进行合并，经过配对处理，挑选出44个病人的44个肿瘤样本与44个正常样本。
原始数据"exp_pair"使用count格式，然后经过 dat <- log(exp_pair + 1) 处理形成下图中的数据"dat"。 前44个是正常样本，后44个是肿瘤样本。
dat <- log(exp_pair + 1) # 因为原始数据呈泊松分布，无法画PCA图，故取其log值。
PCA分析函数
做PCA的函数有很多，常用的函数有：princomp，prcomp及rda。
我们在R中输入的数据类型有两类，分别为R mode和Q mode。一般来说数据每一列为一个变量（variable），每一行为一个数据（observation）。其中R mode的数据列数小于行数，是基于列变量的分析；Q mode数据行数小于列数，是基于行数据的分析。如上述表达数据矩阵中样本数小于基因数，属于R mode型数据。
princomp只能用于R mode，它基于协方差(covariance) 或者相关矩阵(correlation) 提取的特征(eigen)并进行特征值分解。默认用法为 x.princomp = princomp(x, cor = FALSE, scores = TRUE)。
Cor是逻辑值的参数，默认cor = FALSE用协方差矩阵计算。cor = TRUE就会用相关矩阵计算特征值。
Princomp 函数的说明文档中推荐 prcomp 更好： The calculation is done using eigen on the correlation or covariance matrix, as determined by cor.
This is done for compatibility with the S-PLUS result.
A preferred method of calculation is to use svd on x, as is done in prcomp.
默认用法为 x.prcomp = prcomp(x, center = TRUE, scale. = FALSE)。
Center为逻辑值，决定是否要以0为中心化，Scale为逻辑值，决定是否要以单位方差进行标准化。
Prcomp 函数的说明文档中说这种方法在数值上更准确：
This is generally the preferred method fornumerical accuracy。
Princomp和prcomp在算法上略有差异。除了分别为特征值分解和奇异值分解外，两者在之前计算协方差的时候标准化的过程存在差异：princomp 计算时分母为N，而 prcomp 分母为N-1。
冗余分析（redundancy analysis，RDA）是一种回归分析结合主成分分析的排序方法，也是多响应变量（multiresponse）回归分析的拓展。从概念上讲，RDA是响应变量矩阵与解释变量之间多元多重线性回归的拟合值矩阵的PCA分析。rda是vegan包的一个函数，可以用rda这个函数来做PCA。虽然简单，但是功能强大。只输入OTU表时做PCA，如果再加上环境因子就做RDA。函数的说明文档中没有专门提做PCA时的方法。但是做RDA采用的是奇异值分解。
princomp和prcomp都是R自带的stats包中的函数。一般用 prcomp 函数即可。
其次，使用 FactoMineR::PCA 也可以做PCA分析。
1. 使用 stats::prcomp 函数进行PCA分析
# 数据转置，转换成行为样本，列为基因的矩阵。prcomp函数要求分析矩阵的列包含特征向量。 data <- t(dat)
# 使用stats::prcomp函数进行PCA分析 prcomp.pca <- stats::prcomp(data) # 查看PCA的分析结果 summary(prcomp.pca) # 结果中PC的数量为nrow和ncol两数中比较小的那个数 # Standard deviation：标准差，其平方为方差 = 特征值 # Proportion of Variance：方差贡献率 # Cumulative Proportion：方差累计贡献率
# 绘制主成分的碎石图 screeplot(prcomp.pca, npcs = 10, type = "lines") # 折线图 screeplot(prcomp.pca, npcs = 6, type = "barplot") # 条形图
# 查看主成分的结果 head(prcomp.pca$x)
2. 使用 FactoMineR::PCA 函数进行PCA分析
# 数据转置，转换成行为样本，列为基因的矩阵。PCA 函数要求分析矩阵的列包含特征向量。 data <- t(dat)
转置结果同上。
# 使用 FactoMineR::PCA 函数进行PCA分析 PCA.pca <- FactoMineR::PCA(data, graph = FALSE) # 查看PCA的分析结果 summary(PCA.pca) # 结果中DIM的数量为nrow和ncol两数中比较小的那个数-1 # Variance：方差， = 特征值 # % of var.：方差贡献率 # Cumulative % of var.：方差累计贡献率 # Dim：Dimension，维度，类似于主成分PC
除了主成分的结果之外，还提供了对于“个体”和“变量”的汇总，即 Individuals，Variables。
Individuals 中，Dist：distance。测量每个个体在多维空间中从原点到样本点的距离。Dim.X：Dimension。测量每个个体在多维空间中从原点到主成分X的距离投影。
ctr：contribution。以百分比形式表示每个个体对给定主成分的贡献。详见 PCA.pca$ind$contrib，每列总和为100%。
cos2：是每个主成分的平方余弦，计算为（Dim.X/Dist）^2。对于给定的主成分，它越接近1，该主成分在捕捉该个体的所有特征方面就越好。
Variables 中，“Dim.X”/“ctr”/“cos2”的解释类似。
# 绘制主成分的碎石图： # FactoMineR::PCA 函数分析的结果用screeplot函数画不出碎石图，总是报错，具体原因不详。但是可以用 factoextra::fviz_screeplot 函数画出来： factoextra::fviz_screeplot(PCA.pca, addlabels = TRUE)
# 查看主成分的结果 head(PCA.pca[["ind"]][["coord"]]) # 结果只显示前5个Dim
3. 使用 graphics 包绘制PCA图
# 使用 graphics 包绘制PCA图
data <- t(dat) # 转置，转换成行为样本，列为基因的矩阵
# 使用stats::prcomp函数进行PCA分析
prcomp.pca <- prcomp(data)
# 使用FactoMineR::PCA函数进行PCA分析
PCA.pca <- FactoMineR::PCA(data, graph = F)
# 绘制主成分的碎石图
screeplot(prcomp.pca, npcs = 10, type = "lines")
factoextra::fviz_screeplot(PCA.pca, addlabels = TRUE)
# 绘制stats::prcomp函数分析的PCA图
plot(prcomp.pca$x, cex = 2, main = "PCA analysis", # cex：点大小。main：标题。prcomp.pca$x是 matrix 形式，plot函数会自动取 matrix 形式数据的前两列，即PC1和PC2。
col = c(rep("blue",44), rep("red",44)), # 点颜色，重复次数根据样本数量改动。这里正常样本
pch = c(rep(16,44), rep(17,44))) # 点形状，见?graphics::points。重复次数根据样本数量改动。
# 绘制FactoMineR::PCA函数分析的PCA图，只修改输入数据即可
plot(PCA.pca[["ind"]][["coord"]], cex = 2, main = "PCA analysis", # cex：点大小。main：标题。prcomp.pca$x是 matrix 形式，plot函数会自动取 matrix 形式数据的前两列，即PC1和PC2。
col = c(rep("blue",44), rep("red",44)), # 点颜色，重复次数根据样本数量改动。这里正常样本
pch = c(rep(16,44), rep(17,44))) # 点形状，见?graphics::points。重复次数根据样本数量改动。
# 添加分隔线
abline(h = 0, v = 0, # h：horizontal line，水平线的y值。v：vertical line，垂直线的x值。
lty = 2, col = "gray") # lty：line type，线类型。2表示虚线，详见graphics::par。col：color，线颜色。
# 给每个点添加标签
text(prcomp.pca$x, labels = rownames(prcomp.pca$x),
pos = 4, # pos：position，1、2、3和4 分别表示指定(x，y)坐标下、左、上和右侧的位置。
offset = 0.6, cex = 0.5) # offset：偏移量。cex：character expansion，字符大小。
# 添加图例
legend("topleft", title = "Sample", inset = 0, cex = 1, # inset：图例距离图中心的距离。cex：图例字体大小
legend = c("normal", "tumor"),
col = c("blue", "red"),
pch = c(16, 17))
prcomp函数分析画的PCA图横纵坐标为PC1、PC2。PCA函数分析画的PCA图横纵坐标为Dim1、Dim2。
但是该基础包的函数画出的图比较粗糙。
4. 使用 ggplot2 包绘制PCA图
# 使用ggplot2包绘制PCA图
data <- t(dat) # 转置，转换成行为样本，列为基因的矩阵
library(ggplot2)
color = c("#2874C5", "#f87669", "#e6b707", "#868686", "#92C5DE",
"#F4A582", "#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3",
"#A6D854", "#FFD92F", "#E5C494", "#B3B3B3") # 设置颜色备选项
# 使用stats::prcomp函数进行PCA分析
prcomp.pca <- prcomp(data)
# 使用FactoMineR::PCA函数进行PCA分析
PCA.pca <- FactoMineR::PCA(data, graph = F)
# 绘制主成分的碎石图
screeplot(prcomp.pca, npcs = 6, type = "barplot")
factoextra::fviz_screeplot(PCA.pca, addlabels = TRUE)
# 绘制PCA图
group_list <- as.factor(c(rep(0,44),rep(1,44))) # 设置分组变量，1为实验组，0为对照组
pdat = data.frame(PCA.pca[["ind"]][["coord"]], Group = group_list) # 准备画图输入数据
p = ggplot(pdat, aes(Dim.1, Dim.2)) + # 画背景、x轴、y轴
geom_point(aes(Dim.1, Dim.2, fill = Group), shape = 21, colour = "black") + # 画点图
scale_color_manual(values = color[1:nlevels(group_list)]) + # 改变点颜色
scale_fill_manual(values = color[1:nlevels(group_list)]) + # 改变点填充颜色
theme_classic() + # 修改外观主题，此处为经典主题，即具有x、y轴线，没有网格线
theme(legend.position = "top") + # 调整图例位置，默认在右侧
labs(fill= "Group") + # 修改图例标题，默认为分组变量的列名
stat_ellipse(aes(color = Group, fill = Group),
geom = "polygon", alpha = 0.3, linetype = 2) # 添加置信椭圆线，color=线条颜色，fill=填充颜色，alpha=椭圆透明度，linetype=线类型
return(p)
5. 使用 factoextra 包绘制PCA图
# 使用 factoextra 包绘制PCA图
data <- t(dat) # 转置，转换成行为样本，列为基因的矩阵
library(factoextra)
color = c("#2874C5", "#f87669", "#e6b707", "#868686", "#92C5DE",
"#F4A582", "#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3",
"#A6D854", "#FFD92F", "#E5C494", "#B3B3B3") # 设置颜色备选项
# 使用stats::prcomp函数进行PCA分析
prcomp.pca <- prcomp(data)
# 使用FactoMineR::PCA函数进行PCA分析
PCA.pca <- FactoMineR::PCA(data, graph = F)
# 绘制主成分的碎石图
screeplot(prcomp.pca, npcs = 6, type = "barplot")
factoextra::fviz_screeplot(PCA.pca, addlabels = TRUE) # 添加百分比标签
# 绘制PCA图
group_list <- as.factor(c(rep(0,44),rep(1,44))) # 设置分组变量，1为实验组，0为对照组
# 样本按group_list分组的点图
factoextra::fviz_pca_ind(PCA.pca, geom.ind = "point", # 样本只展示点图
col.ind = group_list, # 样本颜色按group_list区分
palette = color, # 按组着色或填充的调色板
addEllipses = TRUE, # 添加置信椭圆
legend.title = "Groups") # 图例标题
factoextra::fviz_pca_ind(PCA.pca, label = "none", # 文本指定要标记的元素，ind标记样本名，var标记变量名，all标记样本和变量，none不做标记
habillage = group_list, # 按组对观测值进行着色
addEllipses = TRUE, ellipse.level = 0.95, # 椭圆置信度
palette = c("#999999", "#E69F00", "#56B4E9"))
# 样本按cos^2分组的点图
# 对于给定的主成分，cos^2越接近1，该主成分在捕捉该个体的所有特征方面就越好。
factoextra::fviz_pca_ind(PCA.pca, geom = "point", # 样本只展示点图
col.ind ="cos2", # 样本颜色按cos^2区分
gradient.cols = c("white", "#2E9FDF", "#FC4E07")) # 梯度颜色
# 变量图，展示各变量箭头，代表主成分载荷
factoextra::fviz_pca_var(PCA.pca, col.var = "steelblue")
# 变量图，以百分比形式展示各变量对给定主成分的贡献
factoextra::fviz_pca_var(PCA.pca, col.var = "contrib", # 按贡献度大小区分颜色
gradient.cols = c("white", "blue", "red"), # 梯度颜色
ggtheme = theme_minimal()) # 没有背景注释的简约主题
# 样本和变量的双标图
factoextra::fviz_pca_biplot(PCA.pca, label = "var", # lable：文本指定要标记的元素，ind标记样本名，var标记变量名，all标记样本和变量，none不做标记
habillage = group_list, # 按组对观测值进行着色
addEllipses = TRUE, ellipse.level = 0.95, # 椭圆置信度
ggtheme = theme_minimal()) # 没有背景注释的简约主题
样本按group_list分组的点图：
对于给定的主成分，cos^2越接近1，该主成分在捕捉该个体的所有特征方面就越好。
变量图，展示各变量箭头，代表主成分载荷：
变量图，以百分比形式展示各变量对给定主成分的贡献：
样本和变量的双标图：
6. 使用 scatterplot3d 包绘制PCA图
# 使用 scatterplot3d 包绘制PCA图
data <- t(dat) # 转置，转换成行为样本，列为基因的矩阵
library(scatterplot3d)
color = c("#2874C5", "#f87669", "#e6b707", "#868686", "#92C5DE",
"#F4A582", "#66C2A5", "#FC8D62", "#8DA0CB", "#E78AC3",
"#A6D854", "#FFD92F", "#E5C494", "#B3B3B3")
# 使用stats::prcomp函数进行PCA分析
prcomp.pca <- prcomp(data)
# 使用FactoMineR::PCA函数进行PCA分析
PCA.pca <- FactoMineR::PCA(data, graph = F)
# 绘制主成分的碎石图
screeplot(prcomp.pca, npcs = 6, type = "barplot")
factoextra::fviz_screeplot(PCA.pca, addlabels = TRUE) # 添加百分比标签
group_list <- as.factor(c(rep(0,44),rep(1,44))) # 设置分组变量，1为实验组，0为对照组
colors = color[as.numeric(group_list)] # 为每个样本点设置颜色
pdat = data.frame(PCA.pca[["ind"]][["coord"]], Group = group_list) # 准备画图输入数据
# 画PCA函数的点图
scatterplot3d::scatterplot3d(pdat[, 1:3], main = "PCA", # 取前3位主成分，main：标题
pch = 21, cex.symbols = 1.5, color = "black", bg = colors) # pch：点形状。color：点外圈颜色。bg：点填充颜色
# 给PCA函数的点图添加图例
graphics::legend("topright", # 图例位置，有"bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center"
legend = c("normal", "tumor"), # 图例文本信息来源
pch = 21, col = "black", pt.bg = color[1:nlevels(group_list)], # pch：点形状。col：点外圈颜色。pt.bg：点填充颜色
xpd = TRUE, # FALSE为不允许在作图区域外作图，改为TRUE即可
horiz = FALSE, # 图例放置方向，TRUE为水平放置图例，FALSE为垂直放置
inset = -0.09, # 图例距离图中心的距离
bty = "o") # 图例框是否画出，"o"为画出，"n"为不画出
# 画prcomp函数的点图
scatterplot3d::scatterplot3d(prcomp.pca$x[,1:3], main = "3D PCA plot", # 取前3位主成分。main：标题
pch = c(rep(16,44), rep(17,44)), # 点形状
color = c(rep("blue",44), rep("red",44)), # 点颜色
angle = 45, # angle：观察角度
cex.symbols = 1.5, # cex.symbols：点放大倍率
mar = c(5,4,4,5)) # c(bottom, left, top, right)的数值向量，设置立方体底面、左面、上面、右面到图边缘的距离。
# 给prcomp函数的点图添加图例
graphics::legend("topright", # 图例位置，有"bottomright", "bottom", "bottomleft", "left", "topleft", "top", "topright", "right" and "center"
title = "Sample", # 图例标题
xpd = TRUE, # FALSE为不允许在作图区域外作图，改为TRUE即可
inset= -0.11, # 图例距离图中心的距离
legend = c("normal","tumor"), # 图例文本信息来源
col = c("blue","red"), # 图例点颜色
pch = c(16,17)) # 点形状
