remove HIV samples from pcoa

figure7.R

@@ -16,26 +16,60 @@ require(ellipse)
 source("src/pcoa.functions.R")
 
 # load distance matrix
-dist <- read.table("data/OTU_table/bray_curtis_otu_table.txt")
+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), mapping$V1),]$V1
-mapping.types.ordered <- mapping[match(rownames(dist), mapping$V1),]$V3
+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))
+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)
 
-c <- drawPcoa(obs.pcoa, df, sample.names= mapping.types.ordered, colorvalues=colorvalues)
+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)
 
 
 pdf("results/figure7/figure_7.pdf", width=6, height=4)
-c
+p
 dev.off()
 
+# add label
+p <- p + geom_label(size = 2, aes(fill = mapping.types.ordered), colour = "white", fontface = "bold")
+
+pdf("results/figure7/figure_7_label.pdf", width=6, height=4)
+p
+dev.off()
+
+

figure8_a.R

@@ -71,7 +71,7 @@ colorvalues <- c("green", "red", "blue")
 
 # plot
 g <- ggplot(df, aes(x=as.numeric(as.matrix(CAP1)),
-                    y=as.numeric(as.matrix(CAP2))), color=type.names)
+                    y=as.numeric(as.matrix(CAP2)), label=sample), color=type.names)
 g <- g + geom_point(size= 4, alpha = 0.8, aes(color= type.names))
 g <- g + scale_color_manual(values = colorvalues)
 g <- g + theme_bw() + geom_hline(yintercept = 0,linetype = 3)
@@ -85,9 +85,15 @@ g <- g + ylab(paste("Constrained PCoA 2 (", format(cap.var[2]*100, digits=4), "
 g <- g + theme(plot.title = element_text(size = 8),
                axis.title.x = element_text(size = 8),
                axis.title.y = element_text(size = 8))
+#g <- g + geom_label(size = 2, aes(fill = type.names), colour = "white", fontface = "bold")
 g
 
-
 pdf("results/figure8/figure_8_constrained_by_cohort.pdf", width=6, height=4)
 g
 dev.off()
+
+# add label
+g <- g + geom_label(size = 2, aes(fill = type.names), colour = "white", fontface = "bold")
+pdf("results/figure8/figure_8_constrained_by_cohort_label.pdf", width=6, height=4)
+g
+dev.off()

figure8_b.R

@@ -74,7 +74,7 @@ colorvalues <- c("green", "red", "blue","grey")
 
 # plot
 g <- ggplot(df, aes(x=as.numeric(as.matrix(CAP1)),
-                    y=as.numeric(as.matrix(CAP2))), color=type.names)
+                    y=as.numeric(as.matrix(CAP2)), label=sample), color=type.names)
 g <- g + geom_point(size= 4, alpha = 0.8, aes(color= type.names))
 g <- g + scale_color_manual(values = colorvalues)
 g <- g + theme_bw() + geom_hline(yintercept = 0,linetype = 3)
@@ -90,7 +90,13 @@ g <- g + theme(plot.title = element_text(size = 8),
                axis.title.y = element_text(size = 8))
 g
 
-
 pdf("results/figure8/figure_8_constrained_by_liver.pdf", width=6, height=4)
 g
 dev.off()
+
+# add label
+g <- g + geom_label(size = 2, aes(fill = type.names), colour = "white", fontface = "bold")
+pdf("results/figure8/figure_8_constrained_by_liver_label.pdf", width=6, height=4)
+g
+dev.off()
+

src/pcoa.functions.R

@@ -31,10 +31,10 @@ pcoaProp <- function(obj){
 }
 
 # draw pcoa using ggplot
-drawPcoa <- function(obs.pcoa, df, sample.names=sample.type, colorvalues=colorvalues)
+drawPcoa <- function(obs.pcoa, df, sample.names=sample.type, colorvalues=colorvalues, label=F)
 {
   p <- ggplot(df, aes(x=MDS1,y=MDS2))
-  p <- p + geom_point(aes(size=4, alpha=0.8, colour= sample.names)) # draw points
+  p <- p + geom_point(aes(size=4, alpha=0.8, colour= sample.names, label =sample.names)) # draw points
   #  if (ellipse){
   #   p <- p + geom_segment(aes(x=x.centroid, y=y.centroid, xend=x, yend=y), color="grey80")
   #  p <- p + geom_path(data=conf.rgn) + theme_bw() + geom_hline(yintercept=0,linetype=3)
@@ -54,6 +54,9 @@ drawPcoa <- function(obs.pcoa, df, sample.names=sample.type, colorvalues=colorva
   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)
+  if (label){
+    p <- p + geom_label(size = 2, aes(fill = type.names), colour = "white", fontface = "bold")
+  }
   return(p)
 }