16S downstream analysis in colorectal carcinoma

1. Install and/or load the required packages

# install.packages("BiocManager")
# BiocManager::install(c("phyloseq", "vegan", "microbiome", "Maaslin2"))
# BiocManager::install(c("microbiomeMarker"))
# BiocManager::install(c("ggplot2", "tidyverse", "dplyr", "data.table"))

library("phyloseq")
library("vegan")
library("ggplot2")
library("microbiomeMarker")
library("tidyverse")
library("dplyr")
library("data.table")
library("microbiome")
library("Maaslin2")
library("igraph")
library("RColorBrewer")
library("nlme")

2. Set the working directory

setwd("C://Users//Malina//Desktop//Microbial compositionality")

3. Load the example dataset

data(kostic_crc)
kostic_crc

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 2505 taxa and 177 samples ]
sample_data() Sample Data:       [ 177 samples by 71 sample variables ]
tax_table()   Taxonomy Table:    [ 2505 taxa by 7 taxonomic ranks ]

## Explore the dataset
dim(sample_data(kostic_crc))

[1] 177  71

head(sample_data(kostic_crc))

Sample Data:        [6 samples by 71 sample variables]:

table(sample_data(kostic_crc)$DIAGNOSIS)

Healthy   Tumor 
     91      86 

head(tax_table(kostic_crc))

Taxonomy Table:     [6 taxa by 7 taxonomic ranks]:
       Kingdom    Phylum           Class                 Order           Family               
304309 "Bacteria" "Firmicutes"     "Clostridia"          "Clostridiales" "Lachnospiraceae"    
469478 "Bacteria" "Firmicutes"     "Clostridia"          "Clostridiales" "Lachnospiraceae"    
208196 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"   "Methylobacteriaceae"
358030 "Bacteria" "Firmicutes"     "Clostridia"          "Clostridiales" NA                   
16076  "Bacteria" "Firmicutes"     "Clostridia"          "Clostridiales" "Ruminococcaceae"    
35786  "Bacteria" "Firmicutes"     "Clostridia"          "Clostridiales" "Ruminococcaceae"    
       Genus              Species                        
304309 "Blautia"          "Blautia producta"             
469478 "Catonella"        NA                             
208196 "Methylobacterium" "Methylobacterium organophilum"
358030 NA                 NA                             
16076  "Ruminococcus"     "Ruminococcus bromii"          
35786  NA                 NA                             

head(otu_table(kostic_crc)[,1:20])

OTU Table:          [6 taxa and 20 samples]
                     taxa are rows
       C0333.N.518126 C0333.T.518046 38U4VAHB.518100 XZ33PN7O.518030 GQ6LSNI6.518106
304309             40              4               1               2              30
469478              0              0               0               0               0
208196              0              0               0               0               0
358030              0              0               0               0               0
16076             271             28             110               7              34
35786               0              0               0               0               0
       C0270.N.518041 HZIMMAM6.518095 LRQ9WN6C.518048 A8A34NCW.518023 C0332.N.518027
304309              1               4               8               7              2
469478              0               0               0               0              0
208196              0               0               0               0              0
358030              0               0               0               0              0
16076               0               0               0               0              8
35786               0               0               0               0              0
       C10.S.517995 MQE9ONGV.518038 C0230.N.517992 C0256.N.518170 HZIMMNL5.518062
304309            2               2              1              1               4
469478            0               0              0              0               0
208196            0               0              0              0               0
358030            0               0              0              0               0
16076            21               0              3             11               0
35786             0               0              0              0               0
       GQ6LSABJ.518010 C0282.N.518138 TV28INUO.518004 C0315.N.518021 C0235.N.517993
304309               7              1               1              1              1
469478               0              0               1              0              0
208196               0              0               0              0              0
358030               0              0               0              0              0
16076                9             13               2              7              0
35786                0              0               0              0              0

4. Function to create a phyloseq object from data files

create_phyloseq <- function(otu_file, sample_file, tax_file) {
  # Load OTU table
  otu_data <- read.csv(otu_file, row.names = 1, stringsAsFactors = FALSE)
  otu_matrix <- as.matrix(otu_data)
  
  # Load sample data
  sample_data <- read.csv(sample_file, row.names = 1)
  
  # Load taxonomic table
  tax_data <- read.csv(tax_file, row.names = 1, stringsAsFactors = FALSE)
  tax_table <- tax_data[, 1:7]  # Assuming taxonomic ranks are in columns 1 to 7
  
  # Create phyloseq object
  physeq_obj <- phyloseq(otu_table(otu_matrix, taxa_are_rows = TRUE),
                         sample_data,
                         tax_table)
  
  return(physeq_obj)
}

# this is how the function is run

# # Provide file paths
# otu_file <- "path/to/otu_table.csv"
# sample_file <- "path/to/sample_data.csv"
# tax_file <- "path/to/taxonomic_table.csv"
# 
# # Create phyloseq object
# physeq_obj <- create_phyloseq(otu_file, sample_file, tax_file)

5. Normalization of the sample counts

## check how large the sample sizes (number of counts) are 

max_difference = max(sample_sums(kostic_crc))/min(sample_sums(kostic_crc))
max_difference

[1] 21.34924

max(sample_sums(kostic_crc))

[1] 11187

min(sample_sums(kostic_crc))

[1] 524

## norm is only needed if sample sizes are more than 10x different (and they are)

## distribution of the sample sizes
hist(sample_sums(kostic_crc), breaks = 177)

## get relative abundance with the microbiome package
kostic_crc.compositional <- microbiome::transform(kostic_crc, "compositional")

head(otu_table(kostic_crc.compositional)[,1:20])

OTU Table:          [6 taxa and 20 samples]
                     taxa are rows
       C0333.N.518126 C0333.T.518046 38U4VAHB.518100 XZ33PN7O.518030 GQ6LSNI6.518106
304309    0.007078393     0.00311042     0.000152765    0.0004374453     0.003442736
469478    0.000000000     0.00000000     0.000000000    0.0000000000     0.000000000
208196    0.000000000     0.00000000     0.000000000    0.0000000000     0.000000000
358030    0.000000000     0.00000000     0.000000000    0.0000000000     0.000000000
16076     0.047956114     0.02177294     0.016804155    0.0015310586     0.003901767
35786     0.000000000     0.00000000     0.000000000    0.0000000000     0.000000000
       C0270.N.518041 HZIMMAM6.518095 LRQ9WN6C.518048 A8A34NCW.518023 C0332.N.518027
304309   0.0001247194    0.0004288624     0.001371507      0.00100416    0.001091107
469478   0.0000000000    0.0000000000     0.000000000      0.00000000    0.000000000
208196   0.0000000000    0.0000000000     0.000000000      0.00000000    0.000000000
358030   0.0000000000    0.0000000000     0.000000000      0.00000000    0.000000000
16076    0.0000000000    0.0000000000     0.000000000      0.00000000    0.004364430
35786    0.0000000000    0.0000000000     0.000000000      0.00000000    0.000000000
       C10.S.517995 MQE9ONGV.518038 C0230.N.517992 C0256.N.518170 HZIMMNL5.518062
