manuscript.Rmd

---
title: "Tissue-specific gene expression in obese hyperglycemic mice"
author: "Mahmoud Shaaban"
date: "12/2/2020"
output: 
    pdf_document
---

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

```{r load_libraries}
# loading required libraries
library(tidyverse)
library(DESeq2)
library(variancePartition)
library(ComplexHeatmap)
library(circlize)
library(reshape2)
library(xtable)
```

```{r raw_reads_qc}
# raw reads qc
raw_reads_qc <- list.files('data/1st RNAanalysis/1.Rawdata_information/',
                           full.names = TRUE) %>%
    map(read_tsv) %>%
    bind_rows()

df <- tibble(
    quantile = c('Q30Ratio', 'Q20Ratio'),
    quantile_name = c('> 30', '> 20')
)
(raw_reads_qc %>%
    gather(quantile, ratio, ends_with('Ratio')) %>%
    left_join(df) %>%
    ggplot(aes(x = ratio, fill = quantile_name)) +
    geom_histogram() +
    labs(x = 'Ratio of Reads',
         y = 'Number of Samples',
         fill = '') +
    theme(legend.position = c(.5, .8),
          legend.background = element_blank(),
          panel.grid = element_blank(),
          legend.key.size = unit(1,"line"))) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/raw_reads_qc.png',
           width = 5, height = 5, units = 'cm', dpi = 300)
```

```{r alignment_qc}
alignment_qc <- list.files('data/1st RNAanalysis/2.ReadAlignmentRate/',
           full.names = TRUE) %>%
    map(read_tsv) %>%
    bind_rows()

df <- tibble(
    type = c("OverallAlignmentRate", "ConcordantZero", "ConcordantPairAlignment", "MultipleAlignment" ),
    type_name = c("Overall", "Disconcordant", "Concordant", "Multiple" )
)

(alignment_qc %>%
    gather(type, value, -SampleName) %>%
    left_join(df) %>%
    mutate(value = as.numeric(str_split(value, '%', simplify = TRUE)[, 1])) %>%
    ggplot(aes(x = value, fill = type_name)) +
    geom_histogram(bins = 50) +
    labs(x = 'Ratio of Reads',
         y = 'Number of Samples',
         fill = '') +
    theme(legend.position = c(.5, .8),
          legend.background = element_blank(),
          panel.grid = element_blank(),
          legend.key.size = unit(.8,"line"))) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/alignment_qc.png',
           width = 5, height = 5, units = 'cm', dpi = 300)
```


```{r counting_qc}
df <- tibble(
    Status = c("Assigned", "Unassigned_Ambiguity", "Unassigned_MultiMapping", "Unassigned_NoFeatures", "Unassigned_Unmapped"),
    Status_name = c("Assigned", "Ambiguity", "Multi-mapping", "No Features", "Unmapped")
)
counting_qc <- read_tsv('data/1st RNAanalysis/3.FeatureCount_result/MouseRNA-seq.txt.summary')

df1 <- raw_reads_qc %>%
    dplyr::select(sample = SampleName, total = ReadCount)

(counting_qc %>%
    gather(sample, value, -Status) %>%
    left_join(df) %>%
    mutate(sample = str_split(sample, '\\.', simplify = TRUE)[, 1]) %>%
    filter(value > 0) %>%
    left_join(df1) %>%
    ggplot(aes(x = value/total, fill = Status_name)) +
    geom_histogram() +
    labs(x = 'Ratio of Reads',
         y = 'Number of Samples',
         fill = '') +
    theme(legend.position = c(.6, .7),
          legend.background = element_blank(),
          panel.grid = element_blank(),
          legend.key.size = unit(.9,"line"))) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/counting_qc.png',
           width = 5, height = 5, units = 'cm', dpi = 300)
```

