三维的PCA分析-让您的文章C位出道
百迈客推出年末活动促销,发文有礼、推广有礼以及多种产品钜惠来袭,百迈客为您倾情打造科研福利,您还在等什么?快快行动起来领取您的超级奖励吧!(详情请见“决战2020!品类全,力度大,仅此一次!”)
说起PCA主成分分析,大家都知道这是一种很常见的降维方法,可以做到将多维的数据展示在尽量低的维度上,并且保留尽可能多的原数据的特性。PCA结果的最常见展现方式便是分别以第一主成分和第二主成分为横纵坐标绘制的PCA散点图。当第一和第二主成分的解释度足够高时,二维的PCA图就已经能够很好的展示数据的真实情况,但如果前两个主成分的解释度不够高,二维的PCA图就会损失较多的信息,那么这个时候就可以尝试绘制三维PCA图。例如,对于我们熟悉的R中的内置鸢尾花数据集来讲,PCA分析后前两个主成分的解释度之和高达95%,这种数据简单绘制二维PCA即可。
library(ggplot2)#鸢尾花各列的数据差别较大,为了使不同维度的方差可以放在一起比较,这里先将数据进行了标准化
data <- scale(iris[,1:4],center = T,scale = T)
#PCA分析
pca <- prcomp(data)
#提取坐标值及解释度
pca_result <- pca$x
pcatab <- data.frame(data.frame(pca_result[,1:3]))
pca_summary <- summary(pca)
importance <- pca_summary$importance[2,]
xlab <- paste0("PC1: ",round(as.numeric(importance[1])*100,2),"% variance")
ylab <- paste0("PC2: ",round(as.numeric(importance[2])*100,2),"% variance")
#使用ggplot绘制二维PCA图
ggplot(pcatab,aes(x=PC1,y=PC2,color=iris$Species))+geom_point(size=3)+
#去除背景
theme_bw()+theme(panel.grid = element_blank())+
#设置标题
labs(x=xlab,y=ylab,color="Species",title = "PCA Plot")+
theme(plot.title = element_text(face = "bold",size = 20,hjust = 0.5))+
#添加0刻度虚线
geom_vline(xintercept=0,linetype=4,color="grey")+
geom_hline(yintercept=0,linetype=4,color="grey")+
#拟合分组的椭圆
stat_ellipse(type = "t")+coord_equal
(向左滑动查看更多)#导入数据(物种丰度表和分组信息表均可在咱们百迈客的分析结果中找到哦~)
data <- read.table("treat1.phylum.reabundance",header = T,check.names = F,row.names = 1)
grouptab <- read.table("treat1.group",header = T,check.names = F)
colnames(grouptab) <- c("sample","group")
#调整物种丰度表的样本顺序与分组信息表保持一致
data <- data[,grouptab$sample]
#PCA分析
pca <- prcomp(t(data))
#提取坐标值及解释度
pca_result <- pca$x
pcatab <- data.frame(data.frame(pca_result[,1:3]))
pca_summary <- summary(pca)
importance <- pca_summary$importance[2,]
xlab <- paste0("PC1: ",round(as.numeric(importance[1])*100,2),"% variance")
ylab <- paste0("PC2: ",round(as.numeric(importance[2])*100,2),"% variance")
#使用ggplot绘制二维PCA图
ggplot(pcatab,aes(x=PC1,y=PC2,color=grouptab$group))+geom_point(size=5)+
#去除背景
theme_bw()+theme(panel.grid = element_blank())+
#设置标题
labs(x=xlab,y=ylab,color="Species",title = "PCA Plot")+
theme(plot.title = ele ment_text(face = "bold",size = 20,hjust = 0.5))+
#添加0刻度虚线
geom_vline(xintercept=0,linetype=4,color="grey")+
geom_hline(yintercept=0,linetype=4,color="grey")+
#拟合分组的椭圆
stat_ellipse(type = "t")+coord_equal()
#加载scatterplot3d包,没有的可以通过install.packages("scatterplot3d")命令进行安装library(scatterplot3d)
zlab <- paste0("PC3: ",round(as.numeric(importance[3])*100,2),"% variance")
#绘制三维PCA图,
scatterplot3d(pcatab$PC1,pcatab$PC2,pcatab$PC3,type = "h",color = rep(col,freq),pch = 16,cex.symbols = 3,xlab = xlab,ylab = ylab,zlab = zlab,cex.axis = 1,cex.lab = 1.2,lwd = 4,font.axis = 2,font.lab = 2)
#添加图例library(RColorBrewer)
group <- unique(grouptab$group)
freq <- as.data.frame(table(as.character(grouptab$group)))[,2]
col <- brewer.pal(length(group),"Set1")
legend(x = "topleft",legend=group,fill = col,bg = "white",box.lwd = 2,text.font = 2)
(向左滑动查看更多)
scatterplot3d(pcatab$PC1,pcatab$PC2,pcatab$PC3,type = "h",color = rep(col,freq),pch = 16,cex.symbols = 3,xlab = xlab,ylab = ylab,zlab = zlab,cex.axis = 1,cex.lab = 1.2,lwd = 4,font.axis = 2,font.lab = 2,lty.axis = 0,lty.hide = 1)
legend(x = "topleft",legend=group,fill = col,bg = "white",box.lwd = 2,text.font = 2)
(向左滑动查看更多)
就这样,一个基础的三维PCA图就绘制完成啦,可以看到三维主成分的解释度高达87.65%,高度还原了数据的真实性、丰富性,可见三维比二维不仅仅只是多了一维!心动不如行动起来,百迈客专业的生信团队为您提供优质的PCA三维个性化分析以及其他多种高级别的分析图谱,让您的文章快速C位出道,让您的文章PCA比别人多一维,永远快人一步。如果想要让您的文章与众不同,快快行动起来,年末活动促销各种大奖等您来拿,欢迎各位老师致电锤询!
文:小师
排版:市场部
往期回顾:
单细胞RNA-seq揭开了先天性纤维化肺部中异位和畸变常驻细胞群的神秘面纱
Science Advances | 使用空间转录组构建小鼠的大脑图谱
ATAC-Seq、RNA-Seq和WGS多组学联合揭示PRDM10重组肿瘤发病机制
第七届全国功能基因组学高峰论坛
兑奖说明:
1、支付宝88.88元红包抽取十人已经全部转账完毕!
2、次账号中奖的10人中,已经兑奖七人,请下面三位将姓名、电话、单位、百迈客云账号发送至:zhangyx@biomarker.com.cn
Bio-bd-zgq、天雪、s56458572
3、培训班代金券中奖12人,已经兑奖4人,请下面八位将姓名、电话、单位、百迈客云账号发送至:zhangyx@biomarker.com.cn
芹泽楠,slowness(两人均获得800元代金券)。
ZengHua,方欣,阮阮,对月聆雪弱,牛哞哞-CY,起风了(六人均获得500元代金券)
4,在线时长奖品(U盘+纪念帆布袋):下图为直播在线超过4小时的人(每个人的ID是固定的,已经根据ID合并所有同类项分钟),可以点击放大查看大图,请在上面的老师识别下面二维码填写信息!收集超过一半人数后,我们会立即安排邮寄!
5,下图为直播在线时长为230分钟-239分钟的,感谢各位的支持,我们也有小礼品送上!收集超过一半人数后,我们会立即安排邮寄!
请中奖的小伙伴们识别二维码填写信息!