diff --git a/.Rbuildignore b/.Rbuildignore
index 91114bf..a45eef6 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,2 +1,3 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
+^\.httr-oauth$
diff --git a/.gitignore b/.gitignore
index 807ea25..f427a7f 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
 .Rproj.user
 .Rhistory
 .RData
+.httr-oauth
diff --git a/R/vegan_wrapper.R b/R/vegan_wrapper.R
index 42457e7..adfd90b 100644
--- a/R/vegan_wrapper.R
+++ b/R/vegan_wrapper.R
@@ -1,24 +1,44 @@
 #NMDS wrapper
-vegan_wrapper <- function(RA="relative abundance",meta="metadata",category="category",labels=T){
+vegan_wrapper <- function(RA="relative abundance",meta="metadata",category_1="category",category_2=NA,labels=T){
   div=diversity(RA,index="shannon")
   if(labels){div=div[match(meta[,1],names(div))]}
   #average Shannon index by desired metadata information
-  category=which(colnames(meta)%in%category)
-  div_cat=tapply(div,meta[,category],mean)
-
+  if(!is.na(category_2)){
+    category_1=which(colnames(meta)%in%category_1)
+    category_2=which(colnames(meta)%in%category_2)
+    div_cat=tapply(div,list(meta[,category_1],meta[,category_2]),mean)
+  }
+  else{
+    category=which(colnames(meta)%in%category)
+    div_cat=tapply(div,meta[,category],mean)
+  }
   #average Pielou's evenness
   pielou=div/log(specnumber(RA))
   if(labels){pielou=pielou[match(meta[,1],names(pielou))]}
-  pie_cat=tapply(pielou,meta[,category],mean)
+  if(!is.na(category_2)){
+    pie_cat=tapply(pielou,list(meta[,category_1],meta[,category_2]),mean)
+  }
+  else{
+    pie_cat=tapply(pielou,meta[,category],mean)
+  }
+
 
   #NMDS and plot
   NMDS=metaMDS(RA,trymax=1000,autotransform=F)
   nmds_dat=as.data.frame(NMDS$points)
   if(labels){nmds_dat=nmds_dat[match(meta[,1],rownames(nmds_dat)),]}
-  NMDS_plot <- ggplot(data = nmds_dat,aes(x=MDS1,y=MDS2))+
-    geom_polygon(data=nmds_dat,aes(fill=as.factor(meta[,category])))+
-    geom_point(aes(fill=as.factor(meta[,category])))
 
+  if(!is.na(category_2)){
+    nmds_dat=cbind.data.frame(nmds_dat,"combined"=paste(meta[,category_1],meta[,category_2],sep="_"))
+    NMDS_plot <- ggplot(data = nmds_dat,aes(x=MDS1,y=MDS2))+
+    geom_polygon(data=nmds_dat,aes(fill=as.factor(combined)))+
+    geom_point(aes(fill=as.factor(combined)))
+  }
+  else{
+    NMDS_plot <- ggplot(data = nmds_dat,aes(x=MDS1,y=MDS2))+
+      geom_polygon(data=nmds_dat,aes(fill=as.factor(meta[,category_1])))+
+      geom_point(aes(fill=as.factor(meta[,category_1])))
+  }
 
     return(list("shannon"=div_cat,"pielou"=pie_cat,"NMDS"=NMDS,"NMDS_plot"=NMDS_plot))
 }