【數(shù)據(jù)降維與聚類】 之 主成分分析(PCA)及可視化 “全網(wǎng)最詳細V烧铡蹂空!”教程

首先,必須說明果录,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”汁展,所以就摸索了一下,順便熟悉一下散點圖繪制厌殉,記錄一下方便以后用~

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末食绿,一起剝皮案震驚了整個濱河市,隨后出現(xiàn)的幾起案子公罕,更是在濱河造成了極大的恐慌器紧,老刑警劉巖,帶你破解...
    沈念sama閱讀 206,013評論 6 481
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件楼眷,死亡現(xiàn)場離奇詭異铲汪,居然都是意外死亡熊尉,警方通過查閱死者的電腦和手機,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評論 2 382
  • 文/潘曉璐 我一進店門掌腰,熙熙樓的掌柜王于貴愁眉苦臉地迎上來狰住,“玉大人,你說我怎么就攤上這事齿梁〈咧玻” “怎么了?”我有些...
    開封第一講書人閱讀 152,370評論 0 342
  • 文/不壞的土叔 我叫張陵勺择,是天一觀的道長创南。 經(jīng)常有香客問我,道長省核,這世上最難降的妖魔是什么稿辙? 我笑而不...
    開封第一講書人閱讀 55,168評論 1 278
  • 正文 為了忘掉前任,我火速辦了婚禮气忠,結(jié)果婚禮上邻储,老公的妹妹穿的比我還像新娘。我一直安慰自己笔刹,他們只是感情好芥备,可當(dāng)我...
    茶點故事閱讀 64,153評論 5 371
  • 文/花漫 我一把揭開白布冬耿。 她就那樣靜靜地躺著舌菜,像睡著了一般。 火紅的嫁衣襯著肌膚如雪亦镶。 梳的紋絲不亂的頭發(fā)上日月,一...
    開封第一講書人閱讀 48,954評論 1 283
  • 那天,我揣著相機與錄音缤骨,去河邊找鬼爱咬。 笑死,一個胖子當(dāng)著我的面吹牛绊起,可吹牛的內(nèi)容都是我干的精拟。 我是一名探鬼主播,決...
    沈念sama閱讀 38,271評論 3 399
  • 文/蒼蘭香墨 我猛地睜開眼虱歪,長吁一口氣:“原來是場噩夢啊……” “哼蜂绎!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起笋鄙,我...
    開封第一講書人閱讀 36,916評論 0 259
  • 序言:老撾萬榮一對情侶失蹤师枣,失蹤者是張志新(化名)和其女友劉穎,沒想到半個月后萧落,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體践美,經(jīng)...
    沈念sama閱讀 43,382評論 1 300
  • 正文 獨居荒郊野嶺守林人離奇死亡洗贰,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 35,877評論 2 323
  • 正文 我和宋清朗相戀三年,在試婚紗的時候發(fā)現(xiàn)自己被綠了陨倡。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片敛滋。...
    茶點故事閱讀 37,989評論 1 333
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖兴革,靈堂內(nèi)的尸體忽然破棺而出矛缨,到底是詐尸還是另有隱情,我是刑警寧澤帖旨,帶...
    沈念sama閱讀 33,624評論 4 322
  • 正文 年R本政府宣布箕昭,位于F島的核電站,受9級特大地震影響解阅,放射性物質(zhì)發(fā)生泄漏落竹。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點故事閱讀 39,209評論 3 307
  • 文/蒙蒙 一货抄、第九天 我趴在偏房一處隱蔽的房頂上張望述召。 院中可真熱鬧,春花似錦蟹地、人聲如沸积暖。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,199評論 0 19
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽夺刑。三九已至,卻和暖如春分别,著一層夾襖步出監(jiān)牢的瞬間遍愿,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 31,418評論 1 260
  • 我被黑心中介騙來泰國打工耘斩, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留沼填,地道東北人。 一個月前我還...
    沈念sama閱讀 45,401評論 2 352
  • 正文 我出身青樓括授,卻偏偏與公主長得像坞笙,于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子荚虚,可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 42,700評論 2 345

推薦閱讀更多精彩內(nèi)容