導(dǎo)讀
主成分分析(Principal Components Analysis，PCA)，也稱(chēng)主分量分析或主成分回歸分析法，是一種無(wú)監(jiān)督的數(shù)據(jù)降維方法。PCA通過(guò)線(xiàn)性變換將原始數(shù)據(jù)變換為一組各維度線(xiàn)性無(wú)關(guān)的表示，可用于提取數(shù)據(jù)的主要特征分量，常用于高維數(shù)據(jù)的降維。這種降維的思想首先減少數(shù)據(jù)集的維數(shù)，同時(shí)還保持?jǐn)?shù)據(jù)集的對(duì)方差貢獻(xiàn)最大的特征，最終使數(shù)據(jù)直觀(guān)呈現(xiàn)在二維坐標(biāo)系。
直觀(guān)上，第一主成分軸 優(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圖）
