Set up environment

library(progeny)
library(dplyr)
library(Seurat)
library(ggplot2)
library(tidyr)
library(readr)
library(pheatmap)
library(tibble)
library(openxlsx)
library(sceasy)
library(reticulate)

params = list(
  seu      = "../../data/Preprocessed_data/seu_Time_series_rna.rds",
  out_name = "Time_series_rna"
)

Load data

seu <- readRDS(params$seu)
Idents(seu) <- "cell_type"
seu

An object of class Seurat 
33694 features across 7391 samples within 1 assay 
Active assay: RNA (33694 features, 3000 variable features)
 3 layers present: counts, data, scale.data
 3 dimensional reductions calculated: pca, harmony, umap

Infer pathway activity across all timepoints

seu_t1 <- subset(seu, subset= batch == unique(seu@meta.data$batch)[1])

seu_t2 <- subset(seu, subset= batch == unique(seu@meta.data$batch)[2])

seu_t3 <- subset(seu, subset= batch == unique(seu@meta.data$batch)[3])

seu_t4 <- subset(seu, subset= batch == unique(seu@meta.data$batch)[4])

UMAP

print("scRNA for all time points:")
DimPlot(seu, label = TRUE, pt.size = 0.5) + NoLegend()

print("scRNA for t1 only:")
DimPlot(seu_t1, label = TRUE, pt.size = 0.5) + NoLegend()

print("scRNA for t2 only:")
DimPlot(seu_t2, label = TRUE, pt.size = 0.5) + NoLegend()

print("scRNA for t3 only:")
DimPlot(seu_t3, label = TRUE, pt.size = 0.5) + NoLegend()

print("scRNA for t4 only:")
DimPlot(seu_t4, label = TRUE, pt.size = 0.5) + NoLegend()

[1] "scRNA for all time points:"
[1] "scRNA for t1 only:"

[1] "scRNA for t2 only:"

[1] "scRNA for t3 only:"

[1] "scRNA for t4 only:"

Heatmap

pathway_analysis <- function(seu){
    # Umap for the cell type label
    DimPlot(seu, label = TRUE, pt.size = 0.5) 

    # Assign the cluster name in dataframe
    CellsClusters <- data.frame(Cell = names(Idents(seu)), 
        CellType = as.character(Idents(seu)),
        stringsAsFactors = FALSE)

    ## Finally, we compute PROGENy pathway activity scores on the scRNA-seq data, and we then characterice the different cell populations based on these scores.
    ## We compute the Progeny activity scores and add them to our Seurat object as a new assay called Progeny. 
    seu.progeny <- progeny(seu, scale=FALSE, organism="Human", top=500, perm=1, 
        return_assay = TRUE)

    ## We can now directly apply Seurat functions in our Progeny scores. 
    ## For instance, we scale the pathway activity scores. 
    seu.progeny <- Seurat::ScaleData(seu.progeny, assay = "progeny") 

    ## We transform Progeny scores into a data frame to better handling the results
    progeny_scores_df <- 
        as.data.frame(t(GetAssayData(seu.progeny, slot = "scale.data", 
            assay = "progeny"))) %>%
        rownames_to_column("Cell") %>%
        gather(Pathway, Activity, -Cell) 

    ## We match Progeny scores with the cell clusters.
    progeny_scores_df <- inner_join(progeny_scores_df, CellsClusters)

    ## We summarize the Progeny scores by cellpopulation
    summarized_progeny_scores <- progeny_scores_df %>% 
        group_by(Pathway, CellType) %>%
        summarise(avg = mean(Activity), std = sd(Activity))

    ## We prepare the data for the plot
    summarized_progeny_scores_df <- summarized_progeny_scores %>%
        dplyr::select(-std) %>%   
        spread(Pathway, avg) %>%
        data.frame(row.names = 1, check.names = FALSE, stringsAsFactors = FALSE) 

    # We plot the different pathway activities for the different cell populations
    paletteLength = 100
    myColor = colorRampPalette(c("Darkblue", "white","red"))(paletteLength)

    progenyBreaks = c(seq(min(summarized_progeny_scores_df), 0, 
                        length.out=ceiling(paletteLength/2) + 1),
                    seq(max(summarized_progeny_scores_df)/paletteLength, 
                        max(summarized_progeny_scores_df), 
                        length.out=floor(paletteLength/2)))
                        
    progeny_hmap = pheatmap(t(summarized_progeny_scores_df[,-1]),fontsize=14, 
                            fontsize_row = 10, 
                            color=myColor, breaks = progenyBreaks, 
                            main = "PROGENy (500)", angle_col = 45,
                            treeheight_col = 0,  border_color = NA)
    return (progeny_scores_df)
}

