首先,必須說明果录,PCA總體來說是主觀的上枕,二維不行的時候上三維。
那么先簡單看一下PCA是什么
PCA是使用最廣泛的數(shù)據(jù)降維算法之一弱恒,利用降維思想辨萍,把多個指標(biāo)轉(zhuǎn)換成少數(shù)幾個綜合指標(biāo)。其主要思想是將n維特征映射到k維上返弹,低維代替高維锈玉,但損失信息很少。
在代數(shù)上义起,它是將原隨機向量的協(xié)方差矩陣變換成對角形針拉背;在幾何上,它是一個線性變換默终,把數(shù)據(jù)變換到一個新的坐標(biāo)系統(tǒng)中椅棺,使得任何數(shù)據(jù)投影的第一大方差在第一個坐標(biāo)(第一主成分)上,第二大方差在第二個坐標(biāo)(第二主成分)上齐蔽,以此類推两疚。它減少數(shù)據(jù)集的維數(shù),同時保持?jǐn)?shù)據(jù)集的方差貢獻最大的特征含滴。
PCA的工作就是從原始的空間中順序地找一組相互正交的坐標(biāo)軸诱渤,新的坐標(biāo)軸的選擇與數(shù)據(jù)本身是密切相關(guān)的。其中谈况,第一個新坐標(biāo)軸選擇是原始數(shù)據(jù)中方差最大的方向源哩,第二個新坐標(biāo)軸選取是與第一個坐標(biāo)軸正交的平面中使得方差最大的,第三個軸是與第1,2個軸正交的平面中方差最大的鸦做。依次類推励烦,可以得到n個這樣的坐標(biāo)軸。通過這種方式獲得的新的坐標(biāo)軸泼诱,我們發(fā)現(xiàn)坛掠,大部分方差都包含在前面k個坐標(biāo)軸中,后面的坐標(biāo)軸所含的方差幾乎為0。于是屉栓,我們可以忽略余下的坐標(biāo)軸舷蒲,只保留前面k個含有絕大部分方差的坐標(biāo)軸。事實上友多,這相當(dāng)于只保留包含絕大部分方差的維度特征牲平,而忽略包含方差幾乎為0的特征維度,實現(xiàn)對數(shù)據(jù)特征的降維處理域滥。
然而纵柿,實際問題會涉及眾多變量,每個變量不同程度反映研究問題的特征启绰,并且變量之間存在一定相關(guān)性昂儒。既然變量之間存在相關(guān)性,那么就必然存在起支配作用的變量委可。同時渊跋,變量反映的信息在一定程度上有重疊,如果不去冗余着倾,會增加問題分析的復(fù)雜性拾酝,但更一般的情形是,某一些變量的線性組合才會差異太小卡者,如果要降低維度蒿囤,消減多余的屬性,我們必須要找到這種組合虎眨,即將原來互相相關(guān)的變量重新組合成一組新的相互無關(guān)的幾個綜合變量蟋软。
時刻牢記:差異才是信息!嗽桩!
步驟如下:
去除平均值
計算協(xié)方差矩陣
計算協(xié)方差矩陣的特征值和特征向量
將特征值排序
保留前N個最大的特征值對應(yīng)的特征向量
將原始特征轉(zhuǎn)換到上面得到的N個特征向量構(gòu)建的新空間中
(最后兩步岳守,實現(xiàn)了特征壓縮)
一句話概括,要對一批樣本進行降維碌冶,需要先對所有的屬性進行歸一化的減均值處理湿痢,然后求其協(xié)方差矩陣的特征向量,將特征值按從大到小的順序排列扑庞,特征值越大新基對應(yīng)的新樣本屬性就越重要譬重。最后我們就可以按照需要舍棄最后面特征值較小對應(yīng)的特征向量作為新基下投影的樣本屬性了。
PCA的主要作用
1.降低所研究的數(shù)據(jù)空間的維數(shù)罐氨,刪除多余變量臀规。
2.有時可通過因子負(fù)荷aij的結(jié)論,弄清X變量間的某些關(guān)系栅隐。
3.多維數(shù)據(jù)的可視化塔嬉,由圖形可直觀地看出各樣品在主分量中的地位玩徊,進而還可以對樣本進行分類處理,可以由圖形發(fā)現(xiàn)遠離大多數(shù)樣本點的離群點谨究。
4.構(gòu)造回歸模型恩袱,即把各主成分作為新自變量代替原來自變量x做回歸分析。
5.篩選回歸變量胶哲∨纤回歸變量的選擇有著重的實際意義,為使模型本身易于做結(jié)構(gòu)分析鸯屿、控制和預(yù)報澈吨,好從原始變量所構(gòu)成的子集合中選擇最佳變量,構(gòu)成最佳變量集合碾盟。用主成分分析篩選變量棚辽,可以用較少的計算量來選擇量技竟,獲得選擇最佳變量子集合的效果冰肴。
當(dāng)然對我們生物和醫(yī)學(xué)來說,最重要的還是怎么簡單實操榔组,畢竟拿到rawdata我們想要的是快速獲得table熙尉,而對table的分析才是重點(我導(dǎo)箴言)
這里PCA分析通過prcomp包,繪圖通過ggplot搓扯、ggrepel和scatterplot3d完成
代碼實現(xiàn)
#二維呈現(xiàn)
library(ggplot2)
setwd("數(shù)據(jù)所在路徑")
#讀取數(shù)據(jù)集
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表示分析前對數(shù)據(jù)進行歸一化
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))) +
#給每個散點加標(biāo)簽,也可以直接用geom_text,不需要加載ggrepel包锨推,但是常出現(xiàn)標(biāo)簽重疊不清晰铅歼。
? 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)
#如果重復(fù)數(shù)太少换可,或者一個組別里的樣本量太少椎椰,是沒有辦法計算置信橢圓的,把 stat_ellipse一行去掉
#當(dāng)上述情況發(fā)生時沾鳄,可以用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()
####三維呈現(xiàn)
library(scatterplot3d)
setwd("數(shù)據(jù)所在路徑")
#讀取數(shù)據(jù)集
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表示分析前對數(shù)據(jù)進行歸一化
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中也要相應(yīng)刪除
shapes=c(16,17,18,19,20,21,22,23) #數(shù)字代表不同的形狀
shapes <- shapes[as.factor(rt$sample)] #as.xxx,具體根據(jù)sample列修改,文字是factor,數(shù)字是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")? #添加圖標(biāo)題和x吞歼,y圈膏,z標(biāo)簽
'''
? ? #grid和box:邏輯值,如果均為TRUE篙骡,則在底部和箱體繪制網(wǎng)格和一個框稽坤,想要刪除改為FALSE即可桥帆。
? ? #x、y 和 z 是指定點的 x慎皱、y老虫、z 坐標(biāo)的數(shù)值向量。
? ? x 可以是一個矩陣或一個包含 3 列對應(yīng)于 x茫多、y 和 z 坐標(biāo)的數(shù)據(jù)框祈匙。
? ? 在這種情況下,參數(shù) y 和 z 是可選的
? ? #grid 指定繪制網(wǎng)格的面天揖。
? ? 可能的值是“xy”夺欲、“xz”或“yz”的組合。 示例:grid = c(“xy”, “yz”)今膊。
? ? #默認(rèn)值為 TRUE 以僅在 xy 平面上添加網(wǎng)格些阅。
? ? col.grid, lty.grid: 用于網(wǎng)格的顏色和線型
? ? type="h"可以給每個點添加x-y平面向上的bar
'''
#source('http://www.sthda.com/sthda/RDoc/functions/addgrids3d.r')? #source源函數(shù)addgrids3d
#addgrids3d(df_pcs[, 1:3], grid = c("xy", "xz", "yz"))? #若要添加不同平面網(wǎng)格線,需要先將scatterplot3d中g(shù)rid和box改為FALSE斑唬。
#這兩步其實會導(dǎo)致在散點的上面畫網(wǎng)格市埋,要想網(wǎng)格線在散點下方還需前置以下:
'''
# 1. 源函數(shù)
source('~/hubiC/Documents/R/function/addgrids3d.r')
# 2. 使用 pch="" 清空 3D 散點圖,創(chuàng)建一個空的圖形恕刘,把scatterplot3d的結(jié)果指向s3d缤谎。
s3d <- scatterplot3d(df_pcs[, 1:3], pch = "", grid=FALSE, box=FALSE)
# 3. 添加網(wǎng)格
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用來制定圖例坐標(biāo)褐着,也可以用關(guān)鍵字“top”坷澡,“bottom”,“right”含蓉,“topleft”等指定频敛。
text(s3d$xyz.convert(df_pcs[,1:3]),labels=rownames(rt),cex=0.7,col = "black") #給散點添加標(biāo)簽
#添加回歸平面見https://blog.csdn.net/m0_49960764/article/details/122249790
dev.off()
記錄這篇主要是因為,我的數(shù)據(jù)是多個分組的馅扣,每個分組只有兩個樣本斟赚,因此在繪圖時我希望能夠把所有分組分別標(biāo)注,但是在網(wǎng)上找到的大多數(shù)教程都是兩分組多重復(fù)岂嗓,類似“control”和“treat”汁展,所以就摸索了一下,順便熟悉一下散點圖繪制厌殉,記錄一下方便以后用~