Introductions

Pseudotime trajectory analysis infers the underlying continuous biological process by computationally ordering cells based on their molecular profiles measured in single-cell sequencing studies. Spatial transcriptomics measure both the gene expression and spatial locations of cells, thus providing opportunity to study the spatial heterogeneity of continuous biolgoical processes.

Paella is a computational method that decomposes the pseudotime cell trajectory given by the user into multiple spatial sub-trajectories that are both temporal and spatial continuous. It provides utilities to identify genes with differential expression patterns across spatial sub-trajectories and visualize the gene expressino patterns.

Load packages

library(Paella)
## Warning: replacing previous import 'circlize::degree' by 'igraph::degree' when
## loading 'Paella'

Breast cancer Visium spatial transcriptomics example

In this vignette, we demonstrate Paella using a breast cancer Visium spatial transcriptomics example. We manually annotated six in situ carcinoma and defined cell trajectories by the distance of each spot to the closest center of in situ carcinoma. For demonstration purposes, we only included 100 genes randomly selected from the original full gene expression matrix.

Below we read in the three inputs required by Paella: pseudotime already identified by user, spatial coordinates of cell/spots, and gene expression matrix

# Read in pseudotime data
pt <- readRDS(paste0(system.file('data',package = 'Paella'),'/pt.rds'))
head(pt)
## AAACAGAGCGACTCCT-1 AAATAACCATACGGGA-1 AAATTACACGACTCTG-1 AACCTTTACGACGTCT-1 
##           6.442731           1.176346           2.002196           5.351055 
## AACGTGATGAAGGACA-1 AACTCCTAATCCCATG-1 
##           6.057746           4.431003
# Read in spatial coordinates
cd <- readRDS(paste0(system.file('data',package = 'Paella'),'/cd.rds'))
head(cd)
##                    [,1] [,2]
## AAACAGAGCGACTCCT-1  -94  -14
## AAATAACCATACGGGA-1  -88  -14
## AAATTACACGACTCTG-1  -86  -14
## AACCTTTACGACGTCT-1  -92  -16
## AACGTGATGAAGGACA-1  -83   -9
## AACTCCTAATCCCATG-1  -92  -12
# Read in gene expression matrix
mat <- readRDS(paste0(system.file('data',package = 'Paella'),'/expr.rds'))
str(mat)
##  num [1:100, 1:329] 0.0426 0.4425 0.4179 0.03 0.2975 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:100] "PLCD1" "NDUFV3" "DHTKD1" "MCOLN3" ...
##   ..$ : chr [1:329] "AAACAGAGCGACTCCT-1" "AAATAACCATACGGGA-1" "AAATTACACGACTCTG-1" "AACCTTTACGACGTCT-1" ...

Running Paella

We run Paella to identify spatial sub-trajectories, which returns a list. The first element of the list stores the identities of spatial sub-trajectories.

paellares <- paella(pt,cd)
table(paellares[[1]])
## 
## traj1 traj2 traj3 traj4 traj5 traj6 
##    64    55    49    66    38    57

The spatial sub-trajectories can be visualized using the below command

spatialplot(paellares)

Differential analysis

We next identify genes that have differential temporal gene expression patterns across spatial sub-trajectories. The result is a data.frame containing FDR and pvalues for three types of differential tests: overall, mean change, and trend change. The final column summarizes the type of differential.

diff <- difftest(mat,paellares)
head(diff)

Get a list of all overall differential genes:

siggene <- rownames(diff)[diff$Overall_FDR < 0.05]

Visulize the expression of differential genes

We first get a smooth fit of temporal gene expression dynamics for the differential genes. The fit is done using mgcv package.

genefit <- genefit(mat,paellares,siggene)

Next we assign genes into multiple clusters. If no cluster number is given, findPC is used to automatically determine number of clusters. The first element of geneclu stores the cluster assignments.

geneclu <- clustergene(genefit)         
head(geneclu[[1]])
##        AARD         EVL       PRDX4 PALM2-AKAP2       EDNRA       TRPC6 
##           1           2           3           4           4           4

Finally we visualize the temporal expression patterns of differential genes using a heatmap:

geneheatmap(genefit,geneclu,diff)

We can also visualize the top six genes using scatterplots

genescatter(mat,genefit,paellares)

Session Info

sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] Paella_1.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.7           lattice_0.20-44      circlize_0.4.15     
##  [4] png_0.1-7            assertthat_0.2.1     digest_0.6.27       
##  [7] foreach_1.5.1        utf8_1.2.1           plyr_1.8.6          
## [10] R6_2.5.0             magic_1.6-0          stats4_4.1.0        
## [13] evaluate_0.14        highr_0.9            ggplot2_3.3.5       
## [16] pillar_1.6.1         GlobalOptions_0.1.2  rlang_0.4.11        
## [19] jquerylib_0.1.4      S4Vectors_0.30.2     GetoptLong_1.0.5    
## [22] Matrix_1.3-3         rmarkdown_2.10       labeling_0.4.2      
## [25] splines_4.1.0        stringr_1.4.0        igraph_1.2.6        
## [28] munsell_0.5.0        compiler_4.1.0       xfun_0.24           
## [31] pkgconfig_2.0.3      BiocGenerics_0.38.0  shape_1.4.6         
## [34] mgcv_1.8-35          findPC_1.0           htmltools_0.5.1.1   
## [37] tidyselect_1.1.1     tibble_3.1.2         IRanges_2.26.0      
## [40] codetools_0.2-18     matrixStats_0.59.0   fansi_0.5.0         
## [43] crayon_1.4.1         dplyr_1.0.7          grid_4.1.0          
## [46] nlme_3.1-152         jsonlite_1.7.2       gtable_0.3.0        
## [49] lifecycle_1.0.0      DBI_1.1.1            magrittr_2.0.1      
## [52] scales_1.1.1         stringi_1.6.2        farver_2.1.0        
## [55] reshape2_1.4.4       sp_1.4-5             doParallel_1.0.16   
## [58] geometry_0.4.6.1     bslib_0.2.5.1        ellipsis_0.3.2      
## [61] vctrs_0.3.8          generics_0.1.0       rjson_0.2.21        
## [64] RColorBrewer_1.1-2   iterators_1.0.13     tools_4.1.0         
## [67] Cairo_1.6-0          glue_1.4.2           purrr_0.3.4         
## [70] abind_1.4-5          parallel_4.1.0       yaml_2.2.1          
## [73] clue_0.3-59          colorspace_2.0-3     cluster_2.1.2       
## [76] ComplexHeatmap_2.8.0 knitr_1.33           patchwork_1.1.1     
## [79] sass_0.4.0