03-IntraIndiv-Variation.Rmd

---
title: "03-IntraIndiv Variation"
author: "Tu Hu"
date: "01/07/2022"
output: html_document
---

## Investigate genes showing high intra-individual variation

```{r}
library(dplyr)
library(SummarizedExperiment)
library(tidySummarizedExperiment)
library(tidybulk)
library(purrr)
library(BiocParallel)
library(ggplot2)
register(MulticoreParam(20))
```

### Load data
```{r}
se <- 
  readRDS("data/se.rds") %>% 
  keep_abundant(minimum_counts = 5, 
                minimum_proportion = 0)

```

### Intra Indiv Variation Index

We further analyzed the intra-individual gene expression variation. 
Within our 339 skin samples, 290 samples were paired (48 HC pairs, 52 NL pairs, and 45 LS pairs).

```{r paired samples}
paired_sample <- 
  colData(se) %>% as_tibble() %>% 
  mutate(subject_visit_skintype = paste(subject, visit, skin_type)) %>% 
  group_by(subject_visit_skintype, subject, visit, skin_type) %>% 
  summarise(n = n()) %>% 
  filter(n > 1)

n_paired_sample <- nrow(paired_sample)

paired_sample_st <-  # paired sample skin type
  paired_sample %>% 
  group_by(skin_type) %>% 
  summarise(n = n())
```

### Variation
```{r}
IntraIndivIndexCalc <- function(index){
  IVI <- 
    se[index,] %>% as_tibble() %>% 
    inner_join(paired_sample) %>% 
    # mutate(expr = log2(counts_scaled + 1)) %>% 
    group_by(skin_type) %>% 
    nest %>% 
    mutate(mean_expr = map_dbl(data, ~ .x$counts_scaled %>% mean),
           mean_var = map_dbl(data, ~ .x %>% group_by(subject, visit) %>% 
                                summarise(diff = max(counts_scaled) - min(counts_scaled)) %>% 
                                pull(diff) %>% 
                                mean),
           IVI = mean_var / mean_expr,
           gene_name = map_chr(data, ~ .x$.feature[1]))
  IVI_res <- IVI %>% select(gene_name, skin_type, mean_expr, mean_var, IVI)
  return(IVI_res)
}


```


```{r }
IntraIndivIndex <-
  runonce::save_run(
    suppressMessages(bplapply(1:nrow(se), IntraIndivIndexCalc)),
  "data/IntraIndivIndex.rds"
  )

IntraIndivIndex_t <- 
  IntraIndivIndex %>% purrr::reduce(bind_rows)

IntraIndivIndex_t <- 
  do.call(rbind.data.frame, IntraIndivIndex)

III_group_g <- 
  IntraIndivIndex_t %>% 
  ggplot(aes(IVI, color = skin_type)) +
  stat_ecdf(geom = "step", pad = FALSE) +
  coord_flip() +
  ylab("Density") +
  scale_color_manual(values = c("LS" = "#eb2d0c", 
                                "NL" = "#eb8b9b", 
                                "HC" = "#91cf60"), 
                     name = "Tissue type") + 
  theme_classic()

ggsave("figure/figure_s4a.png", 
       III_group_g, 
       width = 7 / 1.2, height = 6 / 1.5, dpi = 600)
```

```{r}
IVI <- 
  IntraIndivIndex_t %>% 
  group_by(gene_name) %>% 
  summarise(IVI_g = mean(IVI)) %>% 
  arrange(-IVI_g)

g_high_indiv_var <- 
  c("KRT25", "KRT27", "KRT28", "KRT71", "KRT85",
    "KRTAP5-7", "KRTAP5-8", "KRTAP5-10", "KRTAP9-7", "KRTAP10-4", 
    "SCYGR2", "SCYGR4")

g_high_indiv_var_g <- 
  se[g_high_indiv_var,] %>% 
  as_tibble() %>% 
  inner_join(paired_sample) %>% 
  mutate(expr = log2(counts_scaled + 1)) %>% 
  select(.feature, skin_type, subject_visit_skintype, replicate_ID, expr) %>% 
  ggplot(aes(replicate_ID, expr)) +
  geom_line(aes(group = subject_visit_skintype), color = "gray") +
  geom_violin() +
  geom_point() +
  facet_wrap(~ .feature) +
  theme_classic() +
  xlab("Biological replicate") +
  ylab("Gene expression")

ggsave("figure/figure_s4b.png", 
       g_high_indiv_var_g, 
       width = 7 , height = 6 , dpi = 600)
```