304309 0.0004409171    0.0007501875   0.0009823183   0.0001157273    0.0006255865
469478 0.0000000000    0.0000000000   0.0000000000   0.0000000000    0.0000000000
208196 0.0000000000    0.0000000000   0.0000000000   0.0000000000    0.0000000000
358030 0.0000000000    0.0000000000   0.0000000000   0.0000000000    0.0000000000
16076  0.0046296296    0.0000000000   0.0029469548   0.0012730008    0.0000000000
35786  0.0000000000    0.0000000000   0.0000000000   0.0000000000    0.0000000000
       GQ6LSABJ.518010 C0282.N.518138 TV28INUO.518004 C0315.N.518021 C0235.N.517993
304309         0.00112   0.0005980861    0.0002352941   0.0001344267   0.0006609385
469478         0.00000   0.0000000000    0.0002352941   0.0000000000   0.0000000000
208196         0.00000   0.0000000000    0.0000000000   0.0000000000   0.0000000000
358030         0.00000   0.0000000000    0.0000000000   0.0000000000   0.0000000000
16076          0.00144   0.0077751196    0.0004705882   0.0009409867   0.0000000000
35786          0.00000   0.0000000000    0.0000000000   0.0000000000   0.0000000000

## get relative abundance with the microbiomeMarker package
kostic_crc_TSS = microbiomeMarker::normalize(kostic_crc, method="TSS")
head(otu_table(kostic_crc_TSS)[,1:20])

OTU Table:          [6 taxa and 20 samples]
                     taxa are rows
       C0333.N.518126 C0333.T.518046 38U4VAHB.518100 XZ33PN7O.518030 GQ6LSNI6.518106
304309    0.007078393     0.00311042     0.000152765    0.0004374453     0.003442736
469478    0.000000000     0.00000000     0.000000000    0.0000000000     0.000000000
208196    0.000000000     0.00000000     0.000000000    0.0000000000     0.000000000
358030    0.000000000     0.00000000     0.000000000    0.0000000000     0.000000000
16076     0.047956114     0.02177294     0.016804155    0.0015310586     0.003901767
35786     0.000000000     0.00000000     0.000000000    0.0000000000     0.000000000
       C0270.N.518041 HZIMMAM6.518095 LRQ9WN6C.518048 A8A34NCW.518023 C0332.N.518027
304309   0.0001247194    0.0004288624     0.001371507      0.00100416    0.001091107
469478   0.0000000000    0.0000000000     0.000000000      0.00000000    0.000000000
208196   0.0000000000    0.0000000000     0.000000000      0.00000000    0.000000000
358030   0.0000000000    0.0000000000     0.000000000      0.00000000    0.000000000
16076    0.0000000000    0.0000000000     0.000000000      0.00000000    0.004364430
35786    0.0000000000    0.0000000000     0.000000000      0.00000000    0.000000000
       C10.S.517995 MQE9ONGV.518038 C0230.N.517992 C0256.N.518170 HZIMMNL5.518062
304309 0.0004409171    0.0007501875   0.0009823183   0.0001157273    0.0006255865
469478 0.0000000000    0.0000000000   0.0000000000   0.0000000000    0.0000000000
208196 0.0000000000    0.0000000000   0.0000000000   0.0000000000    0.0000000000
358030 0.0000000000    0.0000000000   0.0000000000   0.0000000000    0.0000000000
16076  0.0046296296    0.0000000000   0.0029469548   0.0012730008    0.0000000000
35786  0.0000000000    0.0000000000   0.0000000000   0.0000000000    0.0000000000
       GQ6LSABJ.518010 C0282.N.518138 TV28INUO.518004 C0315.N.518021 C0235.N.517993
304309         0.00112   0.0005980861    0.0002352941   0.0001344267   0.0006609385
469478         0.00000   0.0000000000    0.0002352941   0.0000000000   0.0000000000
208196         0.00000   0.0000000000    0.0000000000   0.0000000000   0.0000000000
358030         0.00000   0.0000000000    0.0000000000   0.0000000000   0.0000000000
16076          0.00144   0.0077751196    0.0004705882   0.0009409867   0.0000000000
35786          0.00000   0.0000000000    0.0000000000   0.0000000000   0.0000000000

## output of both from above are relative abundances

## css norm (originally from metagenomeSeq packege) ## 
# Cumulative Sum Scaling (CSS) is a median-like quantile normalization which corrects differences in sampling depth (library size).
# While standard relative abundance (fraction/percentage) normalization re-scales all samples to the same total sum (100%), CSS keeps a variation in total counts between samples. 
# CSS re-scales the samples based on a subset (quartile) of lower abundant taxa 
#(relatively constant and independent), 
# thereby excluding the impact of (study dominating) high abundant taxa.

kostic_crc_CSS = microbiomeMarker::normalize(kostic_crc, method="CSS")

Default value being used.

head(otu_table(kostic_crc_CSS)[,1:20])

OTU Table:          [6 taxa and 20 samples]
                     taxa are rows
       C0333.N.518126 C0333.T.518046 38U4VAHB.518100 XZ33PN7O.518030 GQ6LSNI6.518106
304309             40              4               1               2              30
469478              0              0               0               0               0
208196              0              0               0               0               0
358030              0              0               0               0               0
16076             271             28             110               7              34
35786               0              0               0               0               0
       C0270.N.518041 HZIMMAM6.518095 LRQ9WN6C.518048 A8A34NCW.518023 C0332.N.518027
304309              1               4               8               7              2
469478              0               0               0               0              0
208196              0               0               0               0              0
358030              0               0               0               0              0
16076               0               0               0               0              8
35786               0               0               0               0              0
       C10.S.517995 MQE9ONGV.518038 C0230.N.517992 C0256.N.518170 HZIMMNL5.518062
304309            2               2              1              1               4
469478            0               0              0              0               0
208196            0               0              0              0               0
358030            0               0              0              0               0
16076            21               0              3             11               0
35786             0               0              0              0               0
       GQ6LSABJ.518010 C0282.N.518138 TV28INUO.518004 C0315.N.518021 C0235.N.517993
304309               7              1               1              1              1
469478               0              0               1              0              0
208196               0              0               0              0              0
358030               0              0               0              0              0
16076                9             13               2              7              0
35786                0              0               0              0              0

## cpm norm originally lefse package ###
kostic_crc_CPM = microbiomeMarker::normalize(kostic_crc, method="CPM")

head(otu_table(kostic_crc_CPM)[,1:20])

OTU Table:          [6 taxa and 20 samples]
                     taxa are rows
       C0333.N.518126 C0333.T.518046 38U4VAHB.518100 XZ33PN7O.518030 GQ6LSNI6.518106
304309       7078.393        3110.42         152.765        437.4453        3442.736
469478          0.000           0.00           0.000          0.0000           0.000
208196          0.000           0.00           0.000          0.0000           0.000
358030          0.000           0.00           0.000          0.0000           0.000
16076       47956.114       21772.94       16804.155       1531.0586        3901.767
35786           0.000           0.00           0.000          0.0000           0.000
       C0270.N.518041 HZIMMAM6.518095 LRQ9WN6C.518048 A8A34NCW.518023 C0332.N.518027
