Presentation.Rmd

---
title: "Presentation"
author: "s210136, s200177, s202754"
date: "5/10/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Define functions ------------------------------------------------------------

```{r}
source(file = "R/99_project_functions.R")
```

```{r libraries, include=FALSE}
library("tidyverse")
library("tidyseurat")
library("viridis")
library("broom")
library("nnet")
library("patchwork")
library("dendextend")
library("heatmapply")
library("ggdendro")
library("shiny")
library("shinythemes")
```

# Loading data of all patients ------------------------------------------------

```{r patient_load}
patient_1 <- read_rds(file = "Data/_raw/GSM4058963_025I.dgecounts.rds")
```
Get gene identifiers from first patient

```{r gene_extract}
Genes <- patient_1 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>%   
  pluck("Dimnames") %>% 
  pluck(1) %>% 
  as_tibble()

write_csv(Genes,
          file = "Data/Genes_1")
rm(Genes)
```

We downsample in load as we cannot work with such big data on PCs
```{r downsample}
patient_1 <- patient_1 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>% 
  as_tibble() %>% 
  select(c(1:2000))
```

Repeat operations for 5 remaining subjects (not shown)
```{r, echo=FALSE}
patient_2 <- read_rds(file = "Data/_raw/GSM4058907_465C.dgecounts.rds")

Genes <- patient_2 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>%   
  pluck("Dimnames") %>% 
  pluck(1) %>% 
  as_tibble()

write_csv(Genes,
          file = "Data/Genes_2")
rm(Genes)

patient_2 <- patient_2 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>% 
  as_tibble() %>% 
  select(c(1:2000))

patient_3 <- read_rds(file = "Data/_raw/GSM4058921_003C.dgecounts.rds")

Genes <- patient_3 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>%   
  pluck("Dimnames") %>% 
  pluck(1) %>% 
  as_tibble()

write_csv(Genes,
          file = "Data/Genes_3")
rm(Genes)

patient_3 <- patient_3 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>% 
  as_tibble() %>% 
  select(c(1:2000))

patient_4 <- read_rds(file = "Data/_raw/GSM4058936_207CO.dgecounts.rds")

Genes <- patient_4 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>%   
  pluck("Dimnames") %>% 
  pluck(1) %>% 
  as_tibble()

write_csv(Genes,
          file = "Data/Genes_4")
rm(Genes)

patient_4 <- patient_4 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>% 
  as_tibble() %>% 
  select(c(1:2000))

patient_5 <- read_rds(file = "Data/_raw/GSM4058944_056CO.dgecounts.rds")

Genes <- patient_5 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>%   
  pluck("Dimnames") %>% 
  pluck(1) %>% 
  as_tibble()

write_csv(Genes,
          file = "Data/Genes_5")
rm(Genes)

patient_5 <- patient_5 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>% 
  as_tibble() %>% 
  select(c(1:2000))

patient_6 <- read_rds(file = "Data/_raw/GSM4058977_177I.dgecounts.rds")

Genes <- patient_6 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>%   
  pluck("Dimnames") %>% 
  pluck(1) %>% 
  as_tibble()

write_csv(Genes,
          file = "Data/Genes_6")
rm(Genes)

patient_6 <- patient_6 %>% 
  pluck("umicount") %>% 
  pluck("exon") %>% 
  pluck("all") %>% 
  as_tibble() %>% 
  select(c(1:2000))

```

After all operations have been carried out we write all the data to seperate files
```{r writing}
# Patient_1
write.csv(patient_1,
          file = "Data/01_patient_1I.csv.gz")

# Patient_2
write.csv(patient_2,
          file = "Data/01_patient_2C.csv.gz")

# Patient_3
write.csv(patient_3,
          file = "Data/01_patient_3C.csv.gz")

# Patient_4
write.csv(patient_4,
          file = "Data/01_patient_4CO.csv.gz")

# Patient_5
write.csv(patient_5,
          file = "Data/01_patient_5CO.csv.gz")

# Patient_6
write.csv(patient_6,
          file = "Data/01_patient_6I.csv.gz")
```


Cleaning the data:
Load data --------------------------------------------------------------------

```{r}
# Patient 1
patient_1 <- read_csv(file = "Data/01_patient_1I.csv.gz")

# Patient_2
patient_2 <- read_csv(file = "Data/01_patient_2C.csv.gz")

# Patient_3
patient_3 <- read_csv(file = "Data/01_patient_3C.csv.gz")

# Patient_4
patient_4 <- read_csv(file = "Data/01_patient_4CO.csv.gz")

# Patient_5
patient_5 <- read_csv(file = "Data/01_patient_5CO.csv.gz")

# Patient_6
patient_6 <- read_csv(file = "Data/01_patient_6I.csv.gz")

```


Combining patients -----------------------------------------------------------

