## ---- echo=FALSE-------------------------------------------------------------- knitr::opts_chunk$set( cache = FALSE, fig.width = 9, message = FALSE, eval=FALSE, warning = FALSE) knitr::opts_knit$set(global.par = TRUE) ## ----deps--------------------------------------------------------------------- # library(ggplot2) # library(vegan) # library(reshape2) # library(miaSim) # library(philentropy) # library(cluster) ## ----pars--------------------------------------------------------------------- # set.seed(42) # n_species <- 5 # n_resources <- 5 # E <- randomE(n_species, n_resources, mean_consumption = 1, mean_production = 3) # growth_rates <- runif(n_species) # monod_constant <- matrix(rbeta(n_species*n_resources, 10,10), # nrow=n_species, # ncol=n_resources) # t_store <- 50 # n.instances <- 1 # no stochastic process: no need to repeat ## ----funcs-------------------------------------------------------------------- # gradient.df.generator <- function(n_row, n_col, density_row, max_gradient, error_interval){ # list_initial <- list() # dissimilarity.gradient <- seq(from = 0, to = max_gradient, length.out = n_row) # for (i in seq_len(n_row)){ # print(i) # if (i == 1){ # row_temp <- rbeta(n_col, 1, 1/n_col) # col_to_remove <- sample(x = seq_len(n_col), size = n_col-n_col*density_row) # row_temp[col_to_remove] <- 0 # list_initial[[i]] <- row_temp # } else { # while (length(list_initial) < i) { # row_temp <- rbeta(n_col, 1, 1/n_col) # col_to_remove <- sample(x = seq_len(n_col), size = n_col-n_col*density_row) # row_temp[col_to_remove] <- 0 # diff_temp <- abs(vegdist(rbind(list_initial[[1]], row_temp), method = "bray") - dissimilarity.gradient[i]) # if (diff_temp < error_interval) { # list_initial[[i]] <- row_temp # } # } # } # } # dataframe_to_return <- as.data.frame(t(matrix(unlist(list_initial), ncol = n_row))) # return(dataframe_to_return) # } ## ----gen---------------------------------------------------------------------- # n.community <- 5 # you can also try 20 or even 50. # density.community <- 0.8 # set.seed(42) # community.initial.df <- gradient.df.generator(n_row = n.community, # n_col = n_species, # density_row = density.community, # max_gradient = 0.7, # error_interval = 0.1) # dist.community.initial.df <- vegdist(community.initial.df, method = "bray") # community.initial.tse <- TreeSummarizedExperiment(assays=SimpleList(abundances=t(as.matrix(community.initial.df)))) ## ----plots-------------------------------------------------------------------- # makePlot <- function(out_matrix, title = "abundance of species by time", obj = "species", y.label = "x.t"){ # df <- as.data.frame(out_matrix) # dft <- melt(df, id="time") # names(dft)[2] = obj # names(dft)[3] = y.label # lgd = ncol(df)<= 20 # ggplot(dft, aes_string(names(dft)[1], names(dft)[3], col = names(dft)[2])) + # geom_line(show.legend = lgd, lwd=0.5) + # ggtitle(title) + # theme_linedraw() + # theme(plot.title = element_text(hjust = 0.5, size = 14)) # } # makePlotRes <- function(out_matrix, title = "quantity of compounds by time"){ # df <- as.data.frame(out_matrix) # dft <- melt(df, id="time") # names(dft)[2] = "resources" # names(dft)[3] = "S.t" # lgd = ncol(df)<= 20 # ggplot(dft, aes(time, S.t, col = resources)) + # geom_line(show.legend = lgd, lwd=0.5) + # ggtitle(title) + # theme_linedraw() + # theme(plot.title = element_text(hjust = 0.5, size = 14)) # } # makeHeatmap <-function(matrix.A, # title = "Consumption/production matrix", # y.label = 'resources', # x.label = 'species', # midpoint_color = NULL, # lowColor = "red", # midColor = "white", # highColor = "blue"){ # df <- melt(t(matrix.A)) # if (is.null(midpoint_color)) { # midpoint_color <- 0 # } # names(df)<- c("x", "y", "strength") # df$y <- factor(df$y, levels=rev(unique(sort(df$y)))) # fig <- ggplot(df, aes(x,y,fill=strength)) + geom_tile() + coord_equal() + # theme(axis.title = element_blank()) + # scale_fill_gradient2('strength', low = lowColor, # mid = midColor, high = highColor, midpoint = midpoint_color)+ # theme_void() + ggtitle(title) # # if (ncol(matrix.A)<=10 & nrow(matrix.A)<=10){ # fig <- fig + geom_text(aes(label = round(strength, 2))) # } else if (ncol(matrix.A)<=15 & nrow(matrix.A)<=15){ # fig <- fig + geom_text(aes(label = round(strength, 1))) # } else { # fig <- fig # } # # fig <- fig + labs(x = x.label, y = y.label)+ # theme_linedraw() + # theme(plot.title = element_text(hjust = 0.5, size = 14), # axis.text.x = element_text( # angle = 90)) # # if (nrow(matrix.A) >= 20){ # # too many species # fig <- fig + theme( # axis.title.y=element_blank(), # axis.text.y=element_blank(), # axis.ticks.y=element_blank(), # ) # } # if (ncol(matrix.A) >= 20){ # # too many resources # fig <- fig + theme( # axis.title.x=element_blank(), # axis.text.x=element_blank(), # axis.ticks.x=element_blank() # ) # } # fig # } # # makeHeatmap(as.matrix(dist.community.initial.df), # title = "dissimilarity matrix", # x.label = "community.1", # y.label = "community.2") ## ----crm---------------------------------------------------------------------- # crm_params <- list( # n_species = n_species, # n_resources = n_resources, # x0 = NULL, # E = E, # resources = rep(1,n_resources), # monod_constant = monod_constant, # migration_p = 0, # stochastic = FALSE, # t_start = 0, # t_end = 50, # t_step = 1, # t_store = t_store, # growth_rates = growth_rates, # norm=FALSE) ## ----grads-------------------------------------------------------------------- # resourceConcentration <- 10^seq(0,4,1) # 1 to 10000 # n.medium <- 5 # density.medium <- 0.8 # n_species <- 5 # set.seed(42) # resource.initial.df <- gradient.df.generator(n_row = n.medium, # n_col = n_resources, density_row = density.medium, # max_gradient = 0.7, error_interval = 0.1) # crmExample <- simulateConsumerResource( # n_species = n_species, # n_resources = n_resources, # E = E, # x0 = as.numeric(community.initial.df[1,]), # resources = as.numeric(resourceConcentration[3]*resource.initial.df[1,]), # growth_rates = growth_rates, # monod_constant = monod_constant, # stochastic = FALSE, # t_end = 50, # t_step = 1, # t_store = 50, # norm = FALSE) # #makePlot(crmExample$matrix) # #makePlotRes(crmExample$resources) # # ## Generate simulations and store the final community in community.simulation # ## In this step, the final relative abundance table is basisComposition_prop # # set.seed(42) # community.simulation <- list() # counter_i <- 1 # resourceConcentration <- 10^seq(0,4,1) # 1 to 10000 # n.medium <- 5 # for (resConc in resourceConcentration) { # for (medium in seq_len(n.medium)){ # crm_params$resources <- as.numeric(resource.initial.df[medium,]*resConc) # paramx0 <- as.list(as.data.frame(t(community.initial.df))) # crm_param_iter <- list(x0 = paramx0) # print(paste("resConc", resConc, "medium", medium)) # crmMoments <- generateSimulations(model = "simulateConsumerResource", # params_list = crm_params, # param_iter = crm_param_iter, # n_instances = n.instances, # t_end = 50) # # pick community composition at the last time point # community.simulation[[counter_i]] <- as.data.frame(do.call(rbind, lapply(crmMoments, function (x) {assay(x, "counts")[, ncol(x)]}))) # # counter_i <- counter_i + 1 # } # } # basisComposition <- do.call(rbind.data.frame, community.simulation) # rm(counter_i, community.simulation) # basisComposition_prop <- basisComposition / rowSums(basisComposition) # # # ## Make UMAP plots # ## In this step, plot result is stored in umap_CRM_gradient_plot, and # ## this is visualized in different facets. # # resourceConcentration <- 10^seq(0,4,1) # 1 to 10000 # n.medium <- 5 # n.community <- 5 # concentration <- as.factor(rep(resourceConcentration, each = n.medium*n.community)) # medium <- as.factor(rep(seq_len(n.medium), each = n.community ,times = length(resourceConcentration) )) # community <- as.factor(rep(seq_len(n.community), times = length(resourceConcentration)*n.medium)) # # # Visualize with UMAP # # ## Provide the community data as TreeSE object # library(scater) # tse <- TreeSummarizedExperiment( # assays=SimpleList(abundances=t(as.matrix(basisComposition))), # colData=DataFrame(Medium=medium, # Concentration=concentration, # Community=community # ) # ) # ## Add UMAP # tse <- runUMAP(tse, name = "UMAP", exprs_values = "abundances") # ## Plot UMAP # plotReducedDim(tse, "UMAP", colour_by="Medium", shape_by="Concentration") # # # Same for compositional abundance data # library(mia) # ## -- add relative abundances; # tse <- transformSamples(tse, assay.type="abundances", method="relabundance") # tse <- runUMAP(tse, name = "UMAP_compositional", exprs_values = "relabundance") # plotReducedDim(tse, "UMAP_compositional", colour_by="Medium", shape_by="Concentration") # # Finally with communities # umap_CRM_gradient_plot <<- plotReducedDim(tse, "UMAP_compositional", # colour_by="Medium", shape_by="Community", size_by="Concentration") ## ----umapplots, eval=FALSE---------------------------------------------------- # # FIXME: the visual output can be polished later. # print(umap_CRM_gradient_plot) # umap_CRM_gradient_plot + facet_grid(size_by ~ ., labeller = label_both) # umap_CRM_gradient_plot + facet_grid(colour_by ~ size_by, labeller = label_both) # umap_CRM_gradient_plot + facet_grid(shape_by ~ size_by, labeller = label_both) # umap_CRM_gradient_plot + facet_grid(shape_by ~ colour_by, labeller = label_both) ## ----misc, eval=FALSE--------------------------------------------------------- # average_distance <- function(df, res_conc_type, com_type, method = "euclidean"){ # sub_df <- df[df$concentration == res_conc_type & df$community == com_type,] # combines <- combn(sub_df$medium, 2) # distances <- NULL # for (i in seq_len(ncol(combines))) { # distances[i] <- dist(sub_df[combines[,i], c(1, 2)]) # } # return(mean(distances)) # } ## ----df, eval=FALSE----------------------------------------------------------- # distance_saturation_data <- data.frame(concentration = integer(), # community = integer(), # average_distance = numeric()) # # #umap_CRM_coor <- cbind(umap_CRM_coor, concentration, medium, community) # umap_CRM_coor <- data.frame(reducedDim(tse), colData(tse)) # cbind(umap_CRM_coor, concentration, medium, community) # for (res_conc_type in unique(concentration)){ # for (com_type in unique(community)){ # ave_dist <- average_distance(umap_CRM_coor, res_conc_type, com_type) # distance_saturation_data[nrow(distance_saturation_data)+1,] <- # c(res_conc_type, com_type, ave_dist) # } # } # # View(distance_saturation_data) # distance_saturation_data$average_distance <- as.numeric(distance_saturation_data$average_distance) # distance_saturation_data$concentration <- as.factor(distance_saturation_data$concentration) # distance_saturation_data$community <- as.factor(distance_saturation_data$community) # # # distance_saturation_data_plot # p <- ggplot(distance_saturation_data, # aes(concentration, average_distance, # color = community, # group = community)) + # geom_line() + # geom_point() + # scale_shape_manual(values = c(0, 1, 2, 5, 6, 8, 15, 16, 17, 18)) + # labs(x = "Resource concentration", # y = "Average distance between communities in UMAP") + # theme_bw() # # print(p)