304309       124.7194        428.8624        1371.507         1004.16       1091.107
469478         0.0000          0.0000           0.000            0.00          0.000
208196         0.0000          0.0000           0.000            0.00          0.000
358030         0.0000          0.0000           0.000            0.00          0.000
16076          0.0000          0.0000           0.000            0.00       4364.430
35786          0.0000          0.0000           0.000            0.00          0.000
       C10.S.517995 MQE9ONGV.518038 C0230.N.517992 C0256.N.518170 HZIMMNL5.518062
304309     440.9171        750.1875       982.3183       115.7273        625.5865
469478       0.0000          0.0000         0.0000         0.0000          0.0000
208196       0.0000          0.0000         0.0000         0.0000          0.0000
358030       0.0000          0.0000         0.0000         0.0000          0.0000
16076     4629.6296          0.0000      2946.9548      1273.0008          0.0000
35786        0.0000          0.0000         0.0000         0.0000          0.0000
       GQ6LSABJ.518010 C0282.N.518138 TV28INUO.518004 C0315.N.518021 C0235.N.517993
304309            1120       598.0861        235.2941       134.4267       660.9385
469478               0         0.0000        235.2941         0.0000         0.0000
208196               0         0.0000          0.0000         0.0000         0.0000
358030               0         0.0000          0.0000         0.0000         0.0000
16076             1440      7775.1196        470.5882       940.9867         0.0000
35786                0         0.0000          0.0000         0.0000         0.0000

#tail(otu_table(kostic_crc_CSS))


## rarefy norm ###
## from microbiomeMarker package ##
# 
# kostic_crc_rare = microbiomeMarker::normalize(kostic_crc, method="rarefy")
# 
# head(otu_table(kostic_crc_rare)[,1:20])

## from phyloseq package 

ps.rarefied = phyloseq::rarefy_even_depth(kostic_crc, rngseed=123, sample.size=524, replace=F)

`set.seed(123)` was used to initialize repeatable random subsampling.
Please record this for your records so others can reproduce.
Try `set.seed(123); .Random.seed` for the full vector
...
1110OTUs were removed because they are no longer 
present in any sample after random subsampling

...

6. Plot stack bar plots of the top 15 taxa based on abund

##  extract the top 15 taxa on this taxa level for CSS 
top_phy = tax_glom(kostic_crc_CSS, "Genus")
top15 = names(sort(taxa_sums(top_phy), decreasing=TRUE)[1:15])
top15_css_ps = prune_taxa(top15, top_phy)


phyloseq::plot_bar(top15_css_ps, "DIAGNOSIS", fill="Genus")+
  theme( axis.ticks.x=element_blank(), panel.background=element_rect(fill=NA), panel.grid.major=element_line(colour="#ebebeb")) + 
  labs(x=NULL)  

## get rid of the black squares and get better colors
plot_bar(top15_css_ps, "DIAGNOSIS", fill="Genus")+
  theme( axis.ticks.x=element_blank(), panel.background=element_rect(fill=NA), panel.grid.major=element_line(colour="#ebebeb")) + 
  labs(x=NULL) + geom_bar(stat="identity") 

## Generate 15 distinct colors from the palette

color = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)]
sampled_colors <- sample(color, 15)

plot_bar(top15_css_ps, "DIAGNOSIS", fill = "Genus") +
  theme(axis.ticks.x = element_blank(), panel.background = element_rect(fill = NA), panel.grid.major = element_line(colour = "#ebebeb")) +
  labs(x = NULL) +
  geom_bar(stat = "identity", aes(fill = Genus)) +  # Specify fill aesthetic here
  scale_fill_manual(values = sampled_colors )  # Use scale_fill_manual for fill colors

## plot each genus in a separate facet

# plot_bar(top15_css_ps, "DIAGNOSIS", fill="DIAGNOSIS", facet_grid=~Genus)  + # # #geom_bar(stat="identity")

## plot top 5 ##

top_phy5 = tax_glom(kostic_crc_CSS, "Genus")
top5 = names(sort(taxa_sums(top_phy5), decreasing=TRUE)[1:5])
top5_css_ps = prune_taxa(top5, top_phy5)

tax_table(top5_css_ps)

Taxonomy Table:     [5 taxa by 7 taxonomic ranks]:
       Kingdom    Phylum          Class                  Order             Family            
180285 "Bacteria" "Firmicutes"    "Clostridia"           "Clostridiales"   "Ruminococcaceae" 
470429 "Bacteria" "Fusobacteria"  "Fusobacteria (class)" "Fusobacteriales" "Fusobacteriaceae"
49274  "Bacteria" "Firmicutes"    "Clostridia"           "Clostridiales"   "Lachnospiraceae" 
332077 "Bacteria" "Bacteroidetes" "Bacteroidia"          "Bacteroidales"   "Prevotellaceae"  
469709 "Bacteria" "Bacteroidetes" "Bacteroidia"          "Bacteroidales"   "Bacteroidaceae"  
       Genus              Species
180285 "Faecalibacterium" NA     
470429 "Fusobacterium"    NA     
49274  "Clostridium"      NA     
332077 "Prevotella"       NA     
469709 "Bacteroides"      NA     

## plot each genus in a separate facet
plot_bar(top5_css_ps, "DIAGNOSIS", fill="DIAGNOSIS", facet_grid=~Genus) + geom_bar(stat="identity")

## changing the appearance 
## https://stackoverflow.com/questions/63133884/removing-the-original-color-outline-in-r-when-using-a-new-pallette-in-a-barplot
## remove black boxes 
## https://github.com/joey711/phyloseq/issues/721

##  extract the top 15 taxa on this taxa level for TSS
top_phy = tax_glom(kostic_crc_TSS, "Genus")
top15 = names(sort(taxa_sums(top_phy), decreasing=TRUE)[1:15])
top15_tss_ps = prune_taxa(top15, top_phy)


# prepare to create a nice plot for saving
top15_tss_G  =plot_bar(top15_tss_ps, "DIAGNOSIS", fill="Genus")+
  theme( axis.ticks.x=element_blank(), panel.background=element_rect(fill=NA), panel.grid.major=element_line(colour="#ebebeb")) + 
  labs(x=NULL)+
  geom_bar(stat = "identity", aes(fill = Genus)) +  
  scale_fill_manual(values = sampled_colors ) +
  theme_bw() +
  theme(
    text = element_text(size = 12, face = "bold"),  # Make all text bold
    axis.title = element_text(size = 14, face = "bold"),  # Make axis titles bold
    axis.text = element_text(size = 12),  # Adjust axis text size
    legend.title = element_text(size = 14, face = "bold"),  # Make legend title bold
    legend.text = element_text(size = 12),  # Adjust legend text size
    plot.title = element_text(size = 16, face = "bold")  # Make plot title bold
  )

  
top15_tss_G

#save the image
#ggsave("top15_Genera_plot_TSS.png", top15_tss_G, width = 6, height = 8, dpi = 300)

# plot rarefaction curve 
# Extract OTU abundance data and transpose
# otu_table_transposed<- as.data.frame(t(otu_table(kostic_crc)))
# head(otu_table_transposed)
# rarecurve(as.matrix(otu_table_transposed), step=500, cex=0.1)

7. Plot alpha diversity between groups

