Community State Types (CST)

Acknowledgement The main idea of the microbial analysis and the phyloseq dataset are coming from this paper Temporal and spatial variation of the human microbiota during pregnancy The idea which is underlying this analysis is to use clustering to discover the dynamic of the gut community via Community State Types (CST). CSTs are discrete set of clusters.

Initialize:

## [1] '1.34.0'

## [1] '3.3.5'

## [1] '2.1.2'

## [1] '1.2.6'

## [1] '0.8.6'

• Phyloseq dataset From the same paper Temporal and spatial variation of the human microbiota during pregnancy

Dataset from stool

da<- load('RDA/PregnancyClosed15.Rdata')
site <- "Stool"
phy <- PSPreg[[site]]
phy

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1271 taxa and 580 samples ]
## sample_data() Sample Data:       [ 580 samples by 64 sample variables ]
## tax_table()   Taxonomy Table:    [ 1271 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 1271 tips and 1270 internal nodes ]

Remove rare taxa Transform the data (proportions):

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 894 taxa and 580 samples ]
## sample_data() Sample Data:       [ 580 samples by 64 sample variables ]
## tax_table()   Taxonomy Table:    [ 894 taxa by 7 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 894 tips and 893 internal nodes ]

We are not doing differential abundance analysis here, so the proportion transformation is used for exploratory analyses only.

Cluster into CSTs

The vaginal community is dominated by closely related, but functionally distinct, Lactobacillus species. Therefore it is better to use a non-phylogenetically aware distance measure so as to be able to separate these species. Start with an MDS (or PCoA) ordination:

##  [1] 42.210864 22.612218 12.083443  9.496256  8.919936  6.947028  6.888458
##  [8]  5.240856  4.622571  4.539161  4.246411  3.870323  3.374326  3.117964
## [15]  2.861045  2.852497  2.468666  2.355515  2.090535  1.991523

## [1] -0.3676292 -0.3957017 -0.4337434 -0.4744555 -0.5301270 -0.8299508

Denoise distance matrix

We would like to clean some of the noise from the data by restricting this to the truly significant dimensions. The top 5 eigenvalues are clearly very significant, but let’s keep all the positive eigenvalues that clearly exceed the magnitude of the smallest negative eigenvalues:

h_sub5 <- hist(evs[6:length(evs)], 100)

plot(h_sub5$mids, h_sub5$count, log="y", type='h', lwd=10, lend=2)

Looks like eigenvalues 6 and 7 still stand out, so we’ll go with 7 MDS dimensions.

Determine number of clusters

We will use the gap statistic to indicate the number of clusters in this data:

The gap statistic strongly suggests at least three clusters, but makes another big jump at K=5 before the slope gets a lot smaller. So, K=5 it is.

Cluster into CSTs

Perform PAM 5-fold clusters:

Evaluate clustering

Inspect the results in MDS and NMDS ordinations:

## Run 0 stress 0.18648 
## Run 1 stress 0.190787 
## Run 2 stress 0.1862701 
## ... New best solution
## ... Procrustes: rmse 0.01168006  max resid 0.09976727 
## Run 3 stress 0.1933471 
## Run 4 stress 0.201474 
## Run 5 stress 0.1954427 
## Run 6 stress 0.1922218 
## Run 7 stress 0.1879763 
## Run 8 stress 0.1921214 
## Run 9 stress 0.1884075 
## Run 10 stress 0.2235917 
## Run 11 stress 0.1919761 
## Run 12 stress 0.198366 
## Run 13 stress 0.1877754 
## Run 14 stress 0.190696 
## Run 15 stress 0.187102 
## Run 16 stress 0.1891335 
## Run 17 stress 0.1956618 
## Run 18 stress 0.1917857 
## Run 19 stress 0.1842801 
## ... New best solution
## ... Procrustes: rmse 0.01003169  max resid 0.09694101 
## Run 20 stress 0.200065 
## *** No convergence -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     11: stress ratio > sratmax
##      8: scale factor of the gradient < sfgrmin

Heatmaps for taxa per clusters

The ordinations offer support for these being legitimate clusters, even if they are not perfect and some samples look like they might be mixtures of two clusters. Let’s take a look at the heatmaps of each cluster for additional perspective:

Boxplots

Boxplots for each Species in the five clusters. We will investigate shannon index with non parametric Kruskal test for each taxa among the five clusters

Post hoc

Post hoc by implementing Conover’s test for the significant p_values using shannon index with adjusted p_value is (Benjamini-Hochberg).