Grades of membership model finalized

Summary

Load data post filtering, normalization and batch correction (nasal_filtered_normed_batchcorrected.rds).

Apply grades of membership model using CountClust. We chose a random seed based on a series of analysis comparing goodness-of-fit of the model under 100 random seeds (https://jhsiao999.github.io/nasalmicrobiome/gom-seeds.html). The best random seed produces the largest Bayes Factor and is used for running the Grades of Membership (GoM) model.

\(~\)

Loading data and packages

library(knitr)
library(kableExtra)
library(dplyr)
library(metagenomeSeq)
library(CountClust)

Read in filtered data.

MRobj = readRDS("../data/nasal_filtered_normed_batchcorrected.rds")
MRobj

MRexperiment (storageMode: environment)
assayData: 553 features, 197 samples 
  element names: counts 
protocolData: none
phenoData
  sampleNames: EM0042 EM0047 ... E0168 (197 total)
  varLabels: StudyID age ... infnone (94 total)
  varMetadata: labelDescription
featureData
  featureNames:
    k__Bacteria;p__[Thermi];c__Deinococci;o__Deinococcales;f__Deinococcaceae;g__CM44;s__
    k__Bacteria;p__[Thermi];c__Deinococci;o__Deinococcales;f__Deinococcaceae;g__Deinococcus;s__
    ... s__zeae (553 total)
  fvarLabels: Kingdom Phylum ... Species (7 total)
  fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:  

\(~\)

Clustering

Apply CountClust to fit grades of membership model.

set.seed(67)

counts <- MRcounts(MRobj,norm=FALSE,log=FALSE)
counts_clust <- FitGoM(t(counts), K=c(2), tol=0.00001)
saveRDS(counts_clust, file = "../data/count-clust.rds")

Load previous computed count_clust results.

counts_clust <- readRDS(file = "../data/count-clust.rds")

cct <- counts_clust[[1]]$omega
rowMaxs = sapply(seq(nrow(cct)),function(i){
    seq(along=cct[i,])[cct[i,]==max(cct[i,])]
  })
rowMaxs = unlist(rowMaxs)
uniqueGroups = unique(unlist(rowMaxs))
groupIds = split(seq(rowMaxs),factor(rowMaxs))

cl = rep("1",sum(sapply(groupIds,length))); cl[groupIds[[2]]] = 2

table(cl)

cl
  1   2 
167  30 

cbind(groupIDs = cl,round(cct,3)) %>% kable %>%
  kable_styling()

\(~\)

Preserve data with GOM

Next we save the dataset with the sample clustering labels

# MRobj2 = MRobj
# assayData(MRobj2)[["counts"]] = counts
# featureData(MRobj2) = AnnotatedDataFrame(data.frame(Species=rownames(counts)))
pData(MRobj)$GOM = cl
#pData(MRobj2)$GOM = cl

saveRDS(MRobj, file="../data/nasal_GOM.rds")

\(~\)

Session Info

R version 3.4.1 (2017-06-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS High Sierra 10.13

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] CountClust_1.4.1     ggplot2_2.2.1        metagenomeSeq_1.20.1
 [4] RColorBrewer_1.1-2   glmnet_2.0-13        foreach_1.4.4       
 [7] Matrix_1.2-12        limma_3.34.9         Biobase_2.38.0      
[10] BiocGenerics_0.24.0  dplyr_0.7.4          kableExtra_0.7.0    
[13] knitr_1.20          

loaded via a namespace (and not attached):
 [1] maptpx_1.9-3       Rcpp_0.12.15       ape_5.0           
 [4] lattice_0.20-35    gtools_3.5.0       assertthat_0.2.0  
 [7] rprojroot_1.3-2    digest_0.6.15      slam_0.1-42       
[10] R6_2.2.2           plyr_1.8.4         backports_1.1.2   
[13] stats4_3.4.1       evaluate_0.10.1    highr_0.6         
[16] httr_1.3.1         pillar_1.2.1       gplots_3.0.1      
[19] rlang_0.2.0        lazyeval_0.2.1     gdata_2.18.0      
[22] vegan_2.4-6        rmarkdown_1.9      readr_1.1.1       
[25] stringr_1.3.0      munsell_0.4.3      compiler_3.4.1    
[28] pkgconfig_2.0.1    SQUAREM_2017.10-1  mgcv_1.8-23       
[31] htmltools_0.3.6    nnet_7.3-12        tibble_1.4.2      
[34] codetools_0.2-15   matrixStats_0.53.1 permute_0.9-4     
[37] viridisLite_0.3.0  MASS_7.3-49        bitops_1.0-6      
[40] grid_3.4.1         nlme_3.1-131.1     gtable_0.2.0      
[43] git2r_0.21.0       magrittr_1.5       scales_0.5.0      
[46] KernSmooth_2.23-15 stringi_1.1.6      reshape2_1.4.3    
[49] flexmix_2.3-14     bindrcpp_0.2       xml2_1.2.0        
[52] boot_1.3-20        cowplot_0.9.2      iterators_1.0.9   
[55] tools_3.4.1        picante_1.6-2      glue_1.2.0        
[58] hms_0.4.1          yaml_2.1.18        colorspace_1.3-2  
[61] cluster_2.0.6      caTools_1.17.1     rvest_0.3.2       
[64] bindr_0.1          modeltools_0.2-21