## estimation of alpha div only works with integers and not rel abund values (float numbers)
## what each of the alpha div means 
## Observed - counts the number of unique species (OTUs or ASVs) present in a sample
## Chao1 -  It takes into account the number of rare species (singletons and doubletons) observed in the sample and extrapolates the total richness based on the frequency of these rare species
## ACE -  Takes into account both the number of rare species observed and their relative abundances. inflates the number of rare taxa and inflates again the number of taxa with abundance 1.
## Shannon's diversity index (also known as Shannon-Wiener index) considers both species richness and evenness in a community. It takes into account the number of different species present as well as the relative abundance of each species.
## Simpson's diversity index measures the probability that two randomly selected sequences are of the same species


## CSS ##
alpha_div = plot_richness(kostic_crc_CSS, x="DIAGNOSIS", color="DIAGNOSIS", title="CSS", measures=c("Observed", "Chao1", "ACE", "Simpson", "Shannon")) + 
  geom_boxplot() +
  theme_bw() +
  theme(
    text = element_text(size = 12, face = "bold"),  # Make all text bold
    axis.title = element_text(size = 14, face = "bold"),  # Make axis titles bold
    axis.text = element_text(size = 12),  # Adjust axis text size
    legend.title = element_text(size = 14, face = "bold"),  # Make legend title bold
    legend.text = element_text(size = 12),  # Adjust legend text size
    plot.title = element_text(size = 16, face = "bold")  # Make plot title bold
  )

alpha_div

#save the image
#ggsave("alpha_diversity_plot_CSS.png", alpha_div, width = 10, height = 6, dpi = 300)

## rarefies ##

plot_richness( ps.rarefied, x="DIAGNOSIS", color="DIAGNOSIS", title="rarefied", measures=c("Observed", "Chao1", "ACE", "Simpson", "Shannon")) + 
  geom_boxplot() + 
  theme_bw() 

## stat sign diffrences in alpha?
##https://microbiome.github.io/course_2021_radboud/alpha-diversity.html

## https://rpubs.com/lconteville/713954
richness <- estimate_richness(kostic_crc_CSS)
head(richness)

## for normally distributed shannon index values 
hist(richness$Shannon, breaks =  25)

shapiro_test <- shapiro.test(richness$Shannon)
shapiro_test

    Shapiro-Wilk normality test

data:  richness$Shannon
W = 0.95795, p-value = 3.802e-05

## almost normally distributed but because of the outliers it is not 
anova.sh = aov(richness$Shannon ~ sample_data(kostic_crc_CSS)$DIAGNOSIS)
summary(anova.sh)

                                       Df Sum Sq Mean Sq F value   Pr(>F)    
sample_data(kostic_crc_CSS)$DIAGNOSIS   1   7.53   7.527   17.31 4.97e-05 ***
Residuals                             175  76.10   0.435                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

kruskal.test(richness$Shannon ~ sample_data(kostic_crc_CSS)$DIAGNOSIS)

    Kruskal-Wallis rank sum test

data:  richness$Shannon by sample_data(kostic_crc_CSS)$DIAGNOSIS
Kruskal-Wallis chi-squared = 16.548, df = 1, p-value = 4.743e-05

## because the alpha div is one number per sample we can append this to the metadata of 
## our phyloseq object 

## accounting for the multiple sampling from the same individial
# Load the nlme package

tail(sample_data(kostic_crc_CSS)$ANONYMIZED_NAME)

[1] "MQE9O" "32I9U" "UTWNW" "UTWNW" "BFJMK" "32I9U"

richness_exp <- cbind(richness, DIAGNOSIS = sample_data(kostic_crc_CSS)$DIAGNOSIS, ANONYMIZED_NAME = as.factor(sample_data(kostic_crc_CSS)$ANONYMIZED_NAME))

mixed_model <- lme(Shannon ~ DIAGNOSIS, data = richness_exp, random = ~ 1 | ANONYMIZED_NAME)

summary(mixed_model)

Linear mixed-effects model fit by REML
  Data: richness_exp 

Random effects:
 Formula: ~1 | ANONYMIZED_NAME
        (Intercept)  Residual
StdDev:   0.3081981 0.5851293

Fixed effects:  Shannon ~ DIAGNOSIS 

 Correlation: 
               (Intr)
DIAGNOSISTumor -0.618

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-4.12130368 -0.46019741  0.05461757  0.57266976  2.52865584 

Number of Observations: 177
Number of Groups: 94 

## if you have experimental design where the group you check has more than two levels
## you will have to perform posthoc test in addition to anova and KW

## here is an online tutorial that can help with this
## https://scienceparkstudygroup.github.io/microbiome-lesson/04-alpha-diversity/index.html


## Faith's diversity or PD (Phylogenetic diversity) are based on phylogenetic information for each sample (which we don't have for this data)

# library(biomeUtils)
# data("FuentesIliGutData")
# # reduce size for example
# ps1 <- subset_samples(FuentesIliGutData, ILI == "C")
# ps1 <- prune_taxa(taxa_sums(ps1) > 0, ps1)
# 
# meta_tib <- calculatePD(ps1, justDF = TRUE)
# # check
# meta_tib[c(1, 2, 3), c("PD", "SR")]
# #>                 PD  SR
# #> sample_1  967.3003 259
# #> sample_2 1035.1305 291
# #> sample_3 1189.0248 336

## calculate Pielou’s evenness
data_evenness <- vegan::diversity(data.frame(otu_table(kostic_crc_CSS))) / log(specnumber(data.frame(otu_table(kostic_crc_CSS))))
head(data_evenness)

   304309    469478    208196    358030     16076     35786 
0.7162901 0.6620351 0.5970948       NaN 0.7280863 0.5032583 

8. Estimate and plot the overall dissimilarity (beta diversity) in the microbial composition between the groups

## beta diversity ##
dist_methods = unlist(distanceMethodList)
dist_methods

    UniFrac1     UniFrac2        DPCoA          JSD     vegdist1     vegdist2     vegdist3 
   "unifrac"   "wunifrac"      "dpcoa"        "jsd"  "manhattan"  "euclidean"   "canberra" 
    vegdist4     vegdist5     vegdist6     vegdist7     vegdist8     vegdist9    vegdist10 
      "bray" "kulczynski"    "jaccard"      "gower"   "altGower"   "morisita"       "horn" 
   vegdist11    vegdist12    vegdist13    vegdist14    vegdist15   betadiver1   betadiver2 
 "mountford"       "raup"   "binomial"       "chao"        "cao"          "w"         "-1" 
  betadiver3   betadiver4   betadiver5   betadiver6   betadiver7   betadiver8   betadiver9 
         "c"         "wb"          "r"          "I"          "e"          "t"         "me" 
 betadiver10  betadiver11  betadiver12  betadiver13  betadiver14  betadiver15  betadiver16 
         "j"        "sor"          "m"         "-2"         "co"         "cc"          "g" 
 betadiver17  betadiver18  betadiver19  betadiver20  betadiver21  betadiver22  betadiver23 
        "-3"          "l"         "19"         "hk"        "rlb"        "sim"         "gl" 
 betadiver24        dist1        dist2        dist3   designdist 
         "z"    "maximum"     "binary"  "minkowski"        "ANY" 

