library(phyloseq)
library(ggplot2)
theme_set(theme_bw())

setwd("Path/to/PhyloseqObject")
load("DropletPlate_subsampled_filtered.RData", verbose=TRUE)

# median normalization

total <- median(sample_sums(DropletPlate_subsampled_filtered))
standf <- function(x, t=total) round(t * (x / sum(x)))
DropletPlate_subsampled_filtered_normed <- transform_sample_counts(DropletPlate_subsampled_filtered, standf)


# ordination

my.colors<-c('#ff7f00','#66a61e','#cab2d6','#7570b3','#e7298a', '#b15928')
PCoA<-plot_ordination(DropletPlate_subsampled_filtered_normed, ordinate(DropletPlate_subsampled_filtered_normed, "MDS", distance(DropletPlate_subsampled_filtered_normed, method="bray")),color = "sampleType") +
  geom_point(size =1, alpha=0.6) +
  scale_colour_manual(values=my.colors, name="sample type",labels=c("droplet PC", "droplet sample", "plate PC", "plate sample", "soil", "ext. community"))


# community structure

DropletPlate_subsampled_filtered_com<-merge_samples(DropletPlate_subsampled_filtered, "sampleType")
sample_names(DropletPlate_subsampled_filtered_com)<-c("droplet PC", "droplet sample", "plate PC", "plate sample", "soil", "ext. community")

phyl_com<-tax_glom(DropletPlate_subsampled_filtered_com, taxrank="phylum")

phyl_com_top10<-names(sort(taxa_sums(phyl_com), TRUE))[1:10]
phyl_com_top10_prune<-prune_taxa(phyl_com_top10, phyl_com)
phyl_com_top10_prune_rel<-transform_sample_counts(phyl_com_top10_prune, function(x) x / sum(x) )
phyl_com_top10_prune_rel_bar<-psmelt(phyl_com_top10_prune_rel)

my.colors<-c('#e7298a','#d95f02','#66a61e','#1b9e77','#7570b3','#e6ab02',"#f768a1", "#CD9BCD","#49006a", "#00441b")

community<-ggplot(phyl_com_top10_prune_rel_bar, aes(x = reorder(Sample, group), y = Abundance, fill = phylum)) + 
  geom_bar(stat = "identity") +
  scale_fill_manual(values = my.colors) +
  ylab("Relative Abundance") +
  theme(axis.text.x  = element_text(angle=45, vjust=0.95, hjust=0.95, size=10), axis.title.x = element_blank()) 

gPCoA <- ggplotGrob(PCoA)
gcommunity <- ggplotGrob(community)
gPCoA$widths <- gcommunity$widths
gPCoA$heights <- gcommunity$heights


library(cowplot) 
beta<-plot_grid(gPCoA, gcommunity, labels = c("A", "B"), ncol = 2, nrow=1)

ggsave(file="Figure2.svg", plot=beta, height=10, width=21, units="cm", pointsize=12)
