-
Notifications
You must be signed in to change notification settings - Fork 4
/
figure_2a.R
78 lines (54 loc) · 2.42 KB
/
figure_2a.R
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
# code adapted from http://enterotype.embl.de/enterotypes.html
# cleanup
rm(list=ls())
#source("https://bioconductor.org/biocLite.R")
#biocLite("phyloseq")
require(phyloseq)
require(ape)
require(cluster)
require(clusterSim)
require(ade4)
require(ggplot2)
require(RColorBrewer)
source("src/pcoa.functions.R")
# load distance matrix
dist_mat <- read.table("data/OTU_table/bray_curtis_otu_table.txt")
# read mapping
mapping <- read.table("data/mapping_cat.txt", sep=";", header=F)
# get sample ID of AIH/FAP
idstoremove <- mapping[which(mapping$V2 == "HIV" |
mapping$V2 == "FAP" ),]
# remove HIV, FAP
dist_mat2 <- dist_mat[-match(idstoremove$V1, colnames(dist_mat))
,-match(idstoremove$V1, colnames(dist_mat))]
# get sample types from mapping file
mapping.names.ordered <- mapping[match(rownames(dist_mat2), mapping$V1),]$V1
mapping.types.ordered <- mapping[match(rownames(dist_mat2), mapping$V1),]$V3
obs.pcoa <- pcoa(as.dist(dist_mat2))
types.ordered <- mapping[match(rownames(obs.pcoa$vectors), mapping$V1),]$V3
colorvalues <- c("green", "red", "blue")
df <- data.frame(MDS1=obs.pcoa$vectors[,1], MDS2=obs.pcoa$vectors[,2], genotype=types.ordered)
#c <- drawPcoa(obs.pcoa, df, sample.names= mapping.types.ordered, colorvalues=colorvalues, label=F)
p <- ggplot(df, aes(x=MDS1,y=MDS2, label=mapping.names.ordered))
p <- p + geom_point(aes(size=4, alpha=0.8, colour= mapping.types.ordered)) # draw points
p <- p + xlab(paste("PCo 1 (", format(pcoaProp(obs.pcoa)[1]*100, digits=4), "%)", sep=""))
p <- p + ylab(paste("PCo 2 (", format(pcoaProp(obs.pcoa)[2]*100, digits=4),"%)", sep=""))
p <- p + theme(legend.position="none")
p <- p + theme(plot.title = element_text(size = 8),
axis.title.x = element_text(size = 8),
axis.title.y = element_text(size = 8))
p <- p + theme(plot.title = element_text(size = 8),
axis.title.x = element_text(size = 8),
axis.title.y = element_text(size = 8))
p <- p + scale_color_manual(values = colorvalues)
p <- p + theme_bw() + geom_hline(yintercept = 0,linetype = 3)
p <- p + geom_vline(xintercept = 0,linetype = 3)
p <- p + theme(legend.position="none") + coord_equal()
pdf("results/figure_2_print.pdf", width=6, height=4)
p
dev.off()
# add label
p <- p + geom_label(size = 2, aes(fill = mapping.types.ordered), colour = "white", fontface = "bold")
pdf("results/figure_2_label.pdf", width=6, height=4)
p
dev.off()