## Jaccard and Bray-Curtis dissimilarities are non-phylogenetic measures that consider the presence/absence and abundance of features in samples, respectively, while UniFrac dissimilarities are phylogenetic measures that 
## incorporate information about the evolutionary relatedness of microbial taxa

physeq2.ord <- phyloseq::ordinate(kostic_crc_CSS, "NMDS", "bray")

Square root transformation
Wisconsin double standardization
Run 0 stress 0.2315217 
Run 1 stress 0.2350943 
Run 2 stress 0.2275925 
... New best solution
... Procrustes: rmse 0.04361084  max resid 0.327775 
Run 3 stress 0.2321461 
Run 4 stress 0.2386185 
Run 5 stress 0.2289021 
Run 6 stress 0.226593 
... New best solution
... Procrustes: rmse 0.03804441  max resid 0.2846506 
Run 7 stress 0.2309735 
Run 8 stress 0.2328498 
Run 9 stress 0.2342875 
Run 10 stress 0.2322772 
Run 11 stress 0.2336188 
Run 12 stress 0.2452596 
Run 13 stress 0.2404151 
Run 14 stress 0.2344747 
Run 15 stress 0.2329664 
Run 16 stress 0.2321521 
Run 17 stress 0.2311227 
Run 18 stress 0.2300795 
Run 19 stress 0.2349975 
Run 20 stress 0.2359982 
*** Best solution was not repeated -- monoMDS stopping criteria:
    10: no. of iterations >= maxit
    10: stress ratio > sratmax

## NMDS is non-linear and focuses on preserving rank order of dissimilarities, PCA is linear and focuses on capturing maximum variance, 
## and PCoA can capture both linear and non-linear relationships and different types of distances measures

#shape is also a parameter we can assign metadata variable to 
p = plot_ordination(kostic_crc_CSS, physeq2.ord, type="samples", color="DIAGNOSIS",
                    title="OTUs")

# Add ellipses
beta_div <- p + stat_ellipse(level = 0.95, type = "norm", geom = "polygon", alpha = 0, aes(color = DIAGNOSIS), size = 1) +
  geom_point(size = 4) + # Adjust the size of points+
  theme_bw() +
  guides(size = "none") + # 
  theme(
    text = element_text(size = 12, face = "bold"),  # Make all text bold
    axis.title = element_text(size = 14, face = "bold"),  # Make axis titles bold
    axis.text = element_text(size = 12),  # Adjust axis text size
    legend.title = element_text(size = 14, face = "bold"),  # Make legend title bold
    legend.text = element_text(size = 12),  # Adjust legend text size
    plot.title = element_text(size = 16, face = "bold")  # Make plot title bold
  )

beta_div

#save image 
#ggsave("beta_diversity_plot_CSS.png", beta_div , width = 10, height = 6, dpi = 300)


## stat of the beta diversity 

metadata = data.frame(sample_data(kostic_crc_CSS))
phy_beta = phyloseq::distance(kostic_crc_CSS, method = "bray")
adonis2(phy_beta ~ DIAGNOSIS, data = metadata, perm=999)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999

adonis2(formula = phy_beta ~ DIAGNOSIS, data = metadata, permutations = 999)
           Df SumOfSqs      R2      F Pr(>F)    
DIAGNOSIS   1    1.434 0.02351 4.2126  0.001 ***
Residual  175   59.579 0.97649                  
Total     176   61.013 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## The adonis2 function fits a multivariate linear model to the distance matrix using the grouping variable (DIAGNOSIS) as a predictor. 
## It then assesses the significance of the relationship between the grouping variable and the structure of the distance matrix through permutation testing.
## it assumes homogeneity of dispersion among groups

## A small p-value suggests that the grouping variable significantly explains the variation in the distance matrix,
## indicating differences in community composition or structure between the groups defined by the DIAGNOSIS variable

## betadisper function checks if the within groups dispersion are homogeneous(compositions vary similarly within the group)

bd = betadisper(phy_beta, metadata$'DIAGNOSIS')
anova(bd)

Analysis of Variance Table

Response: Distances
           Df  Sum Sq  Mean Sq F value    Pr(>F)    
Groups      1 0.05661 0.056606  12.841 0.0004393 ***
Residuals 175 0.77141 0.004408                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## perform the analyses correcting for multiple sampling from the same individual

## change the sample id to be a factor 

# Extract sample data from the phyloseq object
sample_data <- sample_data(kostic_crc_CSS)
# write.csv(data.frame(sample_data), "sample_data.csv")

# Modify the variable of interest to be a factor
sample_data$ANONYMIZED_NAME <- as.factor(sample_data$ANONYMIZED_NAME)
levels(sample_data$ANONYMIZED_NAME)

 [1] "32I9U" "38U4V" "41E1K" "59S9W" "5EKFO" "5TA9V" "6G2KB" "82S3M" "A8A34" "BFJMK" "C0095"
[12] "C0112" "C0133" "C0134" "C0149" "C0154" "C0159" "C0186" "C0195" "C0198" "C0206" "C0209"
[23] "C0211" "C0212" "C0225" "C0230" "C0235" "C0240" "C0241" "C0252" "C0256" "C0258" "C0268"
[34] "C0269" "C0270" "C0271" "C0275" "C0282" "C0285" "C0294" "C0306" "C0308" "C0311" "C0314"
[45] "C0315" "C0318" "C0322" "C0325" "C0330" "C0332" "C0333" "C0334" "C0335" "C0341" "C0342"
[56] "C0344" "C0349" "C0355" "C0362" "C0366" "C0371" "C0374" "C0378" "C0388" "C0394" "C0395"
[67] "C0398" "C0399" "C1"    "C10"   "C4"    "C6"    "G2OD5" "G2OD6" "G3UBQ" "GQ6LS" "HZIMM"
[78] "I7ROL" "JIDZE" "KIXFR" "LRQ9W" "MQE9O" "MQETM" "O436F" "OTGGZ" "QFHRS" "R8J9Z" "TV28I"
[89] "UEL2G" "UTWNW" "UZ65X" "XBS5M" "XZ33P" "YOTV6"

# Replace the modified sample data in the phyloseq object
sample_data(kostic_crc_CSS) <- sample_data


## check relationship between individuals 
plot_ordination(kostic_crc_CSS, physeq2.ord, type="samples", color="ANONYMIZED_NAME",
                title="OTUs")+
  guides(color = FALSE) # no legend 

# Extract the NMDS results from physeq2.ord
nmds_results <- scores(physeq2.ord)

# Combine the NMDS results with the sample IDs
plot_data <- cbind(sample_data(kostic_crc_CSS), nmds_results$sites)

# Plot the ordination
p <- ggplot(plot_data, aes(x = NMDS1, y = NMDS2, color = X.SampleID)) +
  geom_point() +
  labs(title = "OTUs") +
  guides(color = FALSE)

p + geom_text(aes(label = ANONYMIZED_NAME), size =2.5, vjust = -0.5)

adonis2(phy_beta ~ metadata$DIAGNOSIS, data = metadata, perm=999, strata = metadata$ANONYMIZED_NAME)

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Blocks:  strata 
Permutation: free
Number of permutations: 999

