Autoencoder play

library(pagoda2)
library(Matrix)
library(ggplot2)
library(cowplot)
library(dplyr)
require(parallel)
library(ggrastr)

p2 <- readRDS("p2.rds")

library(keras)
K <- keras::backend()
#use_session_with_seed(0)
library(tensorflow)
tf$random$set_seed(42)

Autoencoder

# molecular data - use variance-normalized log scale counts
odgenes <- rownames(p2$misc[["varinfo"]])[(order(p2$misc[["varinfo"]]$lp, decreasing = FALSE)[1:3e3])]
x <- p2$counts[,odgenes]
x@x <- x@x * rep(p2$misc[["varinfo"]][colnames(x), "gsf"], diff(x@p))
set.seed(1)
x_test <- x[sample(1:nrow(x),1e3),]
x_train <- x[sample(1:nrow(x),1e3),]
input_size <- ncol(x_train)

model = keras_model_sequential() 
model %>%
  layer_dense(units = 200, input_shape = c(input_size)) %>%  layer_activation_leaky_relu() %>% layer_dropout(rate=0.1) %>%
  layer_dense(units = 2,   name = "bottleneck") %>% layer_activation_leaky_relu() %>%
  layer_dense(units = 15) %>% layer_activation_relu() %>%
  layer_dense(units = 300, activation='hard_sigmoid') %>%
  layer_dense(units = input_size, activation='sigmoid')


model %>% compile(loss='mean_squared_logarithmic_error', metrics='mae', optimizer = 'adam')

model %>% fit(x_train, x_train, shuffle=TRUE, epochs=150, batch_size=200, validation_data=list(x_test, x_test) )

encoder <- keras_model(inputs=model$input,outputs=get_layer(model,'bottleneck')$output)
emb2 <- predict(encoder,x)
rownames(emb2) <- rownames(x);

a2 <- sccore::embeddingPlot(emb2,groups=ann,alpha=0.4,size=0.2, font.size=c(3,4),plot.theme=theme_bw())
a2

Poisson autoencoder

odgenes <- rownames(p2$misc[["varinfo"]])[(order(p2$misc[["varinfo"]]$lp, decreasing = FALSE)[1:2e3])]
x <- p2$misc$rawCounts[,odgenes]
#x@x <- x@x * rep(p2$misc[["varinfo"]][colnames(x), "gsf"], diff(x@p))

x <- x/rowSums(x)*5e3

set.seed(0)
x_test <- x[sample(1:nrow(x),1e3),]
x_train <- x[sample(1:nrow(x),1e3),]
input_size <- ncol(x_train)

model = keras_model_sequential() 
model %>%
  layer_dense(units = 20, input_shape = c(input_size)) %>%  layer_activation_leaky_relu() %>% layer_dropout(rate=0.1) %>%
  layer_dense(units = 2,   name = "bottleneck") %>% layer_activation_leaky_relu() %>% 
  layer_dense(units = 20) %>% layer_activation_relu() %>%
  layer_dense(units = input_size, activation='sigmoid')


model %>% compile(loss='poisson', metrics='mae', optimizer = 'adam')

model %>% fit(x_train, x_train, shuffle=TRUE, epochs=150, batch_size=200, validation_data=list(x_test, x_test) )

encoder <- keras_model(inputs=model$input,outputs=get_layer(model,'bottleneck')$output)
emb3 <- predict(encoder,x)
rownames(emb3) <- rownames(x);

a3 <- sccore::embeddingPlot(emb3,groups=ann,alpha=0.4,size=0.2, font.size=c(3,4),plot.theme=theme_bw())
a3

Learning embedding

# embedding
emb <- p2$embeddings$PCA$tSNE
# normalize to be within (0,1) range
emb <- apply(emb,2,function(x) (x-min(x))/diff(range(x)))
# molecular data - use variance-normalized log scale counts
odgenes <- rownames(p2$misc[["varinfo"]])[(order(p2$misc[["varinfo"]]$lp, decreasing = FALSE)[1:2e3])]
x <- p2$counts[,odgenes]
x@x <- x@x * rep(p2$misc[["varinfo"]][colnames(x), "gsf"], diff(x@p))
x <- x[rownames(x) %in% rownames(emb), ]


set.seed(0)
x_test <- x[sample(1:nrow(x),3e3),]
x_train <- x[sample(1:nrow(x),3e3),]
input_size <- ncol(x_test)

model = keras_model_sequential() 
model %>%
  layer_dense(units = 256, input_shape = c(input_size)) %>%  layer_activation_leaky_relu() %>% layer_dropout(rate=0.1) %>%
  layer_dense(units = 512, activation='relu') %>%
  layer_dense(units = 256, activation='relu') %>%
  layer_dense(units = 2,   name = "bottleneck", activation='sigmoid')



model %>% compile(loss='mean_absolute_error', metrics='mae', optimizer = 'adam')

model %>% fit(x_train, emb[rownames(x_train),], shuffle=TRUE, epochs=200, batch_size=200, validation_data=list(x_test, emb[rownames(x_test),]) )

emb4 <- predict(model,x)
rownames(emb4) <- rownames(x);

vc <- setdiff(rownames(emb),rownames(x_test))
p1 <- sccore::embeddingPlot(emb4[vc,],groups=ann,alpha=0.4,size=0.2, font.size=c(3,4),plot.theme=theme_bw())
p2 <- sccore::embeddingPlot(emb[vc,],groups=ann,alpha=0.4,size=0.2, font.size=c(3,4),plot.theme=theme_bw())
plot_grid(plotlist=list(p1,p2),nrow=1)

saveRDS(list(emb2,emb3,emb4,emb),file='auto.emb.rds')