200828_R_notebook_git.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## R notebook for producing figures as in \"Minnesota peat viromes reveal terrestrial and aquatic niche partitioning for local and global viral populations\"\n",
    "Anneliek ter Horst\n",
    "\n",
    "- Accumulation plots\n",
    "- Venn diagrams\n",
    "- PCOA plots\n",
    "- Worldmap chart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load libraries and set working directory\n",
    "setwd(\"/Users/Anneliek/R_scripts/\")\n",
    "\n",
    "library(vegan)\n",
    "library(labdsv)\n",
    "library(ggplot2)\n",
    "library(plyr)\n",
    "library(phyloseq)\n",
    "library(pheatmap)  \n",
    "library(tidyr)\n",
    "library(RColorBrewer)\n",
    "library(scales)\n",
    "library(Ecfun)\n",
    "library(waffle)\n",
    "library(grid)\n",
    "library(broom)\n",
    "library(tidyverse)\n",
    "library(ape)\n",
    "library(dplyr)\n",
    "library(ggplot2)\n",
    "library(ggrepel)\n",
    "library(forcats)\n",
    "library(scales)\n",
    "library(ggpmisc)\n",
    "library(maps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accumulation plots\n",
    "- in: viral OTU table from viromes\n",
    "- viral OTU table from metagenomes from transect\n",
    "- viral OTU table from metagenomes from chambers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in the virome OTU data (this table contains integers instead of values with a comma)\n",
    "votu1 <- read.table(\"../Data/virome_ints.txt\",\n",
    "                   sep=\"\\t\", header=T, row.names=1)\n",
    "\n",
    "# transpose the data\n",
    "votu <- t(votu1)\n",
    "\n",
    "# Load in the bulk OTU data \n",
    "otu1 <- read.table(\"../Data/bulk_ints_own.txt\", \n",
    "                  sep=\"\\t\", header=T, row.names=1)\n",
    "\n",
    "# transpose the data\n",
    "otu <- t(otu1)\n",
    "\n",
    "\n",
    "# Load in the bulk OTU data (bulk metaGs JGI)\n",
    "all_otu1 <- read.table(\"../Data/bulk_ints_82.txt\", \n",
    "                  sep=\"\\t\", header=T, row.names=1)\n",
    "\n",
    "# transpose the data\n",
    "all_otu <- t(all_otu1)\n",
    "\n",
    "\n",
    "# Calculate accumulaton curves using the random method with 100 permutations\n",
    "# for vOTUs\n",
    "sp <- specaccum(votu, method = \"random\", permutations = 100)\n",
    "\n",
    "# for 5 metaG\n",
    "sp.otu <- specaccum(otu, method = \"random\", permutations = 100)\n",
    "\n",
    "# for 82 metaG\n",
    "sp.all_otu <- specaccum(all_otu, method = \"random\", permutations = 100)\n",
    "\n",
    "# Extract each of the permutations and reformat\n",
    "perm <- sp$perm\n",
    "\n",
    "perm.tidy <- as.tibble(perm) %>% \n",
    "  mutate(Sites = 1:nrow(.)) %>% \n",
    "  gather(key = \"Permutation\", value = \"Species\", -Sites) \n",
    "\n",
    "# Extract the consensus richness per number of sites sampled and reformat\n",
    "richness <- data.frame(Sites = sp$sites, Species = sp$richness)\n",
    "\n",
    "# Extract each of the permutations and reformat\n",
    "perm.otu <- sp.otu$perm\n",
    "\n",
    "perm.otu.tidy <- as.tibble(perm.otu) %>% \n",
    "  mutate(Sites = 1:nrow(.)) %>% \n",
    "  gather(key = \"Permutation\", value = \"Species\", -Sites) \n",
    "\n",
    "# Extract the consensus richness per number of sites sampled and reformat\n",
    "richness.otu <- data.frame(Sites = sp.otu$sites, Species = sp.otu$richness)\n",
    "\n",
    "# Extract each of the permutations and reformat\n",
    "perm.all_otu <- sp.all_otu$perm\n",
    "\n",
    "perm.all_otu.tidy <- as.tibble(perm.all_otu) %>% \n",
    "  mutate(Sites = 1:nrow(.)) %>% \n",
    "  gather(key = \"Permutation\", value = \"Species\", -Sites) \n",
    "\n",
    "# Extract the consensus richness per number of sites sampled and reformat\n",
    "richness.all_otu <- data.frame(Sites = sp.all_otu$sites, Species = sp.all_otu$richness)\n",
    "\n",
    "all.perm  <- rbind(mutate(perm.tidy, Extraction = \"Viromes\"),\n",
    "                  mutate(perm.otu.tidy, Extraction = \"Bulk_own\"),\n",
    "                   mutate(perm.all_otu.tidy, Extraction = \"Bulk_all\"))\n",
    "\n",
    "all.richness  <- rbind(mutate(richness, Extraction = \"Viromes\"),\n",
    "                     mutate(richness.otu, Extraction = \"Bulk_own\"),\n",
    "                      mutate(richness.all_otu, Extraction = \"Bulk_all\"))\n",
    "\n",
    "# create the actual plot\n",
    "pdf('accumplot_all_permute.pdf')\n",
    "ggplot(all.perm, aes(Sites, Species, color = Extraction)) +\n",
    "  geom_point(alpha = 0.5, shape = 1, size = 2.5) +\n",
    "  #scale_color_brewer(palette = my_palette) +\n",
    "  scale_colour_manual(values = palet) +\n",
    "  geom_line(data = all.richness, size = 1) +\n",
    "  theme_bw() +\n",
    "  theme(text = element_text(size = 18))\n",
    "dev.off()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statistical tests\n",
    "- Mantel using pearson method\n",
    "- PERMANOVA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# mantel test with pearson method.\n",
    "mantel(otu.dist, env.dist, permutations=99999, method=\"pearson\", na.rm = T)\n",
    "\n",
    "\n",
    "# run permanova for continious data\n",
    "pmanova <- adonis(as.dist(otu.dist) ~ depth, data = cat_meta.ordered)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Venn diagrams\n",
    "- in: csv file with 1 and 0 for each of the clusters and the categories present\n",
    "- 0 means that vOTUs of that category arent present, 1 means they are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load in csv with 0 and 1 data \n",
    "data_soil_water <- read.delim('../Data/200402_soil_water_ncbi.csv',header=T, sep=',')\n",
    "\n",
    "# define the sizes of each of the circles in the venn\n",
    "area1 = nrow(subset(data_soil_water, marine == 1))\n",
    "area2 = nrow(subset(data_soil_water, soil == 1))\n",
    "area3 = nrow(subset(data_soil_water, freshwater == 1))\n",
    "n12 = nrow(subset(data_soil_water, marine == 1 & soil ==1))\n",
    "n13 = nrow(subset(data_soil_water, marine == 1 & freshwater ==1))\n",
    "n23 = nrow(subset(data_soil_water, soil == 1 & freshwater ==1))\n",
    "n123 = nrow(subset(data_soil_water, marine == 1 & soil ==1 & freshwater == 1))\n",
    "\n",
    "# draw the venn, hardcoded the number of singleton clusters for each\n",
    "venn.plot <- draw.triple.venn(\n",
    "area1 = area1 + 26018,\n",
    "area2 = area2 + 15661,\n",
    "area3 = area3 + 5396,\n",
    "n12 = n12,\n",
    "n13 = n13,\n",
    "n23 = n23,\n",
    "n123 = n123,\n",
    "category = c(\"marine\", \"soil\", \"freshwater\"),\n",
    "fill = c(\"red\", \"blue\", \"green\"),\n",
    "lty = \"blank\",\n",
    "cex = 2,\n",
    "cat.cex = 2,\n",
    "cat.col = c(\"red\", \"blue\", \"green\"),\n",
    "euler.d = TRUE,\n",
    "scaled = TRUE\n",
    ");\n",
    "\n",
    "ggsave('../figures/200506_3way_marine_freshwater_soil.pdf', venn.plot) \n",
    "dev.off()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boxplots\n",
    "- In: table with the total number of populations from each sample\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data in\n",
    "total_pops <- read.table(\"../Data/total_populations_against_own_db.txt\", sep=\"\\t\", header=T)\n",
    "\n",
    "# print means for both bulk and virome (populations)\n",
    "colMeans(total_pops)\n",
    "\n",
    "#boxplot with a logarithmic scale\n",
    "pdf(\"boxplot.pdf\")\n",
    "boxplot(total_pops, par(cex.axis=1.3), par(cex.lab=1.3), col=c(\"#1c5ec9\",\"#f6b420\",'#fee391', '#df3123'), # box colors\n",
    "        names=c(\"5 bulk\", '2015 bulk', '2016 bulk', \" 5 viromes\"), # box names\n",
    "        ylab=\"viral genomes per sample (log scale)\", # y-axes name\n",
    "        log=\"y\")+  # logaritmic\n",
    "  theme_bw()\n",
    "  dev.off()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCOA plots\n",
    "- In: Viral OTU table\n",
    "- Table with available metadata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the files required and transpose\n",
    "otu <- read.table(\"../Data/190807_JGI_samples_20152016.txt\", \n",
    "              sep=\"\\t\", header=T, row.names=1)\n",
    "\n",
    "# transpose data\n",
    "otu <- t(otu)\n",
    "# remove rows with all zeros\n",
    "otu <- otu[apply(otu[,-1], 1, function(x) !all(x==0)),]\n",
    "\n",
    "# categorical metadata\n",
    "cat_meta <- read.table(\"../Data/191203_all_meta_correct.csv\", \n",
    "                         sep=\",\", header=T, row.names=1)\n",
    "                 \n",
    "# print head of the metadata\n",
    "head(cat_meta)\n",
    "head(otu)\n",
    "                 \n",
    "# transform the data\n",
    "otu.xform <- decostand(otu, method=\"hellinger\")\n",
    "\n",
    "# create dissimilarity matrix for the vOTU table with bray method\n",
    "otu.dist <- as.matrix(vegdist(otu.xform, method='bray'))\n",
    "\n",
    "# nmds\n",
    "otu.nmds <- metaMDS(otu.dist)\n",
    "otu.nmds$stress\n",
    "\n",
    "\n",
    "# R will not automatically bind datapoints with the same name, but randomly bind them\n",
    "# therefore order cat data with this\n",
    "cat_meta.ordered <- cat_meta[match(row.names(otu.nmds$points), row.names(cat_meta)),]   \n",
    "\n",
    "                 \n",
    "# perform pcoa with ape package pcoa\n",
    "pcoa <- pcoa(as.dist(otu.dist))\n",
    "\n",
    "# make a dataframe named axes, put pcoa values in there\n",
    "axes <- as.data.frame(pcoa$vectors)\n",
    "\n",
    "# Give df extra column with the rownames in it \n",
    "axes$SampleID <- rownames(axes)\n",
    "\n",
    "# put the metadata in the same dataframe, with correct sample name\n",
    "cat_meta.ordered$SampleID <- rownames(cat_meta.ordered)\n",
    "\n",
    "\n",
    "# calculate the eigenvalues for each pcoa axes \n",
    "eigval <- round(pcoa$values$Relative_eig * 100, digits = 2)\n",
    "\n",
    "# merge those dfs\n",
    "axes <- merge(cat_meta.ordered, axes, by.x = \"SampleID\", by.y = \"SampleID\")\n",
    "\n",
    "\n",
    "# Put those eigenvalues in a df so they easy to get to. \n",
    "eigval <- data.frame( PC = 1:length(eigval), Eigval = eigval)\n",
    "#head(eigval) # see top eigenvalues\n",
    "eigval[[1,2]] # see first axes percentage\n",
    "eigval[[2,2]] # second axes\n",
    "eigval[[3,2]] # third axes\n",
    "eigval[[4,2]] # fourth axes\n",
    "\n",
    "axes\n",
    "                 \n",
    "# write to pdf\n",
    "pdf(\"depth_year20152016.pdf\")\n",
    "\n",
    "# set plot\n",
    "p <- ggplot(axes, aes(Axis.1, Axis.2), width = 12, height = 8) \n",
    "\n",
    "# set color of the points as the factor depth, shape as year, set size and see-throughness\n",
    "p + geom_point(aes(colour=as.character(depth), shape=as.character(Year)), size = 4,alpha=0.7, stroke=1) +\n",
    "  \n",
    "  # set text for the axis lables\n",
    "  xlab(paste(\"PCo1 (\", eigval$Eigval[1], \" %)\", sep = \"\")) +  # or somthing else\n",
    "  ylab(paste(\"PCo2 (\", eigval$Eigval[2], \" %)\", sep = \"\")) +\n",
    "  \n",
    "  # dotted lines through the plot\n",
    "  #geom_vline(xintercept = 0, linetype = 2) +\n",
    "  #geom_hline(yintercept = 0, linetype = 2) +\n",
    "  \n",
    "  # set the colors of the points with the colorbrewer pallet\n",
    "  scale_color_brewer(name = \"Depth\", palette = \"Dark2\") +\n",
    "\n",
    "  # you can set colors manually by this: \n",
    "  #scale_color_manual(values = c(\"#E78AC3\", \"#A6D854\", \"#FFD92F\", \"#FFD92F\")) +\n",
    "\n",
    "  # set shapes of the points\n",
    "  scale_shape_manual(name = \"Year\", values=c(16,17)) +\n",
    "\n",
    "  # tell where the legend has to be\n",
    "  guides(color = guide_legend(title.position = \"top\", title.hjust = 0.5)) +\n",
    "  theme_bw() +\n",
    "  \n",
    "  # set text size for whole graph. set the background color (white with no lines)\n",
    "  theme(text = element_text(size = 20), panel.grid.major = element_blank(),\n",
    "  panel.grid.minor = element_blank(),\n",
    "        legend.position = \"left\") \n",
    "dev.off()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Worldmap chart\n",
    "- in: Table with for each different place the x and y coordinates, the number of vOTUs found there and the color of the dot you want it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in the datapoints of clusters found at different places\n",
    "VCs <- read.csv(\"Data/200408num_of_virus_on_map.csv\", \n",
    "                    sep=\",\", header=T, row.names=1)\n",
    "# show table\n",
    "head(VCs)\n",
    "\n",
    "# pull the world map\n",
    "world_map = map_data(\"world\")\n",
    "\n",
    "#Creat a base plot with gpplot2\n",
    "p = ggplot() + coord_fixed() +\n",
    "  xlab(\"\") + ylab(\"\") \n",
    "\n",
    "#Add map to base plot\n",
    "base_world_messy = p + geom_polygon(data=world_map, aes(x=long, y=lat, group=group), \n",
    "                                     colour=\"#b8e186\", fill=\"#b8e186\")\n",
    "\n",
    "#Strip the map down of the lat/long lines \n",
    "cleanup = theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), \n",
    "        panel.background = element_rect(fill = '#e7f2fb', colour = '#e7f2fb'), \n",
    "        axis.line = element_line(colour = \"#e7f2fb\"), legend.position=\"#e7f2fb\",\n",
    "        axis.ticks=element_blank(), axis.text.x=element_blank(),\n",
    "        axis.text.y=element_blank())\n",
    "\n",
    "# clean up the messy world\n",
    "base_world = base_world_messy + cleanup\n",
    "\n",
    "# map points to this world\n",
    "map_data = base_world + geom_point(data=VCs, aes(x = lon, y = lat, size = sqrt(tot_viruses), \n",
    "                color=factor(color), alpha = color)) +\n",
    "                scale_color_manual(name = 'color', labels=c(\"SPRUCE\", \"Other soils\"),\n",
    "                                   values = c('#de77ae', '#8e0152')) +\n",
    "                scale_alpha_manual(name = 'color', values= c(0.5, 0.8)) \n",
    "    \n",
    "\n",
    "# write to PDF\n",
    "pdf('200511_viruses_on_map_green.pdf')\n",
    "\n",
    "# print to file   \n",
    "print(map_data)\n",
    "\n",
    "dev.off()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}