Wrangle data -----------------------------------------------------------------

 ------------------- Getting gene names and cell bar codes -------------------

Loading metadata and getting cell types --------------------------------------
Metadata supplied by the authors of the paper("ids.csv")

```{r}
gene_mapping <- read_csv("Data/_raw/ids.csv",
                         col_names=c("value","Gene")) %>%
  tibble()

```

Transpose so that genes are columns and rows are cells

```{r}
patient_1 <- patient_1 %>% 
  rownames_to_column() %>%  
  pivot_longer(-rowname) %>% 
  pivot_wider(names_from = rowname,
              values_from = value)

Genes <- read_csv("Data/Genes_1")
Genes <- Genes %>%
  left_join(gene_mapping) %>%
  pluck("Gene")

```


Rename genes and make nice name for column with barcodes

```{r}
colnames(patient_1) <- c("Cell_Barcode",
                         Genes)

```


-------------- now do the same for 5 other patients --------------


```{r}
# Patient_2
patient_2 <- patient_2 %>%
  rownames_to_column() %>%  
  pivot_longer(-rowname) %>% 
  pivot_wider(names_from = rowname,
              values_from = value)

Genes <- read_csv("Data/Genes_2")
Genes <- Genes %>% 
  left_join(gene_mapping) %>% 
  pluck("Gene")

colnames(patient_2) <- c("Cell_Barcode",
                         Genes)
# Patient_3
patient_3 <- patient_3 %>% 
  rownames_to_column() %>%  
  pivot_longer(-rowname) %>% 
  pivot_wider(names_from = rowname,
              values_from = value)

Genes <- read_csv("Data/Genes_3")
Genes <-Genes %>% 
  left_join(gene_mapping) %>% 
  pluck("Gene")

colnames(patient_3) <- c("Cell_Barcode",
                         Genes)

#Patient_4               
patient_4 <- patient_4 %>% 
  rownames_to_column() %>%  
  pivot_longer(-rowname) %>% 
  pivot_wider(names_from = rowname,
              values_from = value)

Genes <- read_csv("Data/Genes_4")
Genes <- Genes %>% 
  left_join(gene_mapping) %>% 
  pluck("Gene")

colnames(patient_4) <- c("Cell_Barcode",
                         Genes)

#Patient_5               
patient_5 <- patient_5 %>% 
  rownames_to_column() %>%  
  pivot_longer(-rowname) %>% 
  pivot_wider(names_from = rowname,
              values_from = value)

Genes <- read_csv("Data/Genes_5")
Genes <- Genes %>% 
  left_join(gene_mapping) %>% 
  pluck("Gene")

colnames(patient_5) <- c("Cell_Barcode",
                         Genes)

#Patient_6                
patient_6 <- patient_6 %>% 
  rownames_to_column() %>%  
  pivot_longer(-rowname) %>% 
  pivot_wider(names_from = rowname, 
              values_from = value)

Genes <- read_csv("Data/Genes_6")
Genes <- Genes %>% 
  left_join(gene_mapping) %>% 
  pluck("Gene")

colnames(patient_6) <- c("Cell_Barcode",
                         Genes)

```




 Write data -------------------------------------------------------------------
 
 
```{r}
# Patient_1
write.csv(patient_1,
          file = "Data/02_patient_1I.csv.gz")

# Patient_2
write.csv(patient_2,
          file = "Data/02_patient_2C.csv.gz")

# Patient_3
write.csv(patient_3,
          file = "Data/02_patient_3C.csv.gz")

# Patient_4
write.csv(patient_4,
          file = "Data/02_patient_4CO.csv.gz")

# Patient_5
write.csv(patient_5,
          file = "Data/02_patient_5CO.csv.gz")

# Patient_6
write.csv(patient_6,
          file = "Data/02_patient_6I.csv.gz")
```
 
 


Remove Data ------------------------------------------------------------------

```{r}
rm(patient_1,
   patient_2,
   patient_3,
   patient_4,
   patient_5,
   patient_6,
   Genes,
   gene_mapping)

```



Augmenting the data:
Load data --------------------------------------------------------------------

