【数据降维与聚类】 之 主成分分析(PCA)及可视化 “全网最详细!!”教程

首先,必须说明,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”,所以就摸索了一下,顺便熟悉一下散点图绘制,记录一下方便以后用~

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 203,098评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,213评论 2 380
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 149,960评论 0 336
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,519评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,512评论 5 364
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,533评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,914评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,574评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,804评论 1 296
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,563评论 2 319
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,644评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,350评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,933评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,908评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,146评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,847评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,361评论 2 342

推荐阅读更多精彩内容