導(dǎo)讀
主成分分析(Principal Components Analysis拦耐,PCA),也稱主分量分析或主成分回歸分析法狡耻,是一種無(wú)監(jiān)督的數(shù)據(jù)降維方法淮悼。PCA通過(guò)線性變換將原始數(shù)據(jù)變換為一組各維度線性無(wú)關(guān)的表示,可用于提取數(shù)據(jù)的主要特征分量,常用于高維數(shù)據(jù)的降維悉抵。這種降維的思想首先減少數(shù)據(jù)集的維數(shù)肩狂,同時(shí)還保持?jǐn)?shù)據(jù)集的對(duì)方差貢獻(xiàn)最大的特征,最終使數(shù)據(jù)直觀呈現(xiàn)在二維坐標(biāo)系姥饰。
直觀上傻谁,第一主成分軸 優(yōu)于 第二主成分軸,具有最大可分性列粪。
主坐標(biāo)分析(Principal Coordinates Analysis审磁,PCoA),即經(jīng)典多維標(biāo)度(Classical multidimensional scaling)岂座,用于研究數(shù)據(jù)間的相似性态蒂。
【二者差異】
PCA與PCoA都是降低數(shù)據(jù)維度的方法,但是差異在在于PCA是基于原始矩陣费什,而PCoA是基于通過(guò)原始矩陣計(jì)算出的距離矩陣吃媒。因此,PCA是盡力保留數(shù)據(jù)中的變異讓點(diǎn)的位置不改動(dòng)吕喘,而PCoA是盡力保證原本的距離關(guān)系不發(fā)生改變赘那,也就是使得原始數(shù)據(jù)間點(diǎn)的距離與投影中即結(jié)果中各點(diǎn)之間的距離盡可能相關(guān)。
一氯质、整理數(shù)據(jù)(轉(zhuǎn)錄組的基因表達(dá)量數(shù)據(jù))
基因表達(dá)量數(shù)據(jù)通過(guò)RSEM軟件定量后得到
# RSEM構(gòu)建表達(dá)矩陣
rsem-generate-data-matrix Your_Data_Name*.genes.results > output.matrix
# output.matrix是沒(méi)有標(biāo)準(zhǔn)化的表達(dá)量矩陣,原始count值的表達(dá)量矩陣
# 整理表頭和初步篩選數(shù)據(jù)
awk 'BEGIN{printf "geneid\tI0-1\tI0-2\tI0-3\tI0-4\tI3-1\tI3-2\tI3-3\tI3-4\tU0-1\tU0-2\tU0-3\tU0-4\tU3-1\tU3-2\tU3-3\tU3-4\n"}{if($2+$3+$4+$5+$6+$7+$8+$9+$10+$11+$12+$13+$14+$15+$16+$17>10)print $0}' output.matrix > input.txt
# 將input.txt 導(dǎo)出
二募舟、PCA分析及R包繪圖(2D散點(diǎn)圖)
library(ggplot2)
setwd("D:/Your/Working/Directory/")
## 載入數(shù)據(jù)與整理,目的是得到取整后的基因表達(dá)量矩陣(注意:原始矩陣中有0則應(yīng)該整體加1處理)
rsem_counts <- read.table("input.txt",header = T,row.names = 1)
rsem_counts <- rsem_counts+1
rsem_counts <- as.matrix(rsem_counts)
data <- round(rsem_counts,digits = 0)
## 整理分組信息(注意:根據(jù)自己的實(shí)際情況來(lái)添加分組信息闻察,切勿照抄)
period <- factor(c(rep("DPI0",4),rep("DPI3",4),rep("DPI0",4),rep("DPI3",4)))
condition <- factor(c(rep("infection",8),rep("control",8)))
#batch <- factor(c(rep("first",2),rep("second",3),rep("first",1),rep("second",2),rep("first",2),rep("second",2),rep("first",2),rep("second",2)))
##如有批次效應(yīng)拱礁,可添加batch信息
coldata <-data.frame(row.names = colnames(data),condition,period)
## 差異表達(dá)分析的批次校正,但不改變?cè)季仃嚕―ESeqDataSetFormMatrix)
dds <-DESeqDataSetFromMatrix(countData = data,colData = coldata,design = ~period+condition)
## PCA分析
## 樣本量超過(guò)30推薦用vsd,樣本量低于30用rld
#vsd <- vst(dds, blind = FALSE)
#vsd_period <- plotPCA(vsd, intgroup=c("condition","period"),returnData = T)
rld <- rlog(dds)
rld_period <- plotPCA(rld, intgroup=c("condition","period"),returnData = T)
## 繪圖
library(ggrepel)
p <- ggplot(rld_period, aes(x=PC1,y=PC2,color=period,shape=condition))+
geom_point(alpha=.7, size=7)+
geom_text_repel(aes(x=PC1,y=PC2,label=rownames(rld_period),size=14,col="black"))+
#scale_color_manual(values = rev(c("#0072B2","#D55E00","#CC79A7")))+
theme_bw()+
labs(x=paste("PCA 1", sep=""),
y=paste("PCA 2", sep="")) +
scale_colour_discrete(
name="period",
breaks = c("DPI0","DPI3"),
labels = c("DPI 0","DPI 3"))+
# 圖例
theme(panel.background = element_blank(),
panel.grid = element_blank(),
plot.title = element_text(hjust = 0.5),
legend.position="right",
legend.title = element_text(size=14),
legend.text = element_text(size=14),
axis.text.x = element_text(size=14),
axis.text.y = element_text(size=14),
axis.title.x = element_text(size=14),
axis.title.y = element_text(size=14))
#添加標(biāo)簽辕漂,給-log10(padj)>10的基因加上Gene名字
#geom_text_repel(data=subset(dataset, -log10(pval)>10),aes(label=id),col="black", alpha =2)
p
三呢灶、PCoA分析及R包繪圖(3D散點(diǎn)圖)
除轉(zhuǎn)錄組研究以外,在16S微生物的研究中我們會(huì)根據(jù)物種豐度的文件對(duì)數(shù)據(jù)進(jìn)行PCA或者PCoA分析钉嘹,也是我們所說(shuō)的β多樣性分析鸯乃。根據(jù)PCA或者PCoA的結(jié)果看感染組和對(duì)照組能否分開(kāi),以了解微生物組的總體變化情況跋涣。
具體內(nèi)容及繪圖方法可參考下面這篇文章缨睡。
16s—β多樣性分析(R畫(huà)三維PCoA圖)