---

title: "Getting started"

output: rmarkdown::html_vignette

vignette: >

  %\VignetteIndexEntry{Getting started}

  %\VignetteEngine{knitr::rmarkdown}

  %\VignetteEncoding{UTF-8}

---

```{r, echo=FALSE, warning=FALSE, message=FALSE}

if (!reticulate::py_module_available("tensorflow")) {

  knitr::opts_chunk$set(eval = FALSE)

} else {

  knitr::opts_chunk$set(eval = TRUE)

}

```

```{r, message=FALSE}

library(deepG)

library(keras)

library(magrittr)

library(ggplot2)

```

```{r, echo=FALSE, warning=FALSE, message=FALSE}

options(rmarkdown.html_vignette.check_title = FALSE)

```

```{css, echo=FALSE}

mark.in {

  background-color: CornflowerBlue;

}

mark.out {

  background-color: IndianRed;

}

```

## Introduction 

The goal of the deepG package is to speed up the development of bioinformatical tools for sequence classification, homology detection and other bioinformatical tasks. 

The package offers several functions for

+ Data (pre-) processing

+ Deep learning architectures

+ Model training 

+ Model evaluation

+ Visualizing training progress

### Create dummy data

We create two simple dummy training and validation data sets. Both consist of random <tt>ACGT</tt> sequences but the first category has 

a probability of 40% each for drawing <tt>G</tt> or <tt>C</tt> and the second has equal probability for each nucleotide (first category has around 80% <tt>GC</tt> content and second one around 50%).   

```{r warning = FALSE}

set.seed(123)

vocabulary <- c("A", "C", "G", "T")

data_type <- c("train_1", "train_2", "val_1", "val_2")

for (i in 1:length(data_type)) {

  temp_file <- tempfile()

  assign(paste0(data_type[i], "_dir"), temp_file)

  dir.create(temp_file)

  if (i %% 2 == 1) {

    header <- "high_gc"

    prob <- c(0.1, 0.4, 0.4, 0.1)

  } else {

    header <- "equal_dist"

    prob <- rep(0.25, 4)

  }

  fasta_name_start <- paste0(header, "_", data_type[i], "file")

  create_dummy_data(file_path = temp_file,

                    num_files = 1,

                    seq_length = 10000, 

                    num_seq = 1,

                    header = header,

                    prob = prob,

                    fasta_name_start = fasta_name_start,

                    vocabulary = vocabulary)

}

```

### Training

We can now train a model that can differentiate between the two categories. First, we can create our network architecture. We take an input size of 50 nucleotides. The model has one lstm layer with 16 cells and two dense layers with 8 and 2 neurons.

```{r warning = FALSE}

maxlen <- 50

model <- create_model_lstm_cnn(maxlen = maxlen,

                               layer_lstm = 16,

                               layer_dense = c(8, 2))

```

Next we can train the model using the `train_model` function. Function will internally build a data generator for training. 

```{r warning = FALSE}

hist <- train_model(model,

                    train_type = "label_folder",

                    run_name = "gc_model_1",

                    path = c(train_1_dir, train_2_dir),

                    path_val = c(val_1_dir, val_2_dir),

                    epochs = 4,

                    steps_per_epoch = 25, # one epoch = 25 batches

                    batch_size = 64,

                    step = 50, # take a sample every 50 nt

                    vocabulary_label = c("high_gc", "equal_dist"))

plot(hist)

```

### Evaluation 

We can now evaluate the trained model on all the validation data

```{r warning = FALSE}

eval <- evaluate_model(path_input = c(val_1_dir, val_2_dir),

                       model = model,

                       batch_size = 100,

                       step = 25, # take a sample every 25 nt 

                       vocabulary_label = list(c("high_gc", "equal_dist")),

                       mode = "label_folder",

                       evaluate_all_files = TRUE,

                       verbose = FALSE,

                       auc = TRUE,

                       auprc = TRUE)

eval

```

We can check where our model made mistakes for the sequence with high GC content.

```{r warning = FALSE}

high_gc_file <- microseq::readFasta(list.files(val_1_dir, full.names = TRUE)[1])

high_gc_seq <- high_gc_file$Sequence

pred_high_gc <- predict_model(model = model, 

                              sequence = high_gc_seq,

                              filename = NULL, 

                              step = 25,

                              batch_size = 512,

                              verbose = TRUE,

                              return_states = TRUE,

                              mode = "label")

pred_df <- cbind(pred_high_gc$states, pred_high_gc$sample_end_position) %>% 

  as.data.frame()

names(pred_df) <- c("high_gc_conf", "equal_dist_conf", "sample_end_position")

head(pred_df)

wrong_pred <- pred_df %>% dplyr::filter(high_gc_conf < 0.5)

wrong_pred

if (nrow(wrong_pred) == 0) {

  print("All predictions for high GC content class correct")

} else {

  # extract samples where model was wrong

  wrong_pred_seq <- vector("character", nrow(wrong_pred))

  for (i in 1:length(wrong_pred_seq)) {

    sample_end <- wrong_pred$sample_end_position[i]

    sample_start <- sample_end - maxlen + 1

    wrong_pred_seq[i] <- substr(high_gc_seq, sample_start, sample_end)

  }

  wrong_pred_seq

}

```

We can check the nucleotide distribution of those sequences

```{r warning = FALSE}

l <- list()

for (i in 1:length(wrong_pred_seq)) {

  l[[i]] <- stringr::str_split(wrong_pred_seq[i], "") %>% table() %>% prop.table() %>% t() %>% as.matrix()

}

dist_matrix <- do.call(rbind, l)

dist_matrix

df <- data.frame(distribution = as.vector(dist_matrix),

                 nt = factor(rep(vocabulary, each = nrow(dist_matrix))),

                 sample_id = rep(1:nrow(dist_matrix), 4))

ggplot(df, aes(fill=nt, y=distribution, x=nt)) + 

    geom_bar(position="dodge", stat="identity")  + facet_wrap(~sample_id)

```

Finally, we may want to aggregate all predictions, we made for the sequence.

We can do this using the `summarize_states` function. The function returns 

the mean confidence, the maximum prediction and the vote percentages (percentage of predictions per class).

```{r warning = FALSE}

label_names <- c("high_gc", "equal_dist")

pred_summary <- summarize_states(label_names = label_names, df = pred_df[, 1:2])

print(pred_summary)

```