定义：
主成分分析(Principal Component Analysis,PCA)，是一种掌握事物主要矛盾的统计分析方法，它可以从多元事物中解析出主要影响因素，揭示事物的本质，简化复杂的问题。
主成分分析试图在力保数据信息丢失最少的原则下，用较少的综合变量代替原本较多的变量，而且综合变量间互不相关。
PCA目标：
PCA 的目标是寻找 r （ r<n ）个新变量，使它们反映事物的主要特征，压缩原有数据矩阵的规模。每个新变量是原有变量的线性组合，体现原有变量的综合效果，具有一定的实际含义。这 r 个新变量称为“主成分”，它们可以在很大程度上反映原来 n 个变量的影响，并且这些新变量是互不相关的，也是正交的。
通过采用这样的主成分，便可以只选用若干变量而不是上千的变量来对一种样品进行分析了。这样，就可以将样品有关变量绘制成图，使得样品间的相似性和相异之处一目了然，对不同样品是否可以归为一组，也一清二楚。
生命科学中的PCA应用：
使用在分析复杂的多维数据集的时候。例如不同实验条件下的转录组测序数据，表达谱芯片数据，以及蛋白组和代谢组数据。
当变量的数目比样品的数目多时，PCA可以在不损失信息量的情况下将样品的维度最大程度地减少至样品数。它可以被看做复杂实验数据预处理的一个步骤。
大概流程：
PCA计算过程：
输入示例：
特征向量：
作图：
每个样品对应pc的坐标值：
作2d图：
3d图：
具体步骤：
首先就是输入文件（流程里提到过格式）：
代码：
# 安装所需的R包（如未安装）
packages <- c("ggplot2", "plotly", "rgl", "vegan")
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)
# 加载所需库
library(ggplot2)
library(plotly)
library(rgl)
library(vegan)
# 设置工作目录
setwd("C:/Users/admin/Desktop/4_复试所有有关文件/生信复试内容/pca主成分分析")
# 读取数据
data <- read.table("input.txt", header = TRUE, sep = "\t", row.names = 1, check.names = FALSE)
data <- t(as.matrix(data)) # 矩阵转置，横坐标变为样品名，纵坐标变为基因名
# PCA分析
data.pca <- prcomp(data, scale. = TRUE)
# 保存PCA结果
write.table(data.pca$rotation, file = "PC.xls", quote = FALSE, sep = "\t") # 特征向量
write.table(predict(data.pca), file = "newTab.xls", quote = FALSE, sep = "\t") # 新表
pca.sum <- summary(data.pca)# 对pca进行概括，看每个pc所占的权重
write.table(pca.sum$importance, file = "importance.xls", quote = FALSE, sep = "\t") # PC比重
# 画PCA柱状图
pdf(file = "pcaBarplot.pdf", width = 15)
barplot(pca.sum$importance[2, ] * 100, xlab = "PC", ylab = "Percent", col = "skyblue")
dev.off()
# 画PCA碎石图
pdf(file = "pcaPlot.pdf", width = 15)
plot(pca.sum$importance[2, ] * 100, type = "o", col = "red", xlab = "PC", ylab = "Percent")
dev.off()
# 设定分组（需自行修改）
group <- c(rep("con", 5), rep("A", 5), rep("B", 3)) # 修改为你的分组信息
# 计算PCA预测值
pcaPredict <- predict(data.pca)
PCA <- data.frame(PCA1 = pcaPredict[, 1], PCA2 = pcaPredict[, 2], group = group)
PCA.mean <- aggregate(PCA[, 1:2], list(group = PCA$group), mean)
# 定义椭圆函数
veganCovEllipse <- function(cov, center = c(0, 0), scale = 1, npoints = 100) {
theta <- seq(0, 2 * pi, length.out = npoints)
Circle <- cbind(cos(theta), sin(theta))
as.data.frame(t(center + scale * t(Circle %*% chol(cov))))
}
# 计算每组的椭圆区域
df_ell <- data.frame()
for (g in unique(PCA$group)) {
group_data <- PCA[PCA$group == g, 1:2] # 选取该组数据
if (nrow(group_data) > 1) { # 至少有两个点才能计算协方差
ell_data <- veganCovEllipse(
cov = cov.wt(group_data, wt = rep(1 / nrow(group_data), nrow(group_data)))$cov,
center = colMeans(group_data)
)
ell_data <- as.data.frame(ell_data)
ell_data$group <- g
df_ell <- rbind(df_ell, ell_data)
}
}
# 确保df_ell的列名为PCA1和PCA2
colnames(df_ell) <- c("PCA1", "PCA2", "group")
# 画PCA 2D图
pdf(file = "PCA2d.pdf")
ggplot(data = PCA, aes(PCA1, PCA2, color = group)) +
geom_point(size = 3) +
geom_polygon(data = df_ell, aes(x = PCA1, y = PCA2, fill = group), alpha = 0.2) + # 用polygon填充椭圆区域
geom_text(data = PCA.mean, aes(x = PCA1, y = PCA2, label = group), vjust = -1, size = 5) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
dev.off()
# 3D PCA可视化（使用plotly替代pca3d）
pca_df <- data.frame(PC1 = pcaPredict[,1], PC2 = pcaPredict[,2], PC3 = pcaPredict[,3], Group = group)
fig <- plot_ly(pca_df, x = ~PC1, y = ~PC2, z = ~PC3, color = ~Group, colors = c("red", "blue", "green"),
type = "scatter3d", mode = "markers")
# 保存3D PCA图
htmlwidgets::saveWidget(fig, "pca3d.html")
输出文件：
不展示了，跟前面的一样
