PCA主成分分析

在这里我不去介绍原理，着重讲讲如何做PCA分析，需要什么样的数据
1、数据格式如下，很简单就是一个表达谱数据，行名为基因，列名为样本，本次例子中总共有三组信息，有多少组，由你的数据决定

#读取表达谱数据

data=read.table("input.txt",header=T,sep="\t",row.names=1)  
#转至
data=t(as.matrix(data)) 

#将数据的行名单独赋值

data.class <- rownames(data)

#pca分析

data.pca <- prcomp(data, scale. = TRUE) 

#输出特征向量，即基因的PC

write.table(data.pca$rotation,file="PC.xls",quote=F,sep="\t")   

#输出新表，每个样本的PC值

write.table(predict(data.pca),file="newTab.xls",quote=F,sep="\t")  

pca.sum=summary(data.pca)
 #输出PC比重
write.table(pca.sum$importance,file="importance.xls",quote=F,sep="\t")  

pdf(file="pcaBarplot.pdf",width=15)   #柱状图
barplot(pca.sum$importance[2,]*100,xlab="PC",ylab="percent",col="skyblue")
dev.off()

#碎石图

pdf(file="pcaPlot.pdf",width=15)  
plot(pca.sum$importance[2,]*100,type="o",col="red",xlab="PC",ylab="percent")
dev.off()

#pca 2d plot

library(ggplot2)
#group分组信息，需要大家自己输入
group=c(rep("con",5),rep("A",5),rep("B",3))  
pcaPredict=predict(data.pca)
PCA = data.frame(PCA1 = pcaPredict[,1], PCA2 = pcaPredict[,2],group=group)
PCA.mean=aggregate(PCA[,1:2],list(group=PCA$group),mean)
#定义椭圆
veganCovEllipse<-function (cov, center = c(0, 0), scale = 1, npoints = 100) 
  {
    theta <- (0:npoints) * 2 * pi/npoints
    Circle <- cbind(cos(theta), sin(theta))
    t(center + scale * t(Circle %*% chol(cov)))
  }
df_ell <- data.frame()
for(g in levels(PCA$group)){
  df_ell <- rbind(df_ell, cbind(as.data.frame(with(PCA[PCA$group==g,],
                  veganCovEllipse(cov.wt(cbind(PCA1,PCA2),
                  wt=rep(1/length(PCA1),length(PCA1)))$cov,
                  center=c(mean(PCA1),mean(PCA2))))),group=g))
}

pdf(file="PCA2d.pdf")
ggplot(data = PCA, aes(PCA1, PCA2)) + geom_point(aes(color = group)) +
    geom_path(data=df_ell, aes(x=PCA1, y=PCA2,colour=group), size=1, linetype=2)+
    annotate("text",x=PCA.mean$PCA1,y=PCA.mean$PCA2,label=PCA.mean$group)+
    theme_bw()+
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
dev.off()

#pca 3d plot
library(pca3d)
library(rgl)
pca3d(data.pca, components = 1:3, group = group, show.centroids=TRUE, show.group.labels =TRUE) 

分析很简单，难的是如何解释PCA。
以我目前的理解，我认为PCA分析主要有两方面的作用，一是查看样本的分布情况，将离群的样本剔除，以免影响实验结果。二是筛选基因，通过一定的筛选将PC分数较高的基因挑选出来，有过滤数据的作用。后面可以对过滤后的数据进行lasso回归分析。