adonis2(formula = phy_beta ~ metadata$DIAGNOSIS, data = metadata, permutations = 999, strata = metadata$ANONYMIZED_NAME)
                    Df SumOfSqs      R2      F Pr(>F)    
metadata$DIAGNOSIS   1    1.434 0.02351 4.2126  0.001 ***
Residual           175   59.579 0.97649                  
Total              176   61.013 1.00000                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## try weighted unifrac for diagnosis 

# phy_tree(kostic_crc_CSS)
# physeq2.ord1<- ordinate(kostic_crc_CSS, "NMDS", "wunifrac") #or unifrac
# we don't have the phy_tree 
# but if we did, we could have also explored this beta_div measure 

9. Identify DA taxa between groups

## microbiomeMarker is a wrapper for a lot of R packages specifically designed to 
## profile DA taxa from microbial compositional data 

## lefse ##

## KW detect stat sign features
## Wolcoxon and multigrp_strat test for consistency of stat sign across groups of replicates
## LDA analyses are finally performed to identify the exact effect size of the difference 

## CPM is the default norm method for lefse ##

df_kostic_cpm =normalize(kostic_crc, method = "CPM")%>%microbiomeMarker::run_lefse(
                                       wilcoxon_cutoff = 0.05,
                                       norm = "none",
                                       group = "DIAGNOSIS",
                                       kw_cutoff = 0.05,
                                       multigrp_strat = TRUE,
                                       lda_cutoff = 1.5
)

df_kostic_cpm_table  = (marker_table(df_kostic_cpm))
dim(df_kostic_cpm_table)
# 158 DA taxa across taxa levels 
head(df_kostic_cpm_table)

#write the table
#write.csv(data.frame(df_kostic_cpm_table), "df_kostic_cpm_table_DA.csv")

## results are obtained not only for species but for all taxa levels 
## We can't easily correct for co-founders or repeated measures with these analyses 


### prepare the data for Maaslin2
## MaAsLin2 relies on general linear models to accommodate most modern epidemiological study designs, 
## including cross-sectional and longitudinal, and offers a variety of data exploration, normalization, and transformation methods
## it allows correction for co-founders and repeated measures (correcting for intraindividual variability)


kostic_crc_CPM = microbiomeMarker::normalize(kostic_crc, method = "CPM")

crc_CPM_OTU = data.frame(otu_table(kostic_crc_CPM))
dim(crc_CPM_OTU)
head(crc_CPM_OTU)
crc_CPM_TAXA = data.frame((tax_table(kostic_crc_CPM)))
dim(tax_table(kostic_crc_CPM))
head(crc_CPM_TAXA)
crc_CPM_META = data.frame(sample_data(kostic_crc_CPM))
head(crc_CPM_META)
# These are how they should be 
crc_CPM_TAXA = crc_CPM_TAXA %>% dplyr::mutate(Species1 = paste(rownames(crc_CPM_TAXA), Genus, Species, sep = "_"))
crc_CPM_TAXA = crc_CPM_TAXA %>% dplyr::select(Species1)

head(crc_CPM_TAXA)
head(crc_CPM_OTU[, 1:20])

# Add row names as a column in both data frames
crc_CPM_OTU <- rownames_to_column(crc_CPM_OTU, "RowNames")
crc_CPM_TAXA <- rownames_to_column(crc_CPM_TAXA, "RowNames")

# Merge the two data frames by the new column "RowNames"
crc_CPM_OTU_S <- merge(crc_CPM_OTU, crc_CPM_TAXA, by.x = "RowNames", by.y = "RowNames")
head(crc_CPM_OTU_S)
# Optionally, remove the redundant "RowNames" column after merging
crc_CPM_OTU_S <- crc_CPM_OTU_S %>% dplyr::select(-RowNames)
crc_CPM_OTU_S = column_to_rownames(crc_CPM_OTU_S, "Species1")

crc_CPM_OTU_S = data.frame(t(crc_CPM_OTU_S))
head(crc_CPM_OTU_S)
colnames(crc_CPM_OTU_S) <- sub("X", "", colnames(crc_CPM_OTU_S))

crc_CPM_META$ANONYMIZED_NAME = as.factor(crc_CPM_META$ANONYMIZED_NAME)


head(crc_CPM_OTU_S[1:20, 1:20])
head(crc_CPM_META[1:20, 1:20])

## We will need crc_CPM_OTU_S and crc_CPM_META

fit_data <- Maaslin2(
  crc_CPM_OTU_S, crc_CPM_META, 'C://Users//Malina//Desktop//Microbial compositionality/maaslin_results_random_effect/',
  fixed_effects = c('DIAGNOSIS'),
  normalization = "NONE",
  reference = c("DIAGNOSIS,Healthy"),
  random_effects = c("ANONYMIZED_NAME")
  )


## DA analyses of the same dataset with deseq2
## https://www.bioconductor.org/packages/release/bioc/vignettes/phyloseq/inst/doc/phyloseq-mixture-models.html

10. Get co-abundance interaction network for either Tumor or Healthy samples

### will subset the initial CPM norm phyloseq obj to get only the Tumor samples and will
### generate graph from it 
## will subset species that are present at least in 50% of the samples (filtering on prevalence)
## threshold for correlation will be 0.4 (spearman)

kostic_crc_CPM

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 2505 taxa and 177 samples ]
sample_data() Sample Data:       [ 177 samples by 71 sample variables ]
tax_table()   Taxonomy Table:    [ 2505 taxa by 7 taxonomic ranks ]

subset_kostic_crc <- subset_samples(kostic_crc_CPM, DIAGNOSIS == "Tumor")
# Define the prevalence threshold (e.g., taxa present in at least 50% of samples)
prevalence_threshold <- 0.5

# Filter taxa based on prevalence
subset_kostic_crc <- filter_taxa(subset_kostic_crc , function(x) sum(x > 0) / length(x) >= prevalence_threshold, TRUE)

## create OTU input table with species or genus (when NA is present for species) names
crc_CPM_OTU = data.frame(otu_table(subset_kostic_crc))
dim(crc_CPM_OTU)

[1] 70 86

crc_CPM_TAXA = data.frame((tax_table(subset_kostic_crc)))
dim(tax_table(subset_kostic_crc))

[1] 70  7

crc_CPM_META = data.frame(sample_data(subset_kostic_crc))
head(crc_CPM_META)

# These are how they should be 
crc_CPM_TAXA = crc_CPM_TAXA %>% dplyr::mutate(Species1= paste(rownames(crc_CPM_TAXA), Genus, Species, sep = "_"))
crc_CPM_TAXA = crc_CPM_TAXA %>% dplyr::select(Species1)

head(crc_CPM_TAXA)

head(crc_CPM_OTU[, 1:20])

# Add row names as a column in both data frames
crc_CPM_OTU <- rownames_to_column(crc_CPM_OTU, "RowNames")
crc_CPM_TAXA <- rownames_to_column(crc_CPM_TAXA, "RowNames")

# Merge the two data frames by the new column "RowNames"
crc_CPM_OTU_S <- merge(crc_CPM_OTU, crc_CPM_TAXA, by.x = "RowNames", by.y = "RowNames")
head(crc_CPM_OTU_S)