```{r}
# Patient_1
patient_1 <- read_csv(file = "Data/02_patient_1I.csv.gz")

patient_1 <- patient_1 %>% 
  mutate(Patient_ID = "025I",
         .before = Cell_Barcode) %>% 
  select(Patient_ID:ncol(patient_1)) %>% 
  filter(Cell_Barcode != "X1")

# Patient_2                ---------------------------
patient_2 <- read_csv(file = "Data/02_patient_2C.csv.gz")

patient_2 <- patient_2 %>% 
  mutate(Patient_ID = "465C",
         .before = Cell_Barcode) %>% 
  select(Patient_ID:ncol(patient_2)) %>% 
  filter(Cell_Barcode != "X1")

# Patient_3                ---------------------------
patient_3 <- read_csv(file = "Data/02_patient_3C.csv.gz")

patient_3 <- patient_3 %>% 
  mutate(Patient_ID = "003C",
         .before = Cell_Barcode) %>% 
  select(Patient_ID:ncol(patient_3)) %>% 
  filter(Cell_Barcode != "X1")

# Patient_4                ---------------------------
patient_4 <- read_csv(file = "Data/02_patient_4CO.csv.gz")

patient_4 <- patient_4 %>% 
  mutate(Patient_ID = "207CO",
         .before = Cell_Barcode) %>% 
  select(Patient_ID:ncol(patient_4)) %>% 
  filter(Cell_Barcode != "X1")

# Patient_5                ---------------------------
patient_5 <- read_csv(file = "Data/02_patient_5CO.csv.gz")

patient_5 <- patient_5 %>% 
  mutate(Patient_ID = "056CO",
         .before = Cell_Barcode)%>% 
  select(Patient_ID:ncol(patient_5)) %>% 
  filter(Cell_Barcode != "X1")

# Patient_6                ---------------------------
patient_6 <- read_csv(file = "Data/02_patient_6I.csv.gz")

patient_6 <- patient_6 %>% 
  mutate(Patient_ID = "177I",
         .before = Cell_Barcode) %>% 
  select(Patient_ID:ncol(patient_6)) %>% 
  filter(Cell_Barcode != "X1")

```


Combining patients------------------------------------------------------------


```{r}
data <- bind_rows(
  patient_1,
  patient_2,
  patient_3,
  patient_4,
  patient_5,
  patient_6)
```

Replacing NAs (gene not present in the patient but in another) by 0

```{r}
data <- data %>%
  mutate(across(everything(),
           ~replace_na(.x,
                       0))) %>% 
  select(-starts_with("ENSG"))

```


Clearing environment

```{r}
rm(patient_1, 
   patient_2,
   patient_3,
   patient_4,
   patient_5,
   patient_6)
```


Loading metadata -------------------------------------------------------------

```{r}
metadata <- read_tsv(
  "Data/_raw/GSE136831_AllCells.Samples.CellType.MetadataTable.txt") %>% 
  tibble() %>% 
  select(CellBarcode_Identity,
         nUMI,
         nGene,
         CellType_Category,
         Subclass_Cell_Identity,
         Subject_Identity) %>% 
  rename(Patient_ID = Subject_Identity,
         Cell_Barcode = CellBarcode_Identity) %>% 
  filter(Patient_ID == "025I" |
           Patient_ID == "465C" |
           Patient_ID == "003C" |
           Patient_ID == "207CO"|
           Patient_ID == "056CO"|
           Patient_ID == "177I") %>% 
  mutate(Cell_Barcode = str_replace_all(Cell_Barcode,
                                        "(.+_)(.+)",
                                        "\\2"))
```


Join metadata and patient data------------------------------------------------

```{r}
data <- left_join(data,
                  metadata,
                  by = c("Patient_ID",
                         "Cell_Barcode")) %>% 
  relocate("nGene",
           "nUMI",
           "CellType_Category",
           "Subclass_Cell_Identity",
           .after = "Cell_Barcode") 
```


Introduce group label and make that and Patient_ID factors


```{r}
data <-
  data %>% 
  mutate(group = 
           factor(
             case_when(
             str_detect(Patient_ID,"I") ~ "IPF",
             str_detect(Patient_ID,"CO") ~ "COPD",
             TRUE ~ "Control")),         
         .after = Patient_ID,
         Patient_ID = factor(Patient_ID))

rm(metadata)

```



Now there are quite some NAs as the metadata apparently contains only cells 
filtered by the authors.
We will filter further to second-check this filtering. Thus if it works well, 
the NAs will mostly be gone

Wrangle data -----------------------------------------------------------------

Combine the patients' gene expression to the metadata data of the cells.
Run patient_slicer, meta_slicer, combiner and binder

If there are cells with fewer than 2000 transcripts recorded, these cells are
filtered out from the 10000 starting cell count per patient

```{r}
data <- data %>% 
  filter(nUMI > 2000)
```


Filtering out the cells where the transcripts of the mitochondrial genes 
represent more than 20% of the sum of the transcripts for a cell
```{r}
mt_selection <- select(data, 
                       starts_with("MT"))

mt_sum <- mt_selection %>% 
  mutate(mito_sum = rowSums(mt_selection))

data <- data %>%
  mutate(select(mt_sum,
                mito_sum),
         .after = Cell_Barcode) %>%
  filter(mito_sum / nUMI < 0.2)
```


Now check if there are any NAs remaining

```{r}
data %>% 
  filter(
  across(.cols = everything(),
         .fns = ~ is.na(.x)))

```



Write data -------------------------------------------------------------------

