fixed some big bugs

twbattaglia · Mar 23, 2016 · 383b396 · 383b396
1 parent 98a5142
commit 383b396
Show file tree

Hide file tree

Showing 12 changed files with 164 additions and 142 deletions.
diff --git a/.gitignore b/.gitignore
@@ -13,3 +13,6 @@
 vignettes/*.html
 vignettes/*.pdf
 .Rproj.user
+
+test/
+dev/
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -3,7 +3,7 @@ Title: A suite of R function for all types of microbial diversity analyses
 Version: 0.0.1
 Authors@R: person("Thomas", "Battaglia", email = "tb1280@nyu.edu",
                   role = c("aut", "cre"))
-Description: An delicious assortment of R functions!
+Description: An assortment of R functions!
 Depends:
     R (>= 3.1.0)
 License: MIT

diff --git a/NAMESPACE b/NAMESPACE
@@ -1,3 +1,4 @@
 # Generated by roxygen2: do not edit by hand
 
-export(metagenomic_contributions)
+export(analyze_contributions)
+export(compare_alpha_diversity)
diff --git a/R/analyze_contributions.R b/R/analyze_contributions.R
@@ -7,14 +7,12 @@
 #' Any rows with duplicated row names will be dropped with the first one being
 #' kepted.
 #'
-#' @param metagenomic_contributions.py output file
-#' @param mapping file with sample metadata
-#' @param list of column variables to add to new table
-#' @return A dataframe with taxa information and sample metadata
+#' @param contributions_fp File location of the metagenomic_contributions.py output file.
+#' @param mappingfile File location of your sample metadata.
+#' @param metadata a vector of metadata identfiers which you would like to add as columns.
+#' @return A dataframe with taxa information and sample metadata for each observed OTU.
 #' @export
-analyze_contributions <- function(input_table = table_input,
-                                  mappingfile = mappingfile_input,
-                                  metadata_col = metadata_col_input){
+analyze_contributions <- function(contributions_fp, mappingfile_fp, metadata = c("Description")){
 
   # - - - - - - - - - - - - -
   # Error handling
@@ -24,22 +22,24 @@ analyze_contributions <- function(input_table = table_input,
   # mapping file is a string
   # mapping file is a dataframe
   # column id's are not strings
+  # First column is not #SampleID
 
 
   # - - - - - - - - - - - - -
   # Import files
   # - - - - - - - - - - - - -
-  # Import Neccessary Datatables from given input locations.
-  input <- read.delim(input_table, header = TRUE)
+  message("Importing files...")
+  input <- read.delim(contributions_fp, header = TRUE)
   gg_db <- otu_taxonomy_97
-  mappingfile <- read.delim(mappingfile, header = TRUE)
+  mappingfile <- read.delim(mappingfile_fp, header = TRUE)
   input.df <- dplyr::tbl_df(input)
+  input.df$OTU <- as.factor(input.df$OTU)
 
 
   # - - - - - - - - - - - - -
   # Process strings
   # - - - - - - - - - - - - -
-  message("Converting OTU-IDs's to GreenGenes Taxanomy...")
+  message("Converting identifiers to names...")
   otuid_name <- gg_db[match(input.df$OTU, table = gg_db$V1), ]$V2
   otuid_name <- gsub(pattern = ';', replacement = "", x = otuid_name, fixed = T)
   otuid_name_sep <- sapply(otuid_name, function(x) stringr::str_split(x, " "))
@@ -48,7 +48,7 @@ analyze_contributions <- function(input_table = table_input,
   # - - - - - - - - - - - - -
   # Adding Names to Table
   # - - - - - - - - - - - - -
-  message("Adding Names to Table..")
+  message("Adding names to table...")
   input.df$kingdom <- gsub(pattern = 'k__', replacement = "", x = as.character(lapply(otuid_name_sep, function(x) x[[1]])), fixed = T)
   input.df$phylum <- gsub(pattern = 'p__', replacement = "", x = as.character(lapply(otuid_name_sep, function(x) x[[2]])), fixed = T)
   input.df$class <- gsub(pattern = 'c__', replacement = "", x = as.character(lapply(otuid_name_sep, function(x) x[[3]])), fixed = T)
@@ -60,9 +60,9 @@ analyze_contributions <- function(input_table = table_input,
   # - - - - - - - - - - - - -
   # Add Metadata infomation
   # - - - - - - - - - - - - -
-  message("Collecting Metadata Information and Adding to Table..")
-  metadata_names<- unlist(stringr::str_split(metadata_col, ","))
-  metadata_tbl <- data.frame(mappingfile[match(input.df$Sample, mappingfile$X.SampleID), metadata_names])
+  message("Add metadata to table...")
+  metadata_names <- unlist(stringr::str_split(metadata, ","))
+  metadata_tbl <- mappingfile[match(input.df$Sample, mappingfile$X.SampleID), metadata]
 
 
   # - - - - - - - - - - - - -
@@ -76,5 +76,6 @@ analyze_contributions <- function(input_table = table_input,
   else{
     dplyr::bind_cols(input.df, metadata_tbl) -> input.df
   }
+  message("Completed")
   return(input.df)
 }
diff --git a/R/compare_categories.R b/R/compare_categories.R
@@ -1,79 +1,91 @@
 # @ Thomas W. Battaglia
-# @ tb1280@nyu.edu
-
-compare_beta_diversity <- function(phylo, 
-                                   x = "Day", 
-                                   group = "Treatment", 
-                                   test = "adonis", 
-                                   bdiv = "unweighted", 
-                                   write = F, 
+
+#' Calculate alpha diversity statistics
+#'
+#' Compute p-values and multiple comparisons adjusted q-values for
+#' two-group comparisons across multiple timepoints.
+#'
+#' @param phylo phyloseq object
+#' @param x
+#' @param group
+#' @param test
+#' @param bdiv
+#' @param write
+#' @param filename
+#' @return A dataframe for an PERMANOVA test over each timepoint from each two group comparison.
+compare_beta_diversity <- function(phylo,
+                                   x = "Day",
+                                   group = "Treatment",
+                                   test = "adonis",
+                                   bdiv = "unweighted",
+                                   write = F,
                                    filename = "results", ...){
-  
+
   # Load libraries
   suppressPackageStartupMessages(library(ape))
   suppressPackageStartupMessages(library(vegan))
   suppressPackageStartupMessages(library(phyloseq))
-  
-  
+
+
   # Sample metadata
   metadata <- as(sample_data(phylo), "data.frame")
-  
+
   # Get the different levels (time) for comparing
   x_levels <- levels(metadata[ ,x])
-  
+
   # Get comparing groups
   comparing_groups <- levels(metadata[ ,group])
-  
+
   # Make dataframes to store results
   if(test == "adonis"){
     results_df <- data.frame(Group1 = as.character(),
                              Group2 = as.character(),
-                             x = as.character(), 
-                             n = as.numeric(), 
+                             x = as.character(),
+                             n = as.numeric(),
                              SumsOfSqs = as.numeric(),
                              MeanSqs = as.numeric(),
                              F.Model = as.numeric(),
                              R2 = as.numeric(),
                              pvalue = as.numeric())
   }
-  
+
   if(test == "anosim"){
     results_df <- data.frame(Group1 = as.character(),
                              Group2 = as.character(),
-                             x = as.character(), 
+                             x = as.character(),
                              n = as.numeric(),
                              R_value = as.numeric(),
                              pvalue = as.numeric())
   }
 
-  
+
   # Create list of combinations
   comb_list <- combn(comparing_groups, 2, simplify = F)
-  
+
   # Loop over each comparison
   for(comp in comb_list){
-    
+
     # Subset metadata
     metadata_comparison_subset <- metadata[metadata[ ,group] %in% comp,]
     metadata_comparison_subset[ ,group] <- droplevels(metadata_comparison_subset[ ,group])
-    
+
     # Subset taxadata
     taxadata_comparison_subset <- phylo
     sample_data(taxadata_comparison_subset) <- sample_data(metadata_comparison_subset)
-    
-    
+
+
     # Loop over each level (x)
     for(i in x_levels){
-      
+
       # Subset metadata by x
       metadata_comparison_subset_time <- subset(metadata_comparison_subset, metadata_comparison_subset[ ,x] == i)
-      
+
       # Subset taxonomy by x
       taxadata_comparison_subset_time <- taxadata_comparison_subset
       sample_data(taxadata_comparison_subset_time) <- sample_data(metadata_comparison_subset_time)
-      
-      
-      
+
+
+
       ### Get Unifrac Distances ------------
       # Calculate unweighted values
       if(bdiv == "unweighted"){
@@ -84,13 +96,13 @@ compare_beta_diversity <- function(phylo,
       if(bdiv == "weighted"){
         unifrac <- UniFrac(taxadata_comparison_subset_time, weighted = TRUE, normalized = TRUE, parallel = FALSE)
       }
-      
-      
-      
+
+
+
       ### Get Stastical Values ------------
       # perform adonis test
       if(test == "adonis"){
-        
+
         # Perform Adonis test
         results <- suppressWarnings(try(adonis(formula = as.dist(unifrac) ~ metadata_comparison_subset_time[ ,group], permutations = 999)))
 
@@ -102,7 +114,7 @@ compare_beta_diversity <- function(phylo,
           F.Model <- "NA"
           R2 <- "NA"
         }
-        
+
         # If no error, assign results to variables
         if(class(results) == "adonis") {
           pval <- results$aov.tab$`Pr(>F)`[1]
@@ -111,69 +123,69 @@ compare_beta_diversity <- function(phylo,
           F.Model <- results$aov.tab$F.Model[1]
           R2 <- results$aov.tab$R2[1]
         }
-        
+
         # Find which groups are being compared.
         compared <- levels(metadata_comparison_subset_time[ ,group])
-        
+
         # Place results into dataframe
-        results_mat <- data.frame(Group1 = compared[1], 
+        results_mat <- data.frame(Group1 = compared[1],
                                   Group2 = compared[2],
-                                  x = i, 
-                                  n = nrow(metadata_comparison_subset_time), 
+                                  x = i,
+                                  n = nrow(metadata_comparison_subset_time),
                                   SumsOfSqs = SumsOfSqs,
                                   MeanSqs = MeanSqs,
                                   F.Model = F.Model,
                                   R2 = R2,
                                   pvalue = pval)
-        
+
       } # EOF Adonis
-      
+
       # perform anosim test
       if(test == "anosim"){
-        
+
         # test
         results <- suppressWarnings(try(anosim(as.dist(unifrac), metadata_comparison_subset_time[,group], permutations = 999)))
-        
+
         # If too little amount of samples are present for either group, result in None.
         if(class(results) == "try-error"){
           pval <- 'NA'
           R_value <- "NA"
         }
-        
+
         # If no error, assign results to variables
         if(class(results) == "anosim"){
           pval <- results$signif
           R_val <- results$statistic
         }
-        
+
         # Find which groups are being compared.
         compared <- levels(metadata_comparison_subset_time[ ,group])
-        
+
         # Place results into dataframe
-        results_mat <- data.frame(Group1 = compared[1], 
+        results_mat <- data.frame(Group1 = compared[1],
                                   Group2 = compared[2],
-                                  x = i, 
+                                  x = i,
                                   n = nrow(metadata_comparison_subset_time),
                                   R_value = R_val,
                                   pvalue = pval)
       } # EOF Anosim
 
       # Merge results
       results_df <- rbind(results_mat, results_df)
-    
+
     } # EOF levels
-    
+
   } # EOF comparisons
-  
+
   # Write table if chosen
   if (write == T){
     write.csv(x = results_df, file = paste(filename, ".csv", sep = ""))
   }
-  
+
   #if (write_plots == T){
     #return(list_of_plots)
   #}
-  
+
   return(results_df)
 }
 
diff --git a/R/plot_contributions.R b/R/plot_contributions.R