# Optionally, remove the redundant "RowNames" column after merging
crc_CPM_OTU_S <- crc_CPM_OTU_S %>% dplyr::select(-RowNames)
crc_CPM_OTU_S = column_to_rownames(crc_CPM_OTU_S, "Species1")

rownames(crc_CPM_OTU_S) <- sub("X", "", rownames(crc_CPM_OTU_S))
head(crc_CPM_OTU_S[,1:20])

dim(crc_CPM_OTU_S)

[1] 70 86

### build unweighted graph 
correlation_matrix <- cor(data.frame(t(crc_CPM_OTU_S)),method = "spearman")
head(correlation_matrix[,1:5])

                                                        X127309_Bacteroides_NA
X127309_Bacteroides_NA                                              1.00000000
X13825_Dialister_Dialister.pneumosintes                            -0.04489624
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius            -0.12919516
X145149_Veillonella_NA                                              0.18438430
X15054_Granulicatella_Granulicatella.elegans                       -0.26004938
X157566_Alistipes_NA                                                0.33367864
                                                        X13825_Dialister_Dialister.pneumosintes
X127309_Bacteroides_NA                                                              -0.04489624
X13825_Dialister_Dialister.pneumosintes                                              1.00000000
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                              0.27232119
X145149_Veillonella_NA                                                               0.02085916
X15054_Granulicatella_Granulicatella.elegans                                        -0.13407975
X157566_Alistipes_NA                                                                 0.05705999
                                                        X14227_Peptostreptococcus_Peptostreptococcus.anaerobius
X127309_Bacteroides_NA                                                                              -0.12919516
X13825_Dialister_Dialister.pneumosintes                                                              0.27232119
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                                              1.00000000
X145149_Veillonella_NA                                                                               0.18551504
X15054_Granulicatella_Granulicatella.elegans                                                         0.04136805
X157566_Alistipes_NA                                                                                -0.01130812
                                                        X145149_Veillonella_NA
X127309_Bacteroides_NA                                              0.18438430
X13825_Dialister_Dialister.pneumosintes                             0.02085916
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius             0.18551504
X145149_Veillonella_NA                                              1.00000000
X15054_Granulicatella_Granulicatella.elegans                       -0.09595219
X157566_Alistipes_NA                                                0.10021371
                                                        X15054_Granulicatella_Granulicatella.elegans
X127309_Bacteroides_NA                                                                   -0.26004938
X13825_Dialister_Dialister.pneumosintes                                                  -0.13407975
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                                   0.04136805
X145149_Veillonella_NA                                                                   -0.09595219
X15054_Granulicatella_Granulicatella.elegans                                              1.00000000
X157566_Alistipes_NA                                                                     -0.16257067

threshold <- 0.4

# Create an adjacency matrix based on the correlation threshold
adjacency_matrix <- ifelse(abs(correlation_matrix) > threshold, 1, 0)
head(adjacency_matrix[,1:5])

                                                        X127309_Bacteroides_NA
X127309_Bacteroides_NA                                                       1
X13825_Dialister_Dialister.pneumosintes                                      0
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                      0
X145149_Veillonella_NA                                                       0
X15054_Granulicatella_Granulicatella.elegans                                 0
X157566_Alistipes_NA                                                         0
                                                        X13825_Dialister_Dialister.pneumosintes
X127309_Bacteroides_NA                                                                        0
X13825_Dialister_Dialister.pneumosintes                                                       1
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                                       0
X145149_Veillonella_NA                                                                        0
X15054_Granulicatella_Granulicatella.elegans                                                  0
X157566_Alistipes_NA                                                                          0
                                                        X14227_Peptostreptococcus_Peptostreptococcus.anaerobius
X127309_Bacteroides_NA                                                                                        0
X13825_Dialister_Dialister.pneumosintes                                                                       0
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                                                       1
X145149_Veillonella_NA                                                                                        0
X15054_Granulicatella_Granulicatella.elegans                                                                  0
X157566_Alistipes_NA                                                                                          0
                                                        X145149_Veillonella_NA
X127309_Bacteroides_NA                                                       0
X13825_Dialister_Dialister.pneumosintes                                      0
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                      0
X145149_Veillonella_NA                                                       1
X15054_Granulicatella_Granulicatella.elegans                                 0
X157566_Alistipes_NA                                                         0
                                                        X15054_Granulicatella_Granulicatella.elegans
X127309_Bacteroides_NA                                                                             0
X13825_Dialister_Dialister.pneumosintes                                                            0
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                                            0
X145149_Veillonella_NA                                                                             0
X15054_Granulicatella_Granulicatella.elegans                                                       1
X157566_Alistipes_NA                                                                               0

# Remove the diagonal values 
diag(adjacency_matrix) <- 0
head(adjacency_matrix[,1:5])

                                                        X127309_Bacteroides_NA
X127309_Bacteroides_NA                                                       0
X13825_Dialister_Dialister.pneumosintes                                      0
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                      0
X145149_Veillonella_NA                                                       0
X15054_Granulicatella_Granulicatella.elegans                                 0
X157566_Alistipes_NA                                                         0
                                                        X13825_Dialister_Dialister.pneumosintes
X127309_Bacteroides_NA                                                                        0
X13825_Dialister_Dialister.pneumosintes                                                       0
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                                       0
X145149_Veillonella_NA                                                                        0
X15054_Granulicatella_Granulicatella.elegans                                                  0
X157566_Alistipes_NA                                                                          0
                                                        X14227_Peptostreptococcus_Peptostreptococcus.anaerobius
X127309_Bacteroides_NA                                                                                        0
X13825_Dialister_Dialister.pneumosintes                                                                       0
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                                                       0
X145149_Veillonella_NA                                                                                        0
X15054_Granulicatella_Granulicatella.elegans                                                                  0
X157566_Alistipes_NA                                                                                          0
                                                        X145149_Veillonella_NA
X127309_Bacteroides_NA                                                       0
X13825_Dialister_Dialister.pneumosintes                                      0
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                      0
X145149_Veillonella_NA                                                       0
X15054_Granulicatella_Granulicatella.elegans                                 0
X157566_Alistipes_NA                                                         0
                                                        X15054_Granulicatella_Granulicatella.elegans
X127309_Bacteroides_NA                                                                             0
X13825_Dialister_Dialister.pneumosintes                                                            0
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius                                            0
X145149_Veillonella_NA                                                                             0
X15054_Granulicatella_Granulicatella.elegans                                                       0
X157566_Alistipes_NA                                                                               0

graph <- graph_from_adjacency_matrix(adjacency_matrix , 
                                     mode = "undirected")
graph

IGRAPH 99c429d UN-- 70 162 -- 
+ attr: name (v/c)
+ edges from 99c429d (vertex names):
[1] X127309_Bacteroides_NA                 --X203590_NA_NA                                          
[2] X127309_Bacteroides_NA                 --X320768_Oscillospira_NA                                
[3] X127309_Bacteroides_NA                 --X470973_Ruminococcus_Ruminococcus.torques              
[4] X13825_Dialister_Dialister.pneumosintes--X470239_Dialister_Dialister.invisus                    
[5] X13825_Dialister_Dialister.pneumosintes--X470805_NA_Parvimonas.micra                            
[6] X13825_Dialister_Dialister.pneumosintes--X500281_NA_NA                                          
[7] X13825_Dialister_Dialister.pneumosintes--X89770_Peptostreptococcus_Peptostreptococcus.anaerobius
+ ... omitted several edges