```{r}
write_csv(data, 
          file = "Data/03_data.csv")
```


Remove Data ------------------------------------------------------------------

```{r}
rm(data,
   mt_selection,
   mt_sum)
```


## Load data ------------------------------------------------------------------
```{r}
data <- read_csv("Data/03_data.csv")
```

## Model and Plot Data --------------------------------------------------------
Making some summary statistics
First lets reorder the data according to groups
How many cells are left of each patient
```{r}
cells_per_group_and_patient <- data %>% 
  group_by(group,
           Patient_ID) %>% 
  summarise(total = n())

cells_count_plot <- cells_per_group_and_patient %>%
  ggplot(aes(x = group,
             y = total)) +
  geom_col(aes(fill = fct_reorder2(as.factor(Patient_ID),
                                   total,
                                   group))) +
  scale_fill_viridis_d("Patients") +
  theme_minimal() +
  ylab("count") +
  labs(title = "Cell Counts per Group and Patient after Filtering")
```

How many cells are there of each type in each group
```{r}
sum_per_group <- data %>% 
    group_by(group) %>% 
    summarise(total_count = n())
  
cells_type_per_group <- data %>% 
  group_by(group,
           CellType_Category) %>% 
  summarise(cell_count = n()) %>% 
  pivot_wider(names_from = CellType_Category,
              values_from = cell_count) %>% 
  ungroup() %>% 
  left_join(sum_per_group) %>% 
  pivot_longer(Endothelial:Stromal,
               names_to = "CellType_Category",
               values_to = "cell_count")

```

How are cell types distributed
```{r}
cells_type_plot <- cells_type_per_group %>% 
    ggplot(aes(x = group,
               y = cell_count / total_count)) +
    geom_col(aes(fill = fct_reorder2(as.factor(CellType_Category),
                                     cell_count,
                                     group))) +
    scale_fill_viridis_d("Type", 
                         alpha = 0.6) +
    theme_minimal() +
    coord_flip() +
    ylab("relative count") +
    labs(title = "Distribution of Cell Types per Group") +
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank())
```

### Plot distribution of epithelial
```{r}
subtypes_epithelial <- data %>% 
  filter(CellType_Category == "Epithelial") %>%  
  group_by(group,
           Subclass_Cell_Identity) %>% 
  summarise(cell_count = n()) %>% 
  pivot_wider(names_from = Subclass_Cell_Identity,
              values_from = cell_count) %>%
  ungroup() %>% 
  pivot_longer("ATII_High-Surfactants":"Mystery_Disease_Epithelial",
               names_to = "CellType_Subclass",
               values_to = "cell_count") %>% 
  mutate(cell_count = replace_na(cell_count,
                                 0))

total_subtypes_per_group <- subtypes_epithelial %>% 
  group_by(group) %>% 
  summarise(total = sum(cell_count))

subtypes_epithelial <- left_join(subtypes_epithelial,
                                 total_subtypes_per_group)

plot_epithelial_subtypes <- subtypes_epithelial %>% 
    ggplot(aes(x = group,
               y = cell_count / total)) +
    geom_col(aes(fill = fct_reorder2(as.factor(CellType_Subclass),
                                     cell_count,
                                     group))) +
    scale_fill_viridis_d("Type", 
                         alpha = 0.6) +
    theme_minimal() +
    coord_flip() +
    labs(title = "Distribution of Subtypes of Ciliated Cells") +
    theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())
```

remove to save memory
```{r}
rm(cells_per_group_and_patient,
   sum_per_group,
   cells_type_per_group,
   subtypes_epithelial,
   total_subtypes_per_group)
```

# MODELS ---------------------------------------------------------------------
Take only ciliated cells,
Pivot longer to get genes as rows for modelling
Normalise the Counts

Make datasets for models 

```{r}
ciliated <- data %>% 
  filter(Subclass_Cell_Identity == "Ciliated")

ciliated_COPD <- ciliated %>% 
  pivot_longer(cols = TSPAN6:ncol(ciliated),
                 names_to = "Gene",
                 values_to = "Counts") %>% 
  select(group,
         Gene,
         Counts) %>%
  filter(!is.na(Counts)) %>% 
  filter(group != "IPF") %>%
  mutate(Counts = normalise(Counts),
         group = case_when(group == "Control" ~ 0,
                           group == "COPD" ~ 1)) %>%
  group_by(Gene) %>% 
  nest() %>% 
  ungroup()
```

