interstitial_report7_population_timeline.Rmd

---
title: "scRNAseq Interstitial cells : What happens over time?"
author: "Marion Hardy"
date: "`r Sys.Date()`"
output: 
  html_document: 
    toc: yes
    theme: spacelab
    highlight: monochrome
editor_options: 
  chunk_output_type: console
  markdown: 
    wrap: 72
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(message = FALSE, cache = T, echo = FALSE, warning = F, cache.lazy = F)
knitr::opts_chunk$set(fig.width=8, fig.height=6) 

library(tidyverse)
library(openxlsx)
library(Seurat)
library(scCustomize)
library(ggpubr)
```

# Introduction

Data from our collaborator Panagiotis rds file for a
SingleCellExperiment object containing the single cell data for **the
interstitial cells** of *Hydra Vulgaris* during multiples stages of
regeneration after bisection:

<https://www.dropbox.com/scl/fi/dg76sigjnj5u6qr06xk79/sce_interstitial_Juliano.rds?rlkey=f0y3lqt0wdcwq652zisnrpf9t&dl=0>

BUT they mapped it to *Hydra Magnipapillata (102 version of 105)* "Quantification of the generated single cell libraries was performed using the Salmon-Alevin software suite 
(Salmon version 1.6.0) against the ncbi Hydra 102 transcriptome."

The coldata of the object contain cell annotation including

-   quality metrics: nFeature nCount (not MT percentage interestingly)

-   batch info: either 2869 (3162 barcodes), 3113 (10352 barcodes), 3271
    (13279 barcodes), 3357 (3875 barcodes)

-   originating experiments (head or foot regeneration)

-   experimental time points

-   pseudo-axis assignment (vals.axis ranging from 0-1, increasing in
    the foot-tentacle direction)

-   mitotic and apoptotic signatures indices from 0 to 1

The rowdata contains gene annotation, using Entrez-gene identifiers. I
have also noticed that in the sce objects there's

-   PCA, tSNE and UMAP coordinates for reduced dimensions + corrected
    for batch values

-   assay metafeatures hold gene_id, product, gene, is.rib.prot.gene
    (T/F), HypoMarkers (T/F), ccyle (T/F), apopt (T/F) etc

I converted the sce objects 6.6gb into a seurat object 2.2 gb and
checked that all cited parameters could be found in it.

# Summary of 3/06 meeting with Hannah and Celina

-   check myb as a marker for bipotential progenitors
-   relabel some of the clusters (ec3, Gc, endodermal neurons)
-   plot cell cycle scores
-   plot pseudo-axis


# Data loading and upate on the project

This report is a repeat of report 5 but addresses the above-mentioned comments from the meeting with Celina and Hannah.

In this report we're gonna plot transcription factors of interest over the different timepoints.
I also wanted to add cell population counts over time.
Maybe I could learn to plot gene expression over time too?
Anyway, we're exploring gene expression programs and cell populations over time.

```{r}
# seurat objects
seurat_f = readRDS("./data_output/interstitial/seurat_filtered.rds")
# transcription factors of interest
tr = read.xlsx("./data/mcbi_dataset_MH_annotated.xlsx", sheet = "transcription")
# reorder the timepoints

seurat_f$Timepoint <- factor(seurat_f$Timepoint, 
                            levels=c("REG_HEAD_t0", 
                                     "REG_HEAD_t06", 
                                     "REG_HEAD_t12", 
                                     "REG_HEAD_t24", 
                                     "REG_HEAD_t48", 
                                     "REG_HEAD_t96",
                                     "REG_FOOT_t0", 
                                     "REG_FOOT_t06", 
                                     "REG_FOOT_t12", 
                                     "REG_FOOT_t24", 
                                     "REG_FOOT_t48", 
                                     "REG_FOOT_t96"))

seurat_f$vals.axis_k25 <- as.numeric(seurat_f$vals.axis_k25)
seurat_f$vals.ccycle <- as.numeric(seurat_f$vals.ccycle)
seurat_f$vals.apopt <- as.numeric(seurat_f$vals.apopt)

```

## Changing UMAP resolution

```{r}
Idents(seurat_f) = "seurat_clusters"

DimPlot(seurat_f, label = TRUE)+
  labs(title = "Resolution = 0.1")

```

```{r, echo=TRUE}

target <- c("NSP4","midasin","mini-COL8","SP-D-like","FH20-X3","CAII",
            "CELA3B","zinc-carboxypeptidase","ANO39","TYN1","H2A.2.2","nas2-X2",
            "Lwamide-X1","DMRT1","TUBA4A-X1","HTRA3","polRF-X1","ec3A","ec3B",
            "nop58-X1","OTP","MEP1A","rhammosyl-O-methyltransferase","hywnt3",
            "PPOD1","ks1","hyAlx","ELAV2","POU4","MUC2", "ec3", "grm1","myc",
            "myc1","wnt3","ec2","ec1B","ec4","ec1A","en1","en1_NDF1_DANRE",
            "ec5","en1","en2A","en2B","myb") 

