updates

DESCRIPTION

+Package: metastatistics
+Type: Package
+Title: What the Package Does (Title Case)
+Version: 0.1.0
+Author: Who wrote it
+Maintainer: The package maintainer <yourself@somewhere.net>
+Description: More about what it does (maybe more than one line)
+    Use four spaces when indenting paragraphs within the Description.
+License: What license is it under?
+Encoding: UTF-8
+LazyData: true
+RoxygenNote: 6.1.1
+Imports: vegan, ggplot2, magrittr, ggalt

NAMESPACE

+# Generated by roxygen2: do not edit by hand
+
+export(Filter_abundance)
+export(Filter_goods)
+export(Filter_kingdom)
+export(Generate_tree)
+export(Goods_cov)
+export(PCOA)
+export(Remove_kingdom)

R/Filter.R

+#' @export
+Filter_goods <- function(data_raw, cutoff = 0.9){
+  # Goods
+  cov_matrix <- Goods_cov(data_raw)
+  good_samples <- cov_matrix[cov_matrix[, "Good's coverage"] %>% as.numeric %>% `>=` (cutoff) %>% which, 1]
+  bad_samples <- colnames(data_raw@otu_table@.Data)[-which(colnames(data_raw@otu_table@.Data) %in% good_samples)]
+  data_raw@otu_table@.Data <- data_raw@otu_table@.Data[,which(colnames(data_raw@otu_table@.Data) %in% good_samples)]
+  return(data_raw)
+}
+
+#' @export
+Filter_abundance <- function(data_raw, cutoff = 0.005){
+  #Abundance
+  total_abundance <- data_raw@otu_table@.Data %>% sum
+
+  filter_low_abundance <- function(otu, threshold = 0.005){
+    rsum <- sum(otu)
+    perc <- rsum/total_abundance*100
+    if (perc < threshold){
+      return(F)
+    }else{
+      return(T)
+    }
+  }
+
+  abundances <- apply(data_raw@otu_table@.Data, 1, filter_low_abundance)
+  old_table <- data_raw@otu_table@.Data
+  keep <- rownames(old_table[abundances,])
+  data_new <- prune_taxa(keep, data_raw)
+
+  return(data_new)
+}
+
+#' @export
+Filter_kingdom <- function(data_raw, col_name, desired_kingdom){
+
+  good_kingdoms <- apply(data_raw@tax_table@.Data, 1, function(otu){
+    if(is.na(otu[col_name])){
+      return(F)
+    }
+    if (otu[col_name] == desired_kingdom){
+      return(T)
+    } else{
+      return(F)
+    }
+  })
+
+  only_good <- good_kingdoms[which(good_kingdoms == T)] %>% names
+  available_good <- only_good[which(only_good %in% rownames(data_raw@otu_table@.Data))]
+  data_raw@otu_table@.Data <- data_raw@otu_table@.Data[available_good,]
+
+  return(data_raw)
+}
+
+#' @export
+Remove_kingdom <- function(data_raw, col_name, undesired_kingdom){
+
+  keep <- which(data_raw@tax_table@.Data[,col_name] != undesired_kingdom) %>% names
+
+
+  return(prune_taxa(keep, data_raw))
+}
+
+#' @export
+Generate_tree <- function(phy_obj, otu_fasta, n_proc = 4){
+  library(seqinr)
+  library(DECIPHER)
+  library(phangorn)
+
+  if (missing(otu_fasta)){
+    print("True")
+    if (phy_obj@otu_table@taxa_are_rows){
+      seqs_vec <-  rownames(phy_obj@otu_table@.Data)
+    } else{
+      seqs_vec <- colnames(phy_obj@otu_table@.Data)
+    }
+
+    names(seqs_vec) <- seqs_vec
+  } else{
+    seqs <- read.fasta(otu_fasta)
+    seqs_vec <- sapply(names(seqs), function(otu){
+      s <- seqs[[otu]]
+      paste0(s, collapse = '')
+    })
+  }
+
+  alignment <- DECIPHER::AlignSeqs(DNAStringSet(seqs_vec), anchor=NA, processors = n_proc)
+
+  phang.align <- phyDat(as(alignment,"matrix"), type="DNA")
+  dm <- dist.ml(phang.align)
+  treeNJ <- NJ(dm)# Note, tip order != sequence order
+  fit= pml(treeNJ,data=phang.align)
+  itGTR<- update(fit,k=4,inv=0.2)
+  fitGTR<- optim.pml(fitGTR, model="GTR",optInv=TRUE,optGamma=TRUE,rearrangement="stochastic",control=pml.control(trace=0))
+
+  phy_obj@phy_tree <- phy_tree(fitGTR$tree)
+
+  return(phy_obj)
+}
+

R/Satistics.R

+#' @export
+Goods_cov <- function(data_raw){
+  cov_matrix <- matrix(ncol=2,nrow=0)
+  for (sample in colnames(data_raw@otu_table@.Data)){
+    sample_col <- data_raw@otu_table@.Data[, sample]
+    total_individuals <- sample_col %>% sum
+    total_singletons <- sample_col[which(sample_col == 1)] %>% sum
+    goods_cov <- 1 - (total_singletons / total_individuals)
+    cov_matrix <- rbind(cov_matrix, c(sample, goods_cov))
+  }
+
+  colnames(cov_matrix) <- c("Sample", "Good's coverage")
+  return(cov_matrix)
+}
+
+#' @export
+PCOA <- function(data_phylo, title, color, shape){
+  ordination <- ordinate(data_phylo, method="PCoA", distance="bray")
+
+  if(missing(shape)){
+    return(plot_ordination(data_phylo, ordination, color = color, title = title) + stat_ellipse())
+  } else{
+
+    return(plot_ordination(data_phylo, ordination, color = color, shape = shape, title = title) + stat_ellipse())
+  }
+
+}