Same for IPF
```{r}
ciliated_IPF <- ciliated %>% 
  pivot_longer(cols = TSPAN6:ncol(ciliated),
               names_to = "Gene",
               values_to = "Counts") %>% 
  select(group,
         Gene,
         Counts) %>%
  filter(!is.na(Counts)) %>% 
  filter(group != "COPD") %>%
  mutate(Counts = normalise(Counts),
         group = case_when(group == "Control" ~ 0,
                           group == "IPF" ~ 1)) %>% 
  group_by(Gene) %>% 
  nest() %>% 
  ungroup()


COPD_against_IPF <- ciliated %>% 
  pivot_longer(cols = TSPAN6:ncol(ciliated),
               names_to = "Gene",
               values_to = "Counts") %>% 
  select(group,
         Gene,
         Counts) %>%
  filter(!is.na(Counts)) %>% 
  filter(group != "Control") %>%
  mutate(Counts = normalise(Counts),
         group = case_when(group == "COPD" ~ 0,
                           group == "IPF" ~ 1)) %>% 
  group_by(Gene) %>% 
  nest() %>% 
  ungroup()

multinomial <- ciliated %>% 
  pivot_longer(cols = TSPAN6:ncol(ciliated),
               names_to = "Gene",
               values_to = "Counts") %>% 
  select(group,
         Gene,
         Counts) %>%
  filter(!is.na(Counts)) %>% 
  mutate(Counts = normalise(Counts),
         group = case_when(group == "Control" ~ 0,
                           group == "COPD" ~ 1,
                           group == "IPF" ~ 2)) %>% 
  group_by(Gene) %>% 
  nest() %>% 
  ungroup()
```

Now delete data as its not needed (memory) AND there is a column in the 
nested dataframes called data
```{r}
rm(data)
```

Making four models
```{r}
COPD_model <- ciliated_COPD %>% 
  mutate(mdl = map(data,
                   ~glm(group ~ Counts,
                        data = .x,
                        family = binomial(link = "logit"))),
         tidy = map(mdl,
                    tidy,
                    conf.int = TRUE)) %>% 
  unnest(tidy) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(adj_p.value = p.adjust(p.value,
                                method = "bonferroni")) %>% 
  mutate(identified_as = case_when(p.value < 0.05 ~ "significant",
                                   TRUE ~ "unsignificant")) %>% 
  mutate(adj_identified_as = case_when(adj_p.value < 0.05 ~ "significant",
                                       TRUE ~ "unsignificant"))


# wihout correction 2,042 "significant"
# with bonferroni/holm/fdr correction: 0


IPF_model <- ciliated_IPF %>% 
  mutate(mdl = map(data,
                   ~glm(group ~ Counts,
                        data = .x,
                        family = binomial(link = "logit"))),
         tidy = map(mdl,
                    tidy,
                    conf.int = TRUE)) %>%
  unnest(tidy) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(adj_p.value = p.adjust(p.value,
                                method = "bonferroni")) %>% 
  mutate(identified_as = case_when(p.value < 0.05 ~ "significant",
                                   TRUE ~ "unsignificant")) %>% 
  mutate(adj_identified_as = case_when(adj_p.value < 0.05 ~ "significant",
                                       TRUE ~ "unsignificant"))

COPD_against_IPF_model <- COPD_against_IPF %>% 
  mutate(mdl = map(data,
                   ~glm(group ~ Counts,
                        data = .x,
                        family = binomial(link = "logit"))),
         tidy = map(mdl,
                    tidy,
                    conf.int = TRUE)) %>%
  unnest(tidy) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(adj_p.value = p.adjust(p.value,
                                method = "bonferroni")) %>% 
  mutate(identified_as = case_when(p.value < 0.05 ~ "significant",
                                   TRUE ~ "unsignificant")) %>% 
  mutate(adj_identified_as = case_when(adj_p.value < 0.05 ~ "significant",
                                       TRUE ~ "unsignificant"))


multinomial_model <- multinomial %>% 
  mutate(mdl = map(data,
                   ~multinom(group ~ Counts,
                             data = .x)),
         tidy = map(mdl,
                    tidy,
                    conf.int = TRUE)) %>%
  unnest(tidy) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(adj_p.value = p.adjust(p.value,
                                method="bonferroni")) %>% 
  mutate(identified_as = case_when(p.value < 0.05 ~ "significant",
                                   TRUE ~ "unsignificant")) %>% 
  mutate(adj_identified_as = case_when(adj_p.value < 0.05 ~ "significant",
                                       TRUE ~ "unsignificant"))


sig_without_correction <- c(
  COPD_model %>% 
      filter(identified_as == "significant") %>% 
      nrow(),
  IPF_model %>% 
      filter(identified_as == "significant") %>% 
      nrow(),
  COPD_against_IPF_model %>% 
      filter(identified_as == "significant") %>% 
      nrow(),
  multinomial_model %>% 
      filter(identified_as == "significant") %>% 
      nrow()
  )


sig_with_bonferroni_correction <- c(
  COPD_model %>% 
      filter(adj_identified_as == "significant") %>% 
      nrow(),
  IPF_model %>% 
      filter(adj_identified_as == "significant") %>% 
      nrow(),
  COPD_against_IPF_model %>% 
      filter(adj_identified_as == "significant") %>% 
      nrow(),
  multinomial_model %>% 
      filter(adj_identified_as == "significant") %>% 
      nrow()
  )

significant_genes_without_correction <- tibble(sig_without_correction) %>% 
  ggplot() +
  geom_bar(aes(c("COPD_model",
                 "IPF_model",
                 "COPD_against_IPF",
                 "multinomial_model"),
               sig_without_correction),
           stat="identity",
           fill = "lightblue") +
  xlab("Model") +
  ylab("Count") +
  theme_gray() +
  ggtitle("With Boniferri Correction") +
  scale_y_continuous(expand = c(0,
                                0),
                     limits = c(0,
                                2200))
  

significant_genes_with_correction <- tibble(sig_with_bonferroni_correction) %>% 
  ggplot() +
  geom_bar(aes(c("COPD_model",
                 "IPF_model",
                 "COPD_against_IPF",
                 "multinomial_model"),
               sig_with_bonferroni_correction),
           stat ="identity",
           fill = "red") +
  xlab("Model") +
  ylab("Count") +
  theme_gray() +
  ggtitle("Without Boniferri Correction") +
  scale_y_continuous(expand = c(0,
                                0),
                     limits = c(0,
                                2200))
  

the_power_of_boniferri <- significant_genes_without_correction + 
  theme(axis.text.x = element_text(angle = 45,
                                   hjust = 1)) +
  significant_genes_with_correction +
  theme(axis.text.x = element_text(angle = 45,
                                   hjust = 1)) +
  plot_annotation(title = "Genes Identified As Significant")
```

