---
title: "Generate CPM by developmental stage"
author: "Ioannis Sarropoulos"
date: '`r format(Sys.time(), "%d %B, %Y")`'
output: html_document
---

```{r}
suppressPackageStartupMessages({
  library(tidyverse)
  library(Matrix)
  library(SummarizedExperiment)
})
```


```{r}
## Read summarised experiment object
peak_se <- readRDS(url("https://apps.kaessmannlab.org/mouse_cereb_atac/Mouse_Cereb_ATAC_CREs_SE.rds"))

## Access peak x cell matrix
peak_mat <- assay(peak_se)

## Access peak info
peak_info <- rowData(peak_se)
dim(peak_info)

## Access cell info
cell_info <- colData(peak_se)
dim(cell_info)
```

We are aggregating CRE accessibility within each sample:

```{r}
table(cell_info$Sample_short)

## Iterate over all samples
pseudobulk_mat <- do.call(cbind, lapply(sort(unique(cell_info$Sample_short)), function(x) {
  ## Get cells from that sample
  cells <- cell_info$Cell[cell_info$Sample_short==x]
  ## Aggregate CRE accessibility across all cells
  return(Matrix::rowSums(peak_mat[, cells]))
}))

dim(pseudobulk_mat)
## Seting column names to samples
colnames(pseudobulk_mat) <- sort(unique(cell_info$Sample_short))
```
CPM scaling:

```{r}
pseudobulk_mat_cpm <- t(t(pseudobulk_mat)/Matrix::colSums(pseudobulk_mat)) * 1e6
```
Export:

```{r}
write.table(pseudobulk_mat_cpm, "/home/public_data/mouse_cereb_atac/Mouse_pseudobulk_bySample_cpm.txt", col.names = T, row.names = T, quote = F, sep = "\t")
```

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

