PCA画图学习

昨天,想用R语言画个PCA图,来确认一批样本重复测序的吻合度,虽然已经单个数据审视和检查过了,还是用PCA看一下整体的情况,于是决定学画PCA图。使用搜狗微信搜索,找到一篇百迈客微信公众号的视频讲解,讲的不错,跟着以1.5倍速听完就可以画出自己的图了,挺开心的。
把这个公众号的视频地址放在这:https://v.qq.com/x/page/p0384ninmc2.html
话不多说,开始画图:

1.使用基本给力函数来画

当然,这样的图画出来有一点点丑,不过依然能说明问题,经过修饰应该也能变漂亮。

library(ape)
library(ggplot2)
#download.file("https://github.com/waasp/plane.R/raw/master/plane.R",cwd)
otu_table<-read.table('/Users/zd200572/Desktop/v4-repeat/2otu_table3_L6.txt', row.names=1, header = TRUE, sep = '\t', check.names = FALSE)
dim(otu_table)
pca <- prcomp(t(otu_table[,-1]))
head(pca$x) #这是用来画图的数据
plot(pca$x[,1:2])
group <- factor(c(rep('first', 10),rep('second', 10)))
colour <- rainbow(length(unique(group)))
colour <- colour[as.numeric(factor(group))]
plot(pca$x[,1:2], col = colour, pch = c(21,22)[group])
legend('topleft', legend = levels(group), col = unique(colour), pch = c(21,22))
title('两次测序结果吻合情况图')

画出来的图便是这个样子的:

ggplot2大法

ggplot2的名字真是声震四方,画出来的图自然也是漂亮的没话说,看代码:

library(ggplot2)
otu_table<-read.table('/Users/zd200572/Desktop/v4-repeat/3otu_table3_L6.txt', row.names=1, header = TRUE, sep = '\t', check.names = FALSE)
#group2 <- factor(c(rep('first', 10),rep('second', 10)))c(1,1,2,2,1,2,1,1,2,2,1,2,2,1,1,1,2)c(1,1,1,2,1,2,1,1,2,1,2,2,2,1,1,2,2,1,2,2)
group2 <- factor(c(2,1,2,1,2,2,1,1,1,2,2,1,1,1,2,2)) #这是我的分组信息,样本不多,就这样加上去了
#length(group2)
pca <- prcomp(t(otu_table))
pca_result <- as.data.frame(pca$x)
dim(pca_result)
pca_result <- cbind(pca_result, group2)
p <- ggplot(pca_result) + geom_point(aes(x=pca_result[,1], y=pca_result[,2], color = pca_result$group2, shape = pca_result$group2, size =1 ))
p <- p + theme(legend.title = element_blank()) + labs(x = "PCA1", y="PCA2")
p <- p + theme_bw()
p

图画出来就是不一样:

吻合度还是不错的,基本重合,说明方法还是比较可靠的。

3D PCA分析

平面图只是两个主成分进行分析,3D图可以画出3个主成分,一个图顶3个。

library(scatterplot3d)
par(mar = c(5.1, 2.1, 4.1, 8.1), xpd=T)
length(colour)
scatterplot3d(pca_result[, 1:3],pch = 20, color = pca_result$group2, angle = 60, main = "test_3d",cex.symbols = 2)
legend("topright", legend = group[c(1,11)], col = colour,pch=20)


虽然还是需要修饰才能漂亮点,但有没有感受到一点厉害,R语言给力呀!

另外一个画图方式

这种方式可以圈起来同样分组的数据,数据量大时展示会更清晰,想起了美国肠道计划那个图。

library(ggfortify)
data2 <- cbind(t(otu_table[,-1]), group2)
autoplot(pca,data = data2, colour = 'group2', label = T, frame = F)

发表评论