### Save Plots
```{r}
ggsave("results/04_i_cells_count_plot.png", 
        plot = cells_count_plot)
ggsave("results/04_i_cells_type_plot.png",
        plot = cells_type_plot)
ggsave("results/04_i_plot_epithelial_subtypes.png",
        plot = plot_epithelial_subtypes)
ggsave("results/04_i_the_power_of_boniferri.png",
        plot = the_power_of_boniferri)

```

```{r ,echo=FALSE}
rm(ciliated,
   cells_count_plot,
   cells_type_plot,
   plot_epithelial_subtypes,
   the_power_of_boniferri,
   ciliated_COPD,
   ciliated_IPF,
   COPD_against_IPF,
   multinomial,
   COPD_model,
   IPF_model,
   COPD_against_IPF_model,
   multinomial_model,
   sig_without_correction,
   sig_with_bonferroni_correction,
   significant_genes_without_correction,
   significant_genes_with_correction)
```


```{r}



# Define functions -------------------------------------------------------------


source(file = "R/99_project_functions.R")

# Load data --------------------------------------------------------------------

data <- read_csv("Data/03_data.csv")

gois <- read_csv("Data/_raw/disease_genes.csv")

# Wrangling --------------------------------------------------------------------
ciliated <- data %>% 
  filter(Subclass_Cell_Identity == "Ciliated") %>% 
  pivot_longer(cols = TSPAN6:ncol(data),
               names_to = "gene",
               values_to = "Counts")

ciliated_gois <- right_join(ciliated,
                            gois, 
                            by = c("gene")) 


# Preparing the COPD model vs Control
ciliated_COPD <- ciliated_gois %>% 
  select(group,
         gene,
         Counts) %>%
  filter(!is.na(Counts)) %>% 
  filter(group != "IPF") %>% 
  mutate(Counts = normalise(Counts),
         group = case_when(group == "Control" ~ 0,
                           group == "COPD" ~ 1)) %>% 
  group_by(gene) %>% 
  nest() %>% 
  ungroup()


# Modelling the COPD vs Control 

COPD_model <- ciliated_COPD %>% 
  mutate(mdl = map(data,
                   ~glm(group ~ Counts,
                        data = .x,
                        family = binomial(link = "logit"))),
         tidy = map(mdl,
                    tidy,
                    conf.int = TRUE)) %>% 
  unnest(tidy) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(adj_p.value = 
           p.adjust(p.value,
                    method = "bonferroni")) %>% 
  mutate(identified_as = case_when(p.value < 0.05 ~ "significant",
                                   TRUE ~ "unsignificant")) %>% 
  mutate(adj_identified_as = case_when(adj_p.value < 0.05 ~ "significant",
                                       TRUE ~ "unsignificant"))


# Preparing the IPF model vs Control -------------------------------------------
ciliated_IPF <- ciliated_gois %>% 
  select(group,
         gene,
         Counts) %>%
  filter(!is.na(Counts)) %>% 
  mutate(Counts = normalise(Counts),
         group = case_when(group == "Control" ~ 0,
                           group == "IPF" ~ 1)) %>% 
  filter(group != "COPD") %>% 
  group_by(gene) %>% 
  nest() %>% 
  ungroup()

# Modelling the IPF vs Control 

IPF_model <- ciliated_IPF %>% 
  mutate(mdl = map(data,
                   ~glm(group ~ Counts,
                        data = .x,
                        family = binomial(link = "logit"))),
         tidy = map(mdl,
                    tidy,
                    conf.int = TRUE)) %>% 
  unnest(tidy) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(adj_p.value = p.adjust(p.value,
                    method="bonferroni")) %>% 
  mutate(identified_as = case_when(p.value < 0.05 ~ "significant",
                                   TRUE ~ "unsignificant")) %>% 
  mutate(adj_identified_as = case_when(adj_p.value < 0.05 ~ "significant",
                                       TRUE ~ "unsignificant"))


# Clustering, dendrograms and heatmaps -----------------------------------------

# k-means

# make data ready for kmeans function
# take only group and genes

cluster_data <- ciliated_gois %>% 
  pivot_wider(names_from = gene,
              values_from = Counts) %>% 
  select(c("group",
           CP:AFP)) %>% 
  mutate(across(everything(),
                ~replace_na(.x,
                            0)))

# discard group

data_to_cluster <- cluster_data %>% 
  select(-group)

# make 20 kmeans models with glanced column

kclusts <- 
  tibble(k = 1:20) %>%
  mutate(
    kclust = map(k, ~kmeans(data_to_cluster, .x)),
    glanced = map(kclust, glance),
  )


# unnest glanced column and plot

glanced <- kclusts %>%
  unnest(cols = glanced)

plot_k_clusters <- ggplot(glanced, 
                          aes(k,
                              tot.withinss)) +
  geom_line() +
  geom_point() +
  scale_x_discrete(breaks = 1:20,
                   limits = 1:20) +
  ylab("Total Within Sum of Squares") +
  xlab("Number of Clusters") +
  ggplot(glanced,
         aes(k,
             (betweenss / totss))) +
  geom_line() +
  geom_point() +
  scale_x_discrete(breaks = 1:20,
                    limits = 1:20) +
  ylim(0:1) +
  ylab("Between/Tot Sum of Sq.") +
  xlab("Number of Clusters") +
  plot_annotation(title = 
                    "Evaluation of K-means with Different Number of Centroids")


# shiny ------------------------------------------------------------------------

significant_identification <- function(dataset,p){
  dataset <- dataset %>% 
    mutate(identified_as = case_when(p.value < p ~ "significant",
                                     TRUE ~ "unsignificant"))
}

manhatten_plot <- function(dataset,p){
  dataset %>% 
    mutate(gene = fct_reorder(as.factor(gene),
                              p.value,
                              .desc = TRUE)) %>% 
    ggplot(aes(x = gene,
               y = p.value,
               colour = identified_as)) + 
    geom_point(size = 2) + 
    geom_hline(yintercept = p,
               linetype = "dashed") + 
    labs(x = "Gene",
         y = "p-value") +
    theme(legend.position = "bottom",
          axis.text.x = element_text(angle = 45,
                                     size = 3))
}

filter_sig_genes <- function(dataset){
  data <- dataset %>% 
    filter(identified_as == "significant") %>% 
    select("gene",
           "p.value") %>% 
    mutate(p.value = as.character(p.value))
  return(data)
}

count_sig_genes <- function(dataset){
  data <- dataset %>% 
    filter(identified_as == "significant") %>% 
    select("gene",
           "p.value") %>% 
    mutate(p.value = as.character(p.value))%>% 
    count()
  return(data)
}

ui <- fluidPage(
  theme = shinytheme("cyborg"),
  fluidRow(selectInput("selected_data",
                       label = "Select Model",
                       choice = c("IPF_vs_Control",
                                  "COPD_vs_Control"),
                       )),
  fluidRow(column(8,
                  sliderInput("p",
                              "p-value",
                              1e-3,
                              0.1,
                              value = 0.05,
                              step = 0.01)),
    column(4,
           checkboxInput("bon",
                         "Bonferroni Correction"))),
  fluidRow(plotOutput("plot"),
           h3("Genes Identified as Significant:"),
           tableOutput("sig_genes"),
           ))

server <- function(input,output,session){
  
  data <- eventReactive(
    {input$p
      input$bon
      input$selected_data},{
        if (input$selected_data == "COPD_vs_Control"){
          if (input$bon == FALSE){
            significant_identification(COPD_model,
                                       input$p)
          }
          else{
            significant_identification(COPD_model,
                                       input$p/63)
          }
        }
        else{
          if (input$bon == FALSE){
            significant_identification(IPF_model,
                                       input$p)
          }
          else{
            significant_identification(IPF_model,
                                       input$p/63)
          }
        }
      }
  )
  
  output$plot <- renderPlot(manhatten_plot(data(),
                                           input$p))
  
  output$sig_genes <- renderTable(filter_sig_genes(data()))
  
  output$gene_count <- renderTable(count_sig_genes(data()))
}

shinyApp(ui, server)



# end shiny --------------------------------------------------------------------

# dendrograms

#COPD model

dendroCOPD <- COPD_model %>%
  select(gene,p.value) %>%
  filter(!is.na(p.value))

dendroCOPD <- column_to_rownames(dendroCOPD,
                                 var = "gene")


hcCOPD<- hclust(dist(dendroCOPD), 
                "ave")

COPD_clusters <- dendro_data(hcCOPD,
                             type = "rectangle")

COPD_clusterplot <- ggplot() +
  geom_segment(data = segment(COPD_clusters), 
               aes(x = x,
                   y = y, 
                   xend = xend, 
                   yend = yend)) +
  geom_text(data = label(COPD_clusters), 
            aes(x = x,
                y = y, 
                label = label,
                hjust = 0), 
            size = 3) +
  xlab("genes") +
  ylab("distance") +
  labs(title = "COPD vs Control clustering by significance") +
  coord_flip() +
  scale_y_reverse(expand = c(0.2, 
                             0)) + 
  theme_minimal()


#IPF model

dendroIPF <- IPF_model %>%
  select(gene,p.value) %>%
  filter(!is.na(p.value))

dendroIPF <- column_to_rownames(dendroIPF,
                                var = "gene")

hcIPF<- hclust(dist(dendroIPF),
               "ave")

IPF_clusters <- dendro_data(hcIPF, 
                            type = "rectangle")

IPF_clusterplot <- ggplot() +
  geom_segment(data = segment(IPF_clusters), 
               aes(x = x, 
                   y = y, 
                   xend = xend, 
                   yend = yend)) +
  geom_text(data = label(IPF_clusters), 
            aes(x = x,
                y = y, 
                label = label,
                hjust = 0), 
            size = 3) +
  xlab("genes") +
  ylab("distance") +
  labs(title = "IPF vs Control clustering by significance") +
  coord_flip() +
  scale_y_reverse(expand = c(0.2, 
                             0)) + 
  theme_minimal()

### Start of none-tidyvers part ###

#A non-tidyverse compatible clustering method - NOT USED

d <- dist(sqrt(dendroIPF))     # change to dendroCOPD if wanted  
 
dend_row <- d %>% 
  hclust(method = "average") %>% 
  as.dendrogram 

dendro_plot <- dend_row %>% 
  highlight_branches %>% 
  plot


dend_k <- find_k(dend_row)

n_clusters <- plot(dend_k)


# Plotting the relationship of genes by significance ---------------------------

heat_dendro <- heatmaply(sqrt(dendro),
                         Colv = NULL, 
                         hclust_method = "average", 
                         fontsize_row = 8,fontsize_col = 6,
                         k_row = NA, margins = c(60,
                                                 50,
                                                 70,
                                                 90),
                         xlab = "p value", 
                         ylab = "genes", 
                         main = "Clustering of gene expression by significance",
                         plot_method = "plotly", 
                         row_dend_left = TRUE)


# PCA data preparation --- Dead end(sofus)
prefit <- ciliated_gois %>% 
  pivot_wider(names_from = "gene",
              values_from = "Counts")


prefit_genes <- prefit %>% 
  select(-mito_sum,
         -nGene,
         -nUMI) %>% 
  select(where(is.numeric))


prefit_genes %>% 
  var(na.rm = FALSE)


prefit_genes %>% 
  select(CP) %>% 
  var(na.rm = TRUE)


pca_fit <- prefit %>% 
  select(CP:ncol(prefit)) %>% # retain only numeric columns
  scale() %>% # scale data
  prcomp() # do PCA


# Save plots -------------------------------------------------------------------

ggsave("results/04_ii_plot_k_clusters.png",
       plot = plot_k_clusters)
ggsave("results/04_ii_COPD_clusterplot.png",
       plot = COPD_clusterplot)
ggsave("results/04_ii_IPF_clusterplot.png",
       plot = IPF_clusterplot)
ggsave("results/04_ii_dendro_plot.png",
       plot = dendro_plot)
ggsave("results/04_ii_heat/dendro_plot.png",
       plot = heat_dendro)

# Removing data ----------------------------------------------------------------

rm(data,
   gois,
   ciliated,
   ciliated_gois,
   ciliated_COPD,
   COPD_model,
   ciliated_IPF,
   IPF_model,
   dendro,
   d,
   dend_row,
   dendro_plot,
   dend_k,
   n_clusters,
   heat_dendro,
   prefit,
   prefit_genes,
   pca_fit,
   kclusts,
   plot_k_clusters,
   ui)
```