原理我们已经在前文中讨论过了,这次主要是代码实战
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画置信椭圆。
因为我已经厌倦了调用别人包的生活,用别人写好的东西,永远是技工。
学生信路上偶尔会迷茫,觉得自己像个搬运工。
就酱,仰望星空,脚踏实地。