library(htrsim)
library(tidyverse)
library(reshape2)
library(sva)

## Import & reshape table counts
fn = system.file("extdata/", "public_tablCnts.tsv", package = "htrsim")
tabl_cnts <- read.table(file = fn, header = TRUE)
rownames(tabl_cnts) <- tabl_cnts$gene_id
tabl_cnts <- tabl_cnts %>% select(-gene_id)##suppr colonne GeneID
tabl_cnts <- tabl_cnts %>% select(-gene_name) ##suppr colonne GeneName

## import design of bioProject
fn = system.file("extdata/", "public_bioDesign.csv", package = "htrsim")
bioDesign <- read.table(file = fn, header = T, sep = ';')


dds = run.deseq(tabl_cnts = tabl_cnts, bioDesign = bioDesign)


########## SVA
dat  <- DESeq2::counts(dds, normalized = TRUE)
idx  <- rowMeans(dat) > 1
dat  <- dat[idx, ]
mod <- model.matrix(~genotype + env + genotype:env, bioDesign)
#mod  <- model.matrix(~ dex, colData(dds))
mod0 <- model.matrix(~   1, bioDesign)
svseq <- svaseq(dat, mod, mod0, n.sv = 2)


plot(svseq$sv)


par(mfrow = c(2, 1), mar = c(3,5,3,1))
for (i in 1:2) {
  stripchart(svseq$sv[, i] ~ dds$genotype, vertical = TRUE, main = paste0("SV", i))
  abline(h = 0)
}


## PCA
data2PCA <- dat %>% data.frame() %>%
        mutate(gene_id = rownames(.)) %>% remove_rownames() %>%
        reshape2::melt(.,id = c('gene_id'), variable.name = "sample_id", value.name = "count_norm") %>%
        dcast(., sample_id ~ gene_id)

rownames(data2PCA) = data2PCA$sample_id
data2PCA = data2PCA %>% select(-sample_id)

## PCA processed
pca.obj <- prcomp( data2PCA , scale. = FALSE)
summary(pca.obj)
plot(pca.obj)
## Annotation of the pca object
sample_id = data2PCA %>% rownames()
environment = data2PCA %>% rownames() %>% str_split(pattern = "_", simplify = TRUE)%>% .[,2]
genotype = data2PCA %>% rownames() %>% str_split(pattern = "_", simplify = TRUE)%>% .[,1]

dtp <- data.frame( 'sample_id' = sample_id ,
                   'genotype' = genotype,
                   'environment' = environment,
                   pca.obj$x[,1:2])

## Plot
P <- ggplot(data = dtp) +
  geom_point(aes(x = PC1, y = PC2,
                 col = genotype, ### Possible choice: period, compartments, day, lumen_mucus, ...
                 shape = environment,  ### Possible choice: period, compartments, day, lumen_mucus, ...
                 text = paste("ID:", sample_id))) +
  theme_minimal()
ggplotly(P,  tooltip = c("text"))
P