#plot(graph)

#Cluster based on edge betweenness
ceb = cluster_edge_betweenness(graph)

#Extract community memberships from the clustering
community_membership <- membership(ceb)
head(community_membership)

                                 X127309_Bacteroides_NA 
                                                      1 
                X13825_Dialister_Dialister.pneumosintes 
                                                      2 
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius 
                                                      2 
                                 X145149_Veillonella_NA 
                                                      3 
           X15054_Granulicatella_Granulicatella.elegans 
                                                      4 
                                   X157566_Alistipes_NA 
                                                      5 

length(ceb)

[1] 22

# Count the number of nodes in each cluster
cluster_sizes <- table(community_membership)
cluster_sizes

community_membership
 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 
 2  9  1  1 25  7  1  3  1  1  2  1  2  1  3  1  1  4  1  1  1  1 

# Add cluster memberships as a vertex attribute
igraph::V(graph)$cluster_membership <- community_membership

# Map cluster memberships to colors using a color palette function
cluster_colors <- rainbow(max(community_membership))
head(cluster_colors)

[1] "#FF0000" "#FF4600" "#FF8B00" "#FFD100" "#E8FF00" "#A2FF00"

# Create a color variable based on cluster memberships
V(graph)$color_variable <- cluster_colors[community_membership]

# vetrex degree
vertex_degrees <- igraph::degree(graph)
head(vertex_degrees)

                                 X127309_Bacteroides_NA 
                                                      3 
                X13825_Dialister_Dialister.pneumosintes 
                                                      4 
X14227_Peptostreptococcus_Peptostreptococcus.anaerobius 
                                                      4 
                                 X145149_Veillonella_NA 
                                                      0 
           X15054_Granulicatella_Granulicatella.elegans 
                                                      0 
                                   X157566_Alistipes_NA 
                                                      6 

# Add vertex degree as a vertex attribute
igraph::V(graph)$degree <- vertex_degrees


vertex.attr = list(
  cluster_membership = V(graph)$cluster_membership,
  name = V(graph)$name,
  degrees = V(graph)$degree,
  color = V(graph)$color_variable)

# Specify the file name for saving the GraphML file
graphml_file <- "T1graph_kostic_CPM_with_metadata_unweighted_undirected.graphml"

# Write the igraph object to GraphML format with metadata
write_graph(
  graph,
  graphml_file,
  format  = "graphml")

# Combine vertex attributes into a data frame
vertex_data <- data.frame(vertex.attr)

head(vertex_data)

write.csv(vertex_data, "T1vertex_data_kostic_CPM_with_metadata_unweighted_undirected.csv", row.names = FALSE)

## create a function for building co-abundance network

get_interaction_network = function(norm_phyloseq_obj, group, graphml_file, vertex_data_file){
  
  subset_kostic_crc <- subset_samples(norm_phyloseq_obj, DIAGNOSIS == paste(group))
  # Define the prevalence threshold (e.g., taxa present in at least 50% of samples)
  prevalence_threshold <- 0.5
  
  # Filter taxa based on prevalence
  subset_kostic_crc <- filter_taxa(subset_kostic_crc , function(x) sum(x > 0) / length(x) >= prevalence_threshold, TRUE)
  
  ## create OTU input table with species or genus (when NA is present for species) names
  crc_CPM_OTU = data.frame(otu_table(subset_kostic_crc))

  crc_CPM_TAXA = data.frame((tax_table(subset_kostic_crc)))
 
  
  crc_CPM_META = data.frame(sample_data(subset_kostic_crc))
  
  # These are how they should be 
  crc_CPM_TAXA = crc_CPM_TAXA %>% dplyr::mutate(Species1= paste(rownames(crc_CPM_TAXA), Genus, Species, sep = "_"))
  crc_CPM_TAXA = crc_CPM_TAXA %>% dplyr::select(Species1)
  
  
  # Add row names as a column in both data frames
  crc_CPM_OTU <- rownames_to_column(crc_CPM_OTU, "RowNames")
  crc_CPM_TAXA <- rownames_to_column(crc_CPM_TAXA, "RowNames")
  
  # Merge the two data frames by the new column "RowNames"
  crc_CPM_OTU_S <- merge(crc_CPM_OTU, crc_CPM_TAXA, by.x = "RowNames", by.y = "RowNames")

  # Optionally, remove the redundant "RowNames" column after merging
  crc_CPM_OTU_S <- crc_CPM_OTU_S %>% dplyr::select(-RowNames)
  crc_CPM_OTU_S = column_to_rownames(crc_CPM_OTU_S, "Species1")
  
  rownames(crc_CPM_OTU_S) <- sub("X", "", rownames(crc_CPM_OTU_S))
  
  ### build unweighted graph 
  correlation_matrix <- cor(data.frame(t(crc_CPM_OTU_S)),method = "spearman")
 
  threshold <- 0.4
  
  # Create an adjacency matrix based on the correlation threshold
  adjacency_matrix <- ifelse(abs(correlation_matrix) > threshold, 1, 0)
 
  
  # Remove the diagonal values 
  diag(adjacency_matrix) <- 0
  
  
  graph <- graph_from_adjacency_matrix(adjacency_matrix , 
                                       mode = "undirected")
  
  
  #Cluster based on edge betweenness
  ceb = cluster_edge_betweenness(graph)
  
  #Extract community memberships from the clustering
  community_membership <- membership(ceb)
  
  # Count the number of nodes in each cluster
  cluster_sizes <- table(community_membership)
  
  # Add cluster memberships as a vertex attribute
  igraph::V(graph)$cluster_membership <- community_membership
  
  # Map cluster memberships to colors using a color palette function
  cluster_colors <- rainbow(max(community_membership))
  
  # Create a color variable based on cluster memberships
  V(graph)$color_variable <- cluster_colors[community_membership]
  
  # vetrex degree
  vertex_degrees <- igraph::degree(graph)
 
  # Add vertex degree as a vertex attribute
  igraph::V(graph)$degree <- vertex_degrees
  
  vertex.attr = list(
    cluster_membership = V(graph)$cluster_membership,
    name = V(graph)$name,
    degrees = V(graph)$degree,
    color = V(graph)$color_variable)
  
  # Specify the file name for saving the GraphML file
  graphml_file <- graphml_file
  
  # Write the igraph object to GraphML format with metadata
  write_graph(
    graph,
    graphml_file,
    format  = "graphml")
  
  # Combine vertex attributes into a data frame
  vertex_data <- data.frame(vertex.attr)
  
  write.csv(vertex_data, vertex_data_file, row.names = FALSE)
  

}

## run the function to obtain the graph file and vertex metadata for the Healthy samples

get_interaction_network(
  norm_phyloseq_obj = kostic_crc_CPM,
  group = "Healthy",
  graphml_file = "H1graph_kostic_CPM_with_metadata_unweighted_undirected.graphml",
  vertex_data_file = "H1vertex_data_kostic_CPM_with_metadata_unweighted_undirected.csv")