```{r var_part, fig.height=4.5, fig.width=4.5}
var_part <- read_rds('data/variance_partitioned.rds')

(as_tibble(var_part) %>%
    gather() %>%
    mutate(key = str_to_title(key)) %>%
    ggplot(aes(x = key, y = value)) +
    geom_boxplot(width = .7) +
    labs(x = '', y = 'Variance Explained (%)') +
    theme(panel.grid = element_blank(),
          axis.text.x = element_text(angle = 45, hjust = 1),
          axis.title.x = element_text(size = 0))) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/var_part.png',
           height = 6, width = 6, units = 'cm', dpi = 300)
```

```{r pca, fig.height=4.5, fig.width=4.5}
ob_counts <- read_rds('data/ob_counts_batched.rds')
ob_counts <- ob_counts[, ob_counts$tissue == 'Ep']
mat <- assay(ob_counts, 'counts')

pcs <- prcomp(dist(t(mat)))

(as_tibble(scale(pcs$x[, 1:2])) %>%
    bind_cols(as_tibble(colData(ob_counts))) %>%
    ggplot(aes(x = PC1, y = PC2,
               color = paste(group, '+', diet))) +
    geom_point(size = 2) +
    labs(color = '', shape = '') +
    theme(legend.position = c(.5,.4),
          legend.background = element_blank(),
          legend.spacing.y = unit(-.1, 'in'),
          panel.grid = element_blank())) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/pca.png',
           height = 6, width = 6, units = 'cm', dpi = 300)
```

```{r distances}
ob_counts <- read_rds('data/ob_counts_batched.rds')
mat <- assay(ob_counts, 'transformed')

cc <- cor(scale(assay(ob_counts)))
png(filename = 'manuscript/figures/corr.png',
    height = 7, width = 7, units = 'cm', res = 300)

Heatmap(cc,
        row_split = ob_counts$tissue,
        column_split = ob_counts$tissue,
        column_title_side = 'bottom',
        show_heatmap_legend = FALSE,
        show_row_names = FALSE,
        show_column_names = FALSE,
        show_row_dend = FALSE,
        show_column_dend = FALSE)

dev.off()
```

```{r weights, fig.height=4.5, fig.width=4.5}
weights <- read_tsv('data/mice_measurements/weights.tsv') %>%
    mutate(group = (str_split(mouse, '_', simplify = TRUE)[, 1]),
           group = factor(group, levels = c('WT', 'ob/ob')),
           diet = str_split(mouse, '_', simplify = TRUE)[, 2],
           diet = factor(diet, levels = c('ND', 'HFD')))

(weights %>%
    ggplot(aes(x = group, y = weight, group = diet, color = diet)) +
    geom_point(position = position_dodge(width = .5), size = 2) +
    theme(legend.position = c(.99, .15),
          legend.justification = 1,
          legend.direction = 'horizontal',
          legend.background = element_blank(),
          panel.grid = element_blank()) +
    labs(x = 'Mice Group', y = 'Weight (g)', color = '') +
    lims(y = c(0, 85))) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/weights.png',
           height = 6, width = 6, units = 'cm', dpi = 300)
```

```{r blood_work, fig.height=4.5, fig.width=4.5}
blood_work <- read_tsv('data/mice_measurements/blood_work.tsv')  %>%
    gather(type, value, -mouse) %>%
    mutate(group = (str_split(mouse, '_', simplify = TRUE)[, 1]),
           group = factor(group, levels = c('WT', 'ob/ob')),
           diet = str_split(mouse, '_', simplify = TRUE)[, 2],
           diet = factor(diet, levels = c('ND', 'HFD')))  %>%
    filter(type %in% c('ALT_UL', 'AST_UL', 'Glucose_mgdL', 'T.Chol_mgdL'))

df <- tibble(
    type = c('ALT_UL', 'AST_UL', 'Glucose_mgdL', 'T.Chol_mgdL'),
    type_name = c('ALT (U/L)', 'AST (U/L)', 'Glucose (mg/dL)', 'T.Chol (mg/dL)')
)
(blood_work %>%
    left_join(df) %>%
    ggplot(aes(x = group, y = value, color = diet)) +
    geom_point(position = position_dodge(width = .5)) +
    facet_wrap(~type_name) +
    theme(legend.position = c(0, .95),
          legend.justification = 0,
          legend.direction = 'horizontal',
          legend.key.size = unit(.8,"line"),
          legend.background = element_blank(),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          panel.spacing.y = unit(0, 'mm')) +
    labs(x = 'Mice Group', y = 'Blood Measurement', color = '')) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/bloodwork.png',
           height = 8, width = 8, units = 'cm', dpi = 300)
```

