首先,必须说明,PCA总体来说是主观的,二维不行的时候上三维。
那么先简单看一下PCA是什么
PCA是使用最广泛的数据降维算法之一,利用降维思想,把多个指标转换成少数几个综合指标。其主要思想是将n维特征映射到k维上,低维代替高维,但损失信息很少。
在代数上,它是将原随机向量的协方差矩阵变换成对角形针;在几何上,它是一个线性变换,把数据变换到一个新的坐标系统中,使得任何数据投影的第一大方差在第一个坐标(第一主成分)上,第二大方差在第二个坐标(第二主成分)上,以此类推。它减少数据集的维数,同时保持数据集的方差贡献最大的特征。
PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴,新的坐标轴的选择与数据本身是密切相关的。其中,第一个新坐标轴选择是原始数据中方差最大的方向,第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的,第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推,可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴,我们发现,大部分方差都包含在前面k个坐标轴中,后面的坐标轴所含的方差几乎为0。于是,我们可以忽略余下的坐标轴,只保留前面k个含有绝大部分方差的坐标轴。事实上,这相当于只保留包含绝大部分方差的维度特征,而忽略包含方差几乎为0的特征维度,实现对数据特征的降维处理。
然而,实际问题会涉及众多变量,每个变量不同程度反映研究问题的特征,并且变量之间存在一定相关性。既然变量之间存在相关性,那么就必然存在起支配作用的变量。同时,变量反映的信息在一定程度上有重叠,如果不去冗余,会增加问题分析的复杂性,但更一般的情形是,某一些变量的线性组合才会差异太小,如果要降低维度,消减多余的属性,我们必须要找到这种组合,即将原来互相相关的变量重新组合成一组新的相互无关的几个综合变量。
时刻牢记:差异才是信息!!
步骤如下:
去除平均值
计算协方差矩阵
计算协方差矩阵的特征值和特征向量
将特征值排序
保留前N个最大的特征值对应的特征向量
将原始特征转换到上面得到的N个特征向量构建的新空间中
(最后两步,实现了特征压缩)
一句话概括,要对一批样本进行降维,需要先对所有的属性进行归一化的减均值处理,然后求其协方差矩阵的特征向量,将特征值按从大到小的顺序排列,特征值越大新基对应的新样本属性就越重要。最后我们就可以按照需要舍弃最后面特征值较小对应的特征向量作为新基下投影的样本属性了。
PCA的主要作用
1.降低所研究的数据空间的维数,删除多余变量。
2.有时可通过因子负荷aij的结论,弄清X变量间的某些关系。
3.多维数据的可视化,由图形可直观地看出各样品在主分量中的地位,进而还可以对样本进行分类处理,可以由图形发现远离大多数样本点的离群点。
4.构造回归模型,即把各主成分作为新自变量代替原来自变量x做回归分析。
5.筛选回归变量。回归变量的选择有着重的实际意义,为使模型本身易于做结构分析、控制和预报,好从原始变量所构成的子集合中选择最佳变量,构成最佳变量集合。用主成分分析筛选变量,可以用较少的计算量来选择量,获得选择最佳变量子集合的效果。
当然对我们生物和医学来说,最重要的还是怎么简单实操,毕竟拿到rawdata我们想要的是快速获得table,而对table的分析才是重点(我导箴言)
这里PCA分析通过prcomp包,绘图通过ggplot、ggrepel和scatterplot3d完成
代码实现
#二维呈现
library(ggplot2)
setwd("数据所在路径")
#读取数据集
rt= read.table("expression.txt",header = T,sep = "\t",check.names = F,row.names = 1)
data=rt[,c(1:ncol(rt)-1)]
#进行PCA分析,scale.=TRUE表示分析前对数据进行归一化
df_pca=prcomp(data,scale. = TRUE)
df_pcs<-data.frame(df_pca$x,sample=rt$sample)
head(df_pcs,3)
percentage<-round(df_pca$sdev / sum(df_pca$sdev)*100,2)
percentage<-paste(colnames(df_pcs),"(",paste(as.character(percentage),"%",")",sep = "\t"))
library(ggrepel)
p<-ggplot(df_pcs,aes(x=PC1,y=PC2,color=sample))+
geom_point(size=3.5)+
#stat_ellipse(aes(fill=rt$sample),alpha=0.2,level = 0.95,show.legend = F,geom = "polygon")+
xlab(percentage[1])+
ylab(percentage[2])+
geom_text_repel(aes(x=PC1,y=PC2, label = rownames(data))) +
#给每个散点加标签,也可以直接用geom_text,不需要加载ggrepel包,但是常出现标签重叠不清晰。
theme_classic()+
#annotate('text',label = 'control',x=2.5,y=-2.5,size=5,colour='#F8766D')+
#annotate('text',label = 'treat',x=2.5,y=1,size=5,colour='#28B3B6')+
labs(title = "PCA clustering")+
png(file="PCA2.png",bg="transparent",width = 2000,height = 1500,res = 250,units"px")
#transparent表示透明画布,bg(background)
#如果重复数太少,或者一个组别里的样本量太少,是没有办法计算置信椭圆的,把 stat_ellipse一行去掉
#当上述情况发生时,可以用ggalt包中的geom_encircle把样品包起来,https://blog.csdn.net/qazplm12_3/article/details/120620477
print(p)
dev.off()
pdf(file = paste("PCA.pdf",sep = ""),height = 6,width = 8,onefile = FALSE)
print(p)
dev.off()
####三维呈现
library(scatterplot3d)
setwd("数据所在路径")
#读取数据集
rt= read.table("expression.txt",header = T,sep = "\t",check.names = F,row.names = 1)
data=rt[,c(1:ncol(rt)-1)]
#进行PCA分析,scale.=TRUE表示分析前对数据进行归一化
df_pca=prcomp(data,scale. = TRUE)
df_pcs<-data.frame(df_pca$x,sample=rt$sample)
head(df_pcs,3)
#group=as.vector(rt[,"sample"])
##按组更改点形状和颜色,哪个不要就注释掉,同时注意sactterplot3d中也要相应删除
shapes=c(16,17,18,19,20,21,22,23) #数字代表不同的形状
shapes <- shapes[as.factor(rt$sample)] #as.xxx,具体根据sample列修改,文字是factor,数字是numeric
shapes #用来检查一下组别是否正确
colors<-c("#28B3B6","#F8766D","#BA5952","#55988D","#14147E","#FFB528","#C4FF90","#7079FF")
colors<-colors[as.factor(rt$sample)]
colors
pdf(file=paste("PCA3D.pdf",sep=""),height = 6,width = 6.5,onefile = FALSE)
s3d<-scatterplot3d(df_pcs[,1:3],pch=shapes,angle=45,color=colors,grid=FALSE,box=FALSE,type = "h",
main="3D_PCA_Plot",
xlab="PC1",
ylab="PC2",
zlab="PC3") #添加图标题和x,y,z标签
'''
#grid和box:逻辑值,如果均为TRUE,则在底部和箱体绘制网格和一个框,想要删除改为FALSE即可。
#x、y 和 z 是指定点的 x、y、z 坐标的数值向量。
x 可以是一个矩阵或一个包含 3 列对应于 x、y 和 z 坐标的数据框。
在这种情况下,参数 y 和 z 是可选的
#grid 指定绘制网格的面。
可能的值是“xy”、“xz”或“yz”的组合。 示例:grid = c(“xy”, “yz”)。
#默认值为 TRUE 以仅在 xy 平面上添加网格。
col.grid, lty.grid: 用于网格的颜色和线型
type="h"可以给每个点添加x-y平面向上的bar
'''
#source('http://www.sthda.com/sthda/RDoc/functions/addgrids3d.r') #source源函数addgrids3d
#addgrids3d(df_pcs[, 1:3], grid = c("xy", "xz", "yz")) #若要添加不同平面网格线,需要先将scatterplot3d中grid和box改为FALSE。
#这两步其实会导致在散点的上面画网格,要想网格线在散点下方还需前置以下:
'''
# 1. 源函数
source('~/hubiC/Documents/R/function/addgrids3d.r')
# 2. 使用 pch="" 清空 3D 散点图,创建一个空的图形,把scatterplot3d的结果指向s3d。
s3d <- scatterplot3d(df_pcs[, 1:3], pch = "", grid=FALSE, box=FALSE)
# 3. 添加网格
addgrids3d(df_pcs[, 1:3], grid = c("xy", "xz", "yz"))
# 4. 添加点
s3d$points3d(df_pcs[, 1:3], pch = )
'''
#legend(s3d$xyz.convert(7.5,3,4.5),legend = levels(rt$sample),pch=shapes,col = colors)
#添加图例,xyz.convert用来制定图例坐标,也可以用关键字“top”,“bottom”,“right”,“topleft”等指定。
text(s3d$xyz.convert(df_pcs[,1:3]),labels=rownames(rt),cex=0.7,col = "black") #给散点添加标签
#添加回归平面见https://blog.csdn.net/m0_49960764/article/details/122249790
dev.off()
记录这篇主要是因为,我的数据是多个分组的,每个分组只有两个样本,因此在绘图时我希望能够把所有分组分别标注,但是在网上找到的大多数教程都是两分组多重复,类似“control”和“treat”,所以就摸索了一下,顺便熟悉一下散点图绘制,记录一下方便以后用~