```

 I founda myb by ctrl+f in the annotation file.
 However, in Stefan's supplementary analysis, myb is described as "g27424"
 after blasting, I only got a 32% hit (80% cover) in the 102 annotation file...
 
 Uncharacterized protein LOC105848384 [Hydra vulgaris]	Hydra vulgaris	243	243	82%	2e-68	32.09%	1023	XP_047129678.1

```{r}
Idents(seurat_f) = "seurat_clusters"

p1 = Clustered_DotPlot(seurat_f, features = target, k = 6)
```

```{r}
seurat_f@meta.data <-
  seurat_f@meta.data %>%
  mutate(attr_clusters=
           ifelse(seurat_clusters == "0", "ISC",
           ifelse(seurat_clusters == "1", "Nb", 
           ifelse(seurat_clusters == "2", "GranG/ZymoG",
           ifelse(seurat_clusters == "3", "earlyGc", 
           ifelse(seurat_clusters == "4", "Nb",
           ifelse(seurat_clusters == "5", "earlyNem",
           ifelse(seurat_clusters == "6", "Nb", 
           ifelse(seurat_clusters == "7", "Nb",
           ifelse(seurat_clusters == "8", "earlyNeur", 
           ifelse(seurat_clusters == "9", "Gc",
           ifelse(seurat_clusters == "10", "Nb", 
           ifelse(seurat_clusters == "11", "ec1A",
           ifelse(seurat_clusters == "12", "ec3", 
           ifelse(seurat_clusters == "13", "Nb",
           ifelse(seurat_clusters == "14", "ec2",
           ifelse(seurat_clusters == "15", "Nb", 
           ifelse(seurat_clusters == "16", "ec1A/ec1B",
           ifelse(seurat_clusters == "17", "ec3", 
           ifelse(seurat_clusters == "18", "en",
           ifelse(seurat_clusters == "19", "ec4", 
           ifelse(seurat_clusters == "20", "Nb",
           NA))))))))))))))))))))))
```

```{r}
Idents(seurat_f) = "attr_clusters"

DimPlot(seurat_f, label = TRUE)+
  labs(title = "Labelled UMAP")

```


## Subsetting by head/foot and timepoint

```{r}

Idents(seurat_f) = "Timepoint"
DimPlot(seurat_f, label = FALSE)

# Creating the head and foot data subset
hfootf = subset(seurat_f, subset = EXP == "REG_FOOT")
hheadf = subset(seurat_f, subset = EXP == "REG_HEAD")

```

```{r, fig.width=30, fig.height=6}

# Plotting conditions separately
Idents(hfootf) = "attr_clusters"
DimPlot(hfootf, reduction = "umap", split.by = "Timepoint", pt.size = 0.5)

Idents(hheadf) = "attr_clusters"
DimPlot(hheadf, reduction = "umap", split.by = "Timepoint", pt.size = 0.5)

```

### Pseudo axis values
```{r, fig.height=7, fig.width=35} 

#annot run this for some reason my computer gets a conniption

FeaturePlot(object = hheadf, 
            features = "vals.axis_k25",
            reduction = "umap",
            order = TRUE,
            min.cutoff = 'q10', 
            label = F,
            repel = TRUE,
            pt.size = 1.2, 
            split.by = "Timepoint")

```

```{r, fig.height=7, fig.width=35} 

#annot run this for some reason my computer gets a conniption

FeaturePlot(object = hfootf, 
            features = "vals.axis_k25",
            reduction = "umap",
            order = TRUE,
            min.cutoff = 'q10', 
            label = F,
            repel = TRUE,
            pt.size = 1.2, 
            split.by = "Timepoint")

```

### Cell cycle values
```{r, fig.height=7, fig.width=35} 

#annot run this for some reason my computer gets a conniption

FeaturePlot(object = hheadf, 
            features = "vals.ccycle",
            reduction = "umap",
            order = TRUE,
            min.cutoff = 'q10', 
            label = F,
            repel = TRUE,
            pt.size = 1.2, 
            split.by = "Timepoint")

```

```{r, fig.height=7, fig.width=35} 

#annot run this for some reason my computer gets a conniption

FeaturePlot(object = hfootf, 
            features = "vals.ccycle",
            reduction = "umap",
            order = TRUE,
            min.cutoff = 'q10', 
            label = F,
            repel = TRUE,
            pt.size = 1.2, 
            split.by = "Timepoint")