```{r echo_mri, fig.height=4.5, fig.width=4.5}
(read_tsv('data/mice_measurements/echo_mri.tsv') %>%
    group_by(mouse) %>%
    summarise(`Free/Total Water` = free_water/total_water,
              `Lean/Fat` = lean/fat) %>%
    gather(type, value, -mouse) %>%
    mutate(group = (str_split(mouse, '_', simplify = TRUE)[, 1]),
           group = factor(group, levels = c('WT', 'ob/ob')),
           diet = str_split(mouse, '_', simplify = TRUE)[, 2],
           diet = factor(diet, levels = c('ND', 'HFD'))) %>%
    ggplot(aes(x = group, y = value)) +
    geom_col(width = .7) +
    facet_grid(type ~ diet, scales = 'free_y') +
    labs(x = 'Mice Group', y = 'Echo MRI (Ratio)') +
    theme(panel.grid = element_blank(),
          strip.background = element_blank(),
          panel.spacing = unit(1, 'mm'))) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/echo_mri.png',
           height = 8, width = 8, units = 'cm', dpi = 300)
```

```{r biotypes_ecdf}
ob_counts <- read_rds('data/ob_counts.rds')

fd <- rowRanges(ob_counts)
ind <- as_tibble(fd) %>%
    filter(rowMeans(assay(ob_counts)) > 0) %>%
    group_by(biotype) %>%
    summarise(n = n()) %>%
    filter(grepl('RNA', biotype), n > 100) %>%
    pull(biotype)

types <- map(ind, function(x){
    m <- assay(ob_counts)[rowRanges(ob_counts)$biotype == x,]
    as_tibble(melt(m))
}) %>%
    set_names(ind) %>%
    bind_rows(.id = 'biotype') 

(types %>%
    filter(value > 0) %>%
    mutate(biotype = ifelse(biotype == 'bidirectional_promoter_lncRNA',
                            'lncRNA', biotype)) %>%
    ggplot(aes(x = log2(value), color = biotype)) +
    stat_ecdf() +
    theme(legend.position = c(.7, .45),
          legend.background = element_blank(),
          legend.key.size = unit(.9,"line"),
          panel.grid = element_blank()) +
    labs(x = 'Counts (Log_2)', y = 'ECDF', color = '')) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/biotypes_ecdf.png',
           height = 6, width = 6, units = 'cm', dpi = 300)
```

```{r biotype_count}
(types %>%
    mutate(biotype = ifelse(biotype == 'bidirectional_promoter_lncRNA',
                            'lncRNA', biotype)) %>%
    group_by(biotype) %>%
    summarise(n = length(unique(Var1))) %>%
    mutate(biotype = ifelse(biotype == 'bidirectional_promoter_lncRNA',
                            'lncRNA', biotype),
           biotype =  fct_reorder(biotype, n)) %>%
    ggplot(aes(x = biotype, y = n)) +
    geom_col(width = .7) +
    labs(x = '', y = 'Number of Features') +
    coord_flip() +
    theme(panel.grid = element_blank(),
          axis.title.y = element_text(size = 0),
          axis.text.x = element_text(size = 8))) %>%
    ggsave(plot = .,
           filename = 'manuscript/figures/biotypes_count.png',
           height = 6, width = 6, units = 'cm', dpi = 300)
```