print("pathway for all scRNA:")
all_pathwayscore <- pathway_analysis(seu)

print("pathway for t1 scRNA only:")
t1_pathwayscore <- pathway_analysis(seu_t1)
t1_pathwayscore$Condition <- 't1'

print("pathway for t2 scRNA only:")
t2_pathwayscore <- pathway_analysis(seu_t2)
t2_pathwayscore$Condition <- 't2'

print("pathway for t3 scRNA only:")
t3_pathwayscore <- pathway_analysis(seu_t3)
t3_pathwayscore$Condition <- 't3'

print("pathway for t4 scRNA only:")
t4_pathwayscore <- pathway_analysis(seu_t4)
t4_pathwayscore$Condition <- 't4'

[1] "pathway for all scRNA:"

Warning message in asMethod(object):
"sparse->dense coercion: allocating vector of size 1.9 GiB"
Warning message:
"Layer counts isn't present in the assay object; returning NULL"
Centering and scaling data matrix

Joining with `by = join_by(Cell)`
`summarise()` has grouped output by 'Pathway'. You can override using the
`.groups` argument.

[1] "pathway for t1 scRNA only:"

Warning message:
"Layer counts isn't present in the assay object; returning NULL"
Centering and scaling data matrix

Joining with `by = join_by(Cell)`
`summarise()` has grouped output by 'Pathway'. You can override using the
`.groups` argument.

[1] "pathway for t2 scRNA only:"

Warning message:
"Layer counts isn't present in the assay object; returning NULL"
Centering and scaling data matrix

Joining with `by = join_by(Cell)`
`summarise()` has grouped output by 'Pathway'. You can override using the
`.groups` argument.

[1] "pathway for t3 scRNA only:"

Warning message:
"Layer counts isn't present in the assay object; returning NULL"
Centering and scaling data matrix

Joining with `by = join_by(Cell)`
`summarise()` has grouped output by 'Pathway'. You can override using the
`.groups` argument.

[1] "pathway for t4 scRNA only:"

Warning message:
"Layer counts isn't present in the assay object; returning NULL"
Centering and scaling data matrix

Joining with `by = join_by(Cell)`
`summarise()` has grouped output by 'Pathway'. You can override using the
`.groups` argument.

timeseries_pathwayscore <- rbind(t1_pathwayscore, t2_pathwayscore, t3_pathwayscore, t4_pathwayscore)
head(timeseries_pathwayscore)

A data.frame: 6 x 5

	Cell	Pathway	Activity	CellType	Condition
	<chr>	<chr>	<dbl>	<chr>	<chr>
1	Cell_1	Androgen	0.9584746	Interstitial progenitor cell	t1
2	Cell_2	Androgen	-1.5111976	Interstitial progenitor cell	t1
3	Cell_3	Androgen	2.1025764	Interstitial progenitor cell	t1
4	Cell_4	Androgen	1.5326740	Interstitial progenitor cell	t1
5	Cell_5	Androgen	-0.4777982	Cartilage progenitor cell	t1
6	Cell_6	Androgen	1.9560105	Interstitial progenitor cell	t1