#### Load the libraries ####

library(Seurat)
library(ggplot2)
library(SingleR)
library(dplyr)
library(RColorBrewer)

# set a random seed
set.seed(1234)

#### Load filtered matrices and Create a Seurat Object for each sample ####
# 3 files alternative files are loaded in GEO repository

## E04
E04_counts <- Read10X_h5("E04/outs/filtered_feature_bc_matrix.h5")
E04_obj <- CreateSeuratObject(counts = E04_counts, project = "E04", assay = "RNA")
E04_obj # 36602 features across 226 samples within 1 assay 


## E05
E05_counts <- Read10X_h5("E05/outs/filtered_feature_bc_matrix.h5")
E05_obj <- CreateSeuratObject(counts = E05_counts, project = "E05", assay = "RNA")
E05_obj # 36602 features across 475 samples within 1 assay 

## E06
E06_counts <- Read10X_h5("E06/outs/filtered_feature_bc_matrix.h5")
E06_obj <- CreateSeuratObject(counts = E06_counts, project = "E06", assay = "RNA")
E06_obj # 36602 features across 354 samples within 1 assay

## E07
E07_counts <- Read10X_h5("E07/outs/filtered_feature_bc_matrix.h5")
E07_obj <- CreateSeuratObject(counts = E07_counts, project = "E07", assay = "RNA")
E07_obj # 36602 features across 1142 samples within 1 assay

## E08
E08_counts <- Read10X_h5("E08/outs/filtered_feature_bc_matrix.h5")
E08_obj <- CreateSeuratObject(counts = E08_counts, project = "E08", assay = "RNA")
E08_obj # 36602 features across 742 samples within 1 assay 

## E09
E09_counts <- Read10X_h5("E09/outs/filtered_feature_bc_matrix.h5")
E09_obj <- CreateSeuratObject(counts = E09_counts, project = "E09", assay = "RNA")
E09_obj # 36602 features across 1153 samples within 1 assay 

#### Merge into one Seurat Object ####

merged_obj <- merge(E04_obj, y = c(E05_obj, E06_obj, E07_obj, E08_obj, E09_obj), add.cell.ids = c("E04", "E05", "E06", "E07", "E08", "E09"), project = "MPS")
merged_obj # 36602 features across 4092 samples within 1 assay 

WPRE_raw_counts <- FetchData(object = merged_obj, vars = "WPRE", layer = "counts")
saveRDS(WPRE_raw_counts, file = "WPRE_raw_counts.rds")
merged_obj <- AddMetaData(merged_obj, WPRE_raw_counts)
saveRDS(object = merged_obj, file = "WPRE_obj.rds")