```{r receptors_hy}
ob_counts <- read_rds('data/ob_counts_batched.rds')
mat <- assay(ob_counts, 'removed_tissue')
ind <- c(Lepr = 'ENSMUSG00000057722',
         Npy = 'ENSMUSG00000029819',
         Ghsr = 'ENSMUSG00000051136',
         Agrp = 'ENSMUSG00000005705')

m <- mat[rownames(ob_counts) %in% unlist(ind), ob_counts$tissue == 'Hy']
rownames(m) <- names(ind)[match(rownames(m), ind)]
groups <- as_tibble(colData(ob_counts))[ob_counts$tissue == 'Hy',]
rownames(m)
png(filename = 'manuscript/figures/receptors_hy.png',
    height = 7, width = 7, units = 'cm', res = 300)
Heatmap(m,
        column_split = (groups$group),
        show_row_dend = FALSE,
        row_title = NULL,
        show_heatmap_legend = FALSE,
        column_labels = groups$diet,
        show_column_dend = FALSE)
dev.off()
```

```{r receptors_ep}
ob_counts <- read_rds('data/ob_counts_batched.rds')
mat <- assay(ob_counts, 'removed_tissue')
ind <- c(Ins1 = 'ENSMUSG00000035804',
         Irs1 = 'ENSMUSG00000055980',
         Pik3r1 = 'ENSMUSG00000041417',
         Slc2a4 = 'ENSMUSG00000018566',
         Slc2a2 = 'ENSMUSG00000027690',
         Prkaa1 = 'ENSMUSG00000050697',
         Prkaa2 = 'ENSMUSG00000028518')

tissue <- ob_counts$tissue %in% 'Ep'

m <- mat[rownames(ob_counts) %in% unlist(ind), tissue]
rownames(m) <- names(ind)[match(rownames(m), ind)]
groups <- as_tibble(colData(ob_counts))[tissue,]

png(filename = 'manuscript/figures/receptors_ep.png',
    height = 7, width = 7, units = 'cm', res = 300)
Heatmap((m),
        column_split = (groups$group),
        row_split = c(1, 2, 2, 1, 2, 2, 1),
        row_title = NULL,
        show_row_dend = FALSE,
        show_heatmap_legend = FALSE,
        column_labels = groups$diet,
        show_column_dend = FALSE)
dev.off()
```

```{r over_representation}
ov <- read_rds('data/over_representation.rds')
ov %>%
    filter(p.adjust < .1) %>%
    mutate(n1 = str_split(GeneRatio, '/', simplify = TRUE)[, 1],
           n2 = str_split(GeneRatio, '/', simplify = TRUE)[, 2],
           ratio = round(as.numeric(n1)/as.numeric(n2), 2)) %>%
    select(variable, Description, ratio) %>%
    arrange(desc(ratio)) %>%
    group_by(variable) %>%
    dplyr::slice(1:5) %>%
    ungroup() %>%
    mutate(variable = ifelse(duplicated(variable), '', str_to_title(variable))) %>%
    setNames(c('Variable', 'Term', 'Ratio')) %>%
    xtable(align = 'cllc') %>%
    print(floating = FALSE,
          include.rownames = FALSE,
          booktabs = TRUE,
          sanitize.text.function = identity,
          comment = FALSE,
          add.to.row = list(pos = list(5, 10),
          command = rep('\\midrule ', 2)),
          file = 'manuscript/tables/over-representation.tex')
```

```{r dataset}
ob_counts <- read_rds('data/ob_counts.rds')
pd <- as_tibble(colData(ob_counts))

pd %>%
    group_by(group, diet) %>%
    summarise(n_tissue = length(unique(tissue)),
              n = length(unique(id))) %>%
    ungroup() %>%
    mutate(group = ifelse(duplicated(group), '', (group))) %>%
    setNames(c('Group', 'Diet', 'Tissue (N)', 'Samples (N)')) %>%
    xtable(align = 'cllcc') %>%
    print(floating = FALSE,
          include.rownames = FALSE,
          booktabs = TRUE,
          sanitize.text.function = identity,
          comment = FALSE,
          add.to.row = list(pos = list(2),
          command = '\\midrule '),
          file = 'manuscript/tables/dataset.tex')
```