用R做主成分分析(PCA)及基本参数的选择

科学计算与大数据

浏览数:161

2019-8-29

AD:资源代下载服务

原理我们已经在前文中讨论过了,这次主要是代码实战

1. 定义一个数据集

data.matrix <- matrix(nrow=100, ncol=10) #横行为基因,纵行为样本
colnames(data.matrix) <- c(
  paste("wt", 1:5, sep=""),
  paste("ko", 1:5, sep=""))
rownames(data.matrix) <- paste("gene", 1:100, sep="")
for (i in 1:100) {
  wt.values <- rpois(5, lambda=sample(x=10:1000, size=1))    ##用柏松分布填充数值
  ko.values <- rpois(5, lambda=sample(x=10:1000, size=1))
  data.matrix[i,] <- c(wt.values, ko.values)
}
head(data.matrix)
dim(data.matrix)

2. PCA主函数

pca <- prcomp(t(data.matrix), scale=TRUE)  ###prcomp函数的横行必须是样本,所以倒置一下

3. 选出值最高的两个PC

我们在PCA原理部分已经讲过了,十个样本,会产生10个PC,我们要选值最高的两个作图
pca.var <- pca$sdev^2  ## sdev是标准偏差,十个样本,就有十个标准偏差,平方是避免负数的干扰
pca.var.per <- round(pca.var/sum(pca.var)*100, 1)  ##求每个样本的variation
barplot(pca.var.per, main="Scree Plot", xlab="Principal Component", ylab="Percent Variation")  ##用柱状图可视化

PC1和PC2的variation比较大,可以拿来做图。故提取出PC1和PC2

pca.data <- data.frame(Sample=rownames(pca$x),
                       X=pca$x[,1],
                       Y=pca$x[,2])
pca.data

4 用ggplot2画PCA

ggplot(data=data_pca,aes(x=PC1,y=PC2,color=group,shape=group))+
  geom_point(size=3)+
  theme_bw()+theme(panel.grid=element_blank())+
  xlab(paste("PC1(",pca.var.per[1],"%","variance)",sep=""))+
  ylab(paste("PC2(",pca.var.per[2],"%","variance)",sep=""))+
  theme(legend.position = c(0.80,0.9))+lines(transM+mu)

好丑,等我看到更好看的包会把这篇再更新一下的~

===============我是分割线=============11月30日更新========
一直说要更新却懒了。画PCA有专用的包,他叫ggord,再配上Y叔的YYplot画置信椭圆。
因为我已经厌倦了调用别人包的生活,用别人写好的东西,永远是技工。
学生信路上偶尔会迷茫,觉得自己像个搬运工。
就酱,仰望星空,脚踏实地。

作者:PriscillaBai