-
Notifications
You must be signed in to change notification settings - Fork 6
/
plot_PCA（含3D）
84 lines (66 loc) · 3.52 KB
/
plot_PCA（含3D）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
rm(list=ls())
library(ggplot2)
library(scales)
a = read.table("PCA_157_cattle_snp_geno01_maf005-ld502502_nchr.gcta.out.eigenvec",header=F)
m = as.matrix(read.table("PCA_157_cattle_snp_geno01_maf005-ld502502_nchr.gcta.out.eigenval",header=F))
explainm = m/sum(m)
pc1 = paste("PC1","(",percent(explainm[1,], accuracy=0.01),")", sep="")
pc2 = paste("PC2","(",percent(explainm[2,], accuracy=0.01),")", sep="")
pc3 = paste("PC3","(",percent(explainm[3,], accuracy=0.01),")", sep="")
Breed=a[,1]; PC1=a[,3]; PC2=a[,4]; PC3=a[,5]
a$V1 = factor(a$V1,levels=c("ANG","HOL","SIM","SHO","CHR","HAW","KAZ","MON","YB","KUC","MIS","WAG","WL","BRA"))
m_color = c("ANG"="#6BB93F","HOL"="#6BB93F","SIM"="#6BB93F","SHO"="#6BB93F","CHR"="#6BB93F","HAW"="#873186",
"KAZ"="#873186","MON"="#873181","YB"="#873181","MIS"="#18A2CA","KUC"="#18A2CA","WAG"="#E20593",
"WL"="#DBB71D","BRA"="#b35107")
m_shape = c("ANG"=21,"HOL"=21,"SIM"=21,"SHO"=21,"CHR"=21,"HAW"=21,"KAZ"=21,"MON"=21,"YB"=21,"MIS"=25,"KUC"=24,
"WAG"=22,"WL"=23,"BRA"=23)
p = ggplot(data=a, aes(x = PC1, y = PC2, group = Breed,
shape=Breed,color="black",fill=Breed))+
geom_point(size=5,alpha=0.9,stroke=0.2,color="black") +
scale_color_manual(values=m_color)+
scale_fill_manual(values=m_color)+
scale_shape_manual(values=m_shape)
p + labs(x = pc1, y = pc2)+
geom_hline(yintercept=0,linetype="dashed",color="black")+
geom_vline(xintercept=0,linetype="dashed",color="black")+
theme_classic()+
theme(panel.border = element_rect(fill=NA,color="black",linewidth=0.5,linetype ="solid"))+
guides(color = guide_legend(override.aes = list(size=4, stroke=3)))+
theme(legend.title = element_blank())
###############pca-3D ############
library(scatterplot3d)
## 在做完pca后的数据中将一二列删除，加上对应的品种一列，并在第一行中加上行名
dat = read.table("PCA3D.txt",header=T, stringsAsFactors = T )
# 颜色
####, "blue", "red1", "cyan", "blue",
### "magenta4", "yellowgreen", "darkorange3", "grey60", "black","red4", "lawngreen"
##, "black","red4","orange","red","magenta4", "yellowgreen", "darkorange3","purple3","paleturquoise3","green3"
color = c("cyan","darkred","green","dodgerblue1","blue","grey60")
color <- color[as.numeric(dat$Breed)]
shapes = c(1:6)
shapes <- shapes[as.numeric(dat$Breed)]
scatterplot3d(dat[,1:3], main='PCA', type='p',
highlight.3d=F, angle=45, grid=T, box=T, scale.y=1,
cex.symbols=0.9 , col.grid='lightblue',
xlab = "PC1",ylab = "PC2",zlab = "PC3",
pch= shapes,
color = color)
legend("topleft", legend = levels(dat$Breed) ,
col = c("cyan","darkred","green","dodgerblue1","blue","grey60"),
cex = 0.8, xpd = TRUE,inset = 0.03,ncol = 2, bty = "n",
pch= c(1:6))
pdf("pca-3D.pdf", onefile = TRUE, width = 8, height = 8)
diffangle = function(ang){
scatterplot3d(dat[,1:3], main='PCA', type='p',
highlight.3d=F, angle=ang, grid=T, box=T, scale.y=1,
cex.symbols=0.8 , col.grid='lightblue',
xlab = "PC1",ylab = "PC2",zlab = "PC3",
pch= shapes,
color = color)
legend("topright", legend = levels(dat$Breed) ,
col = c("cyan","darkred","green","dodgerblue1","blue","grey60", "black","red4","red","magenta4", "yellowgreen", "darkorange3", "grey60"),
cex = 0.8, xpd = TRUE,inset = 0.03,ncol = 2, bty = "n",
pch= c(1:13))
}
sapply(seq(-360,260,5),diffangle)
dev.off()