```

### Apoptosis score
```{r, fig.height=7, fig.width=35} 

#annot run this for some reason my computer gets a conniption

FeaturePlot(object = hheadf, 
            features = "vals.apopt",
            reduction = "umap",
            order = TRUE,
            min.cutoff = 'q10', 
            label = F,
            repel = TRUE,
            pt.size = 1.2, 
            split.by = "Timepoint")

```

```{r, fig.height=7, fig.width=35} 

#annot run this for some reason my computer gets a conniption

FeaturePlot(object = hfootf, 
            features = "vals.apopt",
            reduction = "umap",
            order = TRUE,
            min.cutoff = 'q10', 
            label = F,
            repel = TRUE,
            pt.size = 1.2, 
            split.by = "Timepoint")

```

```{r, eval=FALSE}
write_rds(seurat_f, "./data_output/interstitial/seurat_filtered_res01.rds")
```

## Cell population evolution over time
### Regenerating a head
#### "Raw" cell numbers

```{r General statistics levels}

aa = table(hheadf$attr_clusters, hheadf$Timepoint) 
aa = rbind(aa, Sum = colSums(aa))
aa = aa[,1:6]

aa %>% 
  knitr::kable()
```

```{r, fig.height=6, fig.width=12}

aadf =
  aa %>% as.data.frame()

aadf$celltype = rownames(aadf)

aadf = 
aadf %>% 
  pivot_longer(names_to = "Timepoint",
               values_to = "N_cells",
               cols = 1:6)
p1 =   
aadf %>% 
  filter(celltype != "Sum") %>% 
  ggplot(aes(x = Timepoint,
             y = N_cells, group = celltype, color = celltype))+
  geom_line(linewidth = 1)+ 
  theme_bw()+
  theme(axis.text.x = element_text(angle = 70, vjust = 1, hjust=1))+
  labs(title = "Regenerating a head", subtitle = "Cell number per population over time")

p2 =
aadf %>% 
  filter(celltype != "Sum", celltype != "Nb") %>% 
  ggplot(aes(x = Timepoint,
             y = N_cells, group = celltype, color = celltype))+
  geom_line(linewidth = 1)+ 
  theme_bw()+
  theme(axis.text.x = element_text(angle = 70, vjust = 1, hjust=1))+
  labs(title = "Regenerating a head", subtitle = "Removed nematoblasts")

p1+p2
```

#### Relative number of cells

```{r, fig.width=12, fig.height=6}

vec = colSums(aa[1:13,])
vec = vec[vec!=0]

percent = round(sweep(aa,2,vec, '/')*100,2)

percent %>% 
  knitr::kable()

percentdf = as.data.frame(percent)
percentdf$celltype = rownames(percentdf)

percentdf = 
percentdf %>% 
  pivot_longer(names_to = "Timepoint",
               values_to = "N_cells",
               cols = 1:6)

p1 = 
percentdf %>% 
  filter(celltype != "Sum") %>% 
  ggplot(aes(x = Timepoint,
             y = N_cells, group = celltype, color = celltype))+
  geom_line(linewidth = 1)+ 
  theme_bw()+
  theme(axis.text.x = element_text(angle = 70, vjust = 1, hjust=1))+
  ylab("% cell population")+
  labs(title = "Regenerating a head", subtitle = "Relative cell % per population over time")

p2=
percentdf %>% 
  filter(!celltype %in% c("Sum", "Nb", "ec2","ec3n/en1n","ec4","doublets/triplets", 
                          "unknown/ZymoG?","ec1A","ec1A/ec1B")) %>% 
  ggplot(aes(x = Timepoint,
             y = N_cells, group = celltype, color = celltype))+
  geom_line(linewidth = 1)+ 
  theme_bw()+
  theme(axis.text.x = element_text(angle = 70, vjust = 1, hjust=1))+
  ylab("% cell population")+
  labs(title = "Regenerating a head", subtitle = "Removed multiple cell types")

ggarrange(p1,p2,ncol = 2, nrow = 1)

```

Notice how early Gc decrease almost 5x between 12-24 hours while ISC, GranG and early neurons spike up 2 to 3X. I think that population is called upon to create not only gland cells but can help support other progenitors.

### Regenerating a foot
#### "Raw" cell numbers
```{r}
ab = table(hfootf$attr_clusters, hfootf$Timepoint) 
ab = rbind(ab, Sum = colSums(ab))
ab = ab[,7:12]  

ab %>% 
  knitr::kable()
```

```{r, fig.height=6, fig.width=12}
abdf =
  ab %>% as.data.frame()

abdf$celltype = rownames(abdf)

abdf = 
abdf %>% 
  pivot_longer(names_to = "Timepoint",
               values_to = "N_cells",
               cols = 1:6)
p1 =  
abdf %>% 
  filter(celltype != "Sum") %>% 
  ggplot(aes(x = Timepoint,
             y = N_cells, group = celltype, color = celltype))+
  geom_line(linewidth = 1)+ 
  theme_bw()+
  theme(axis.text.x = element_text(angle = 70, vjust = 1, hjust=1))+
  labs(title = "Regenerating a foot", subtitle = "Cell number per population over time")

p2 =
abdf %>% 
  filter(celltype != "Sum", celltype != "Nb") %>% 
  ggplot(aes(x = Timepoint,
             y = N_cells, group = celltype, color = celltype))+
  geom_line(linewidth = 1)+ 
  theme_bw()+
  theme(axis.text.x = element_text(angle = 70, vjust = 1, hjust=1))+
  labs(title = "Regenerating a foot", subtitle = "Removed nematoblasts")

ggarrange(p1,p2,ncol = 2, nrow = 1)

```

#### Relative number of cells

```{r, fig.width=12, fig.height=6}

vec = colSums(ab[1:13,])
vec = vec[vec!=0]

percent = round(sweep(ab,2,vec, '/')*100,2)

percent %>% 
  knitr::kable()

percentdf = as.data.frame(percent)
percentdf$celltype = rownames(percentdf)

percentdf = 
percentdf %>% 
  pivot_longer(names_to = "Timepoint",
               values_to = "N_cells",
               cols = 1:6)

p1 = 
percentdf %>% 
  filter(celltype != "Sum") %>% 
  ggplot(aes(x = Timepoint,
             y = N_cells, group = celltype, color = celltype))+
  geom_line(linewidth = 1)+ 
  theme_bw()+
  ylab("% cell population")+
  theme(axis.text.x = element_text(angle = 70, vjust = 1, hjust=1))+
  labs(title = "Regenerating a foot", subtitle = "Relative cell % per population over time")

p2=
percentdf %>% 
  filter(!celltype %in% c("Sum", "Nb", "ec2","ec3n/en1n","ec4","doublets/triplets", 
                          "unknown/ZymoG?","ec1A","ec1A/ec1B")) %>% 
  ggplot(aes(x = Timepoint,
             y = N_cells, group = celltype, color = celltype))+
  geom_line(linewidth = 1)+ 
  theme_bw()+
  ylab("% cell population")+
  theme(axis.text.x = element_text(angle = 70, vjust = 1, hjust=1))+
  labs(title = "Regenerating a foot", subtitle = "Removed multiple cell types")

ggarrange(p1,p2,ncol = 2, nrow = 1)

```

Notice how early neurons only go up after a spike in earlGc. Is the opposite happening here?

## Transcription factors

We are looking at transcription factors of interest that, based on bulk RNAseq, appeared active only in the injured animal, not the homeostatic one.

### Dotplot

```{r}
Idents(hfootf) = "Timepoint"
Idents(hheadf) = "Timepoint"

DotPlot_scCustom(hfootf, features = tr$Symbol_updated, 
                      colors_use = viridis_plasma_dark_high,
                      x_lab_rotate = T)+
  coord_flip()

```

```{r}
DotPlot_scCustom(hheadf, features = tr$Symbol_updated, 
                      colors_use = viridis_plasma_dark_high,
                      x_lab_rotate = T)+
  coord_flip()
```

### Feature plot
```{r, fig.height=10, fig.width=15}

FeaturePlot(object = seurat_f, 
            features = tr$Symbol_updated,
            reduction = "umap",
            order = TRUE,
            min.cutoff = 'q10', 
            label = F,
            repel = TRUE,
            pt.size = 0.75)

```

### Regenerating a head
```{r, fig.height=45, fig.width=30}

FeaturePlot(object = hheadf, 
            features = tr$Symbol_updated,
            reduction = "umap",
            order = TRUE,
            min.cutoff = 'q10', 
            label = F,
            repel = TRUE,
            pt.size = 1.2,
            split.by = "Timepoint")

```

It looks like the appearing cells might be the one with the transcription factor activity

### Regenerating a foot
```{r, fig.height=45, fig.width=30}

FeaturePlot(object = hfootf, 
            features = tr$Symbol_updated,
            reduction = "umap",
            order = TRUE,
            min.cutoff = 'q10', 
            label = F,
            repel = TRUE,
            pt.size = 1.2,
            split.by = "Timepoint")

```

```{r, eval = FALSE}
write_rds(seurat_f, "./data_output/interstitial/seurat_filtered.rds")
write_rds(neurons, "./data_output/interstitial/neurons.rds")
```

# Session info

```{r}
sessionInfo()
```