<jupyter_start><jupyter_text>t-SNE for species vectors<jupyter_code># Housekeeping

library(caret)
library(gridExtra)
library(reshape2)
library(Rtsne)
library(tidyr)
# Read in data

species_composition = read.table("../../../data/amplicon/species_composition_relative_abundance.txt", 
                                 sep = "\t",
                                 header = T,
                                 row.names = 1)

metadata = read.table("../../../data/amplicon/metadata.txt",
                      sep = "\t",
                      header = T,
                      row.names = 1)

head(species_composition)
head(metadata)
# Extract regime shift data without predation

x = metadata$Experiment == "RegimeShift" & # keep only regime shift data
    metadata$Predation == 0 & # exclude predation
    metadata$Immigration != "stock" # exclude initial strain mix data

# Subset

species_composition = species_composition[x,]
metadata = metadata[x,]

species_composition$TREATMENT = rownames(species_composition)
metadata$TREATMENT = rownames(metadata)

dim(species_composition)

df = merge(species_composition, metadata, by = "TREATMENT")
dim(df)
df2 = melt(df, id.vars = c("TREATMENT", "Replicate", "Streptomycin", "Predation", "Immigration", "Time_point", "Experiment"))
colnames(df2) = c("TREATMENT", "Replicate", "Streptomycin", "Predation", "Immigration", "Time_point", "Experiment", "Species", "Abundance")
df2$TREATMENT = paste(df2$Streptomycin, df2$Immigration, df2$Replicate, sep = "_")
head(df2)
# Create abundance vector matrix

df3 = dcast(data = df2,formula = TREATMENT + Species~Time_point,fun.aggregate = sum,value.var = "Abundance")
rownames(df3) = paste(df3$Species, df3$TREATMENT, sep = "_")
m1 = as.matrix(df3[,-c(1:2)])
head(m1)
# Remove empty rows

dim(m1)
m1 = m1[rowSums(m1) > 0.03,]
dim(m1)
head(m1)
# Remove species with 3 or less data points

df4 = as.data.frame(m1)
df4$labels = rownames(df4)
df4 = separate(df4, labels, c("S1", "S2", "S3", "S4", "Streptomycin", "Immigration", "Replicate"), sep = "_")
df4$Species = paste(df4$S1, df4$S2, df4$S3, df4$S4, sep = " ")
df4 = df4[,-c(4:7)]
df5 = as.data.frame(table(df4$Species))
colnames(df5) = c("Species", "Count")
df6 = merge(df4, df5, all.x = T)
df6 = df6[order(df6$Streptomycin, df6$Immigration, df6$Replicate),]
df4$Species == df6$Species
x = df6$Count >= 3
m1 = m1[x,]
dim(m1)
# Rtsne function may take some minutes to complete...

set.seed(9)
tsne_model_1 = Rtsne(m1, check_duplicates=FALSE, pca=TRUE, perplexity=30, theta=0.5, dims=2)
# Getting the two dimension matrix

d_tsne_1 = as.data.frame(tsne_model_1$Y)
# Plotting the results without clustering

ggplot(d_tsne_1, aes(x=V1, y=V2)) +
  geom_point(size=0.25) +
  guides(colour=guide_legend(override.aes=list(size=6))) +
  xlab("") + ylab("") +
  ggtitle("t-SNE") +
  theme_light(base_size=20) +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank()) +
  scale_colour_brewer(palette = "Set2")
# Setting the cluster model

# Next piece of code will create the k-means and hierarchical cluster models.
# To then assign the cluster number (1-7) to which each input case belongs.

## keeping original data
d_tsne_1_original=d_tsne_1

## Creating k-means clustering model, and assigning the result to the data used to create the tsne
fit_cluster_kmeans=kmeans(scale(d_tsne_1), 6)
d_tsne_1_original$cl_kmeans = factor(fit_cluster_kmeans$cluster)

## Creating hierarchical cluster model, and assigning the result to the data used to create the tsne
fit_cluster_hierarchical=hclust(dist(scale(d_tsne_1)))

## setting 3 clusters as output
d_tsne_1_original$cl_hierarchical = factor(cutree(fit_cluster_hierarchical, k=6))
# Plotting the cluster models onto t-SNE output
# Now time to plot the result of each cluster model, based on the t-SNE map.

plot_cluster=function(data, var_cluster, palette)
{
  ggplot(data, aes_string(x="V1", y="V2", color=var_cluster)) +
  geom_point(size=0.25) +
  guides(colour=guide_legend(override.aes=list(size=6))) +
  xlab("") + ylab("") +
  ggtitle("") +
  theme_light(base_size=20) +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        legend.direction = "horizontal", 
        legend.position = "bottom",
        legend.box = "horizontal") + 
    scale_colour_brewer(palette = palette) 
}


plot_k=plot_cluster(d_tsne_1_original, "cl_kmeans", "Accent")
plot_h=plot_cluster(d_tsne_1_original, "cl_hierarchical", "Set1")

## and finally: putting the plots side by side with gridExtra lib...
grid.arrange(plot_k, plot_h,  ncol=2)
dim(d_tsne_1)
d_tsne_1$labels = rownames(m1)
head(d_tsne_1)
d_tsne_1 = separate(d_tsne_1, labels, c("S1", "S2", "S3", "S4", "Streptomycin", "Immigration", "Replicate"), sep = "_")
d_tsne_1$Species = paste(d_tsne_1$S1, d_tsne_1$S2, d_tsne_1$S3, d_tsne_1$S4, sep = " ")
d_tsne_1 = d_tsne_1[,-c(3:6)]
head(d_tsne_1)
# Plotting the results by adding treatment information

ggplot(d_tsne_1, aes(x=V1, y=V2, colour = Species)) +
  geom_point(size=0.25) +
  guides(colour=guide_legend(override.aes=list(size=6))) +
  xlab("") + ylab("") +
  ggtitle("t-SNE") +
  theme_light(base_size=20) +
  theme(axis.text.x=element_blank(),
        axis.text.y=element_blank(),
        legend.direction = "horizontal", 
        legend.position = "bottom",
        legend.box = "horizontal") + 
    scale_colour_brewer(palette = "Set2")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14625092798812175, 'fraction_numerical': 0.027654046028210838, 'mean_word_length': 4.625260960334029, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '477885', 'license': 'MIT', 'n_tokens_mistral': 2129, 'n_tokens_neox': 2022, 'n_words': 414, 'path': 'src/notebooks/amplicon_data/.ipynb_checkpoints/tSNE_species_vectors-checkpoint.ipynb', 'repo_name': 'johannescairns/repeatable_dynamics', 'size': '5388'}