Day 1 // Hands-on 2 // Explore SingleCellExperiment

Load packages

library( "Seurat" )
library( "scater" )
library( "scran" )

library( "tidyverse" )

Load data

Seurat -> SingleCellExperiment

sce <- readRDS( "/home/idies/workspace/practical_genomics/day1/retina.rds" )
sce <- as.SingleCellExperiment( sce )

Check data

sce

class: SingleCellExperiment 
dim: 21599 8217 
metadata(0):
assays(1): logcounts
rownames(21599): FAM138A RP11-34P13.7 ... AC145205.1
  AC145205.1-1
rowData names(1): hgnc_gene_symbol
colnames(8217): AAACCCAGTCGGCACT-1 AAACGAACAGGGACTA-1 ...
  TTTGGAGCAGTTTGGT-6 TTTGGAGTCGTTCTAT-6
colData names(18): final_cluster_labels libraryID ...
  cell_type ident
reducedDimNames(1): UMAP
altExpNames(0):

1. assay()

Introduce

• How to access the expression data?

help( assay )

m_assay <- assay( sce, "logcounts" )

• What class of object is this?

class( m_assay )

[1] "matrix" "array" 

• How large is it?

dim( m_assay )

[1] 21599  8217

• How to subset first 10 rows, first 3 columns?

m_assay[ 1:10, 1:3 ]

              AAACCCAGTCGGCACT-1 AAACGAACAGGGACTA-1 AAACGAAGTTCTCCAC-1
FAM138A                        0           0.000000                  0
RP11-34P13.7                   0           0.000000                  0
AL627309.1                     0           0.000000                  0
AP006222.2                     0           1.443678                  0
RP4-669L17.10                  0           0.000000                  0
RP5-857K21.4                   0           0.000000                  0
RP11-206L10.3                  0           0.000000                  0
RP11-206L10.2                  0           0.000000                  0
RP11-206L10.9                  0           0.000000                  0
FAM87B                         0           0.000000                  0

Subset RHO

x_rho <- m_assay[ "RHO", ]

class( x_rho )

[1] "numeric"

length( x_rho )

[1] 8217

hist( x_rho )

summary( x_rho )

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.0000  0.0000  0.7655  1.0630  5.9748 

Subset RHO, VIM

goi <- c( "RHO", "VIM" )
x_rhovim <- m_assay[ goi, ]

class( x_rhovim )

[1] "matrix" "array" 

dim( x_rhovim )

[1]    2 8217

hist( x_rhovim )

ggplot RHO, VIM

• t(), data.frame(), rownames_to_column(), pivot_longer()

df_rhovim <- t( x_rhovim )
df_rhovim <- data.frame( df_rhovim )
df_rhovim <- rownames_to_column( df_rhovim )
df_rhovim <- pivot_longer( df_rhovim, cols=c( "RHO", "VIM" ) )
df_rhovim

ggplot( df_rhovim ) +
  geom_histogram( aes( value ) ) +
  facet_wrap( vars(name) )

2. colData()

Introduce

• How to access the cell metadata?

help( colData )

meta <- colData( sce )

• What class of object is this?

class( meta )

[1] "DFrame"
attr(,"package")
[1] "S4Vectors"

• What is DFrame?

meta

DataFrame with 8217 rows and 18 columns
                   final_cluster_labels          libraryID   celltype
                               <factor>           <factor>   <factor>
AAACCCAGTCGGCACT-1                   9       fovea_donor_1 Unknown   
AAACGAACAGGGACTA-1                   15      fovea_donor_1 Glial-F1  
AAACGAAGTTCTCCAC-1                   3       fovea_donor_1 Cones-Fov 
AAACGAATCATAGCAC-1                   5       fovea_donor_1 Bipolar-1 
AAACGCTAGGAGGCAG-1                   8A      fovea_donor_1 Horizontal
...                                 ...                ...        ...
TTTGACTCATGACTGT-6                    5 peripheral_donor_3  Bipolar-1
TTTGACTGTCCGGACT-6                    6 peripheral_donor_3  Bipolar-2
TTTGACTTCTGAGTCA-6                    2 peripheral_donor_3  Rods-2   
TTTGGAGCAGTTTGGT-6                    6 peripheral_donor_3  Bipolar-2
TTTGGAGTCGTTCTAT-6                    6 peripheral_donor_3  Bipolar-2
                       region development_stage_ontology_term_id
                     <factor>                           <factor>
AAACCCAGTCGGCACT-1      fovea                     HsapDv:0000087
AAACGAACAGGGACTA-1      fovea                     HsapDv:0000087
AAACGAAGTTCTCCAC-1      fovea                     HsapDv:0000087
AAACGAATCATAGCAC-1      fovea                     HsapDv:0000087
AAACGCTAGGAGGCAG-1      fovea                     HsapDv:0000087
...                       ...                                ...
TTTGACTCATGACTGT-6 peripheral                     HsapDv:0000087
TTTGACTGTCCGGACT-6 peripheral                     HsapDv:0000087
TTTGACTTCTGAGTCA-6 peripheral                     HsapDv:0000087
TTTGGAGCAGTTTGGT-6 peripheral                     HsapDv:0000087
TTTGGAGTCGTTCTAT-6 peripheral                     HsapDv:0000087
                   development_stage disease_ontology_term_id
                            <factor>                 <factor>
AAACCCAGTCGGCACT-1 human adult stage            MONDO:0005041
AAACGAACAGGGACTA-1 human adult stage            MONDO:0005041
AAACGAAGTTCTCCAC-1 human adult stage            MONDO:0005041
AAACGAATCATAGCAC-1 human adult stage            MONDO:0005041
AAACGCTAGGAGGCAG-1 human adult stage            MONDO:0005041
...                              ...                      ...
TTTGACTCATGACTGT-6 human adult stage             PATO:0000461
TTTGACTGTCCGGACT-6 human adult stage             PATO:0000461
TTTGACTTCTGAGTCA-6 human adult stage             PATO:0000461
TTTGGAGCAGTTTGGT-6 human adult stage             PATO:0000461
TTTGGAGTCGTTCTAT-6 human adult stage             PATO:0000461
                              disease assay_ontology_term_id
                             <factor>               <factor>
AAACCCAGTCGGCACT-1 glaucoma (disease)            EFO:0009898
AAACGAACAGGGACTA-1 glaucoma (disease)            EFO:0009898
AAACGAAGTTCTCCAC-1 glaucoma (disease)            EFO:0009898
AAACGAATCATAGCAC-1 glaucoma (disease)            EFO:0009898
AAACGCTAGGAGGCAG-1 glaucoma (disease)            EFO:0009898
...                               ...                    ...
TTTGACTCATGACTGT-6             normal            EFO:0009898
TTTGACTGTCCGGACT-6             normal            EFO:0009898
TTTGACTTCTGAGTCA-6             normal            EFO:0009898
TTTGGAGCAGTTTGGT-6             normal            EFO:0009898
TTTGGAGTCGTTCTAT-6             normal            EFO:0009898
                               assay ethnicity_ontology_term_id
                            <factor>                   <factor>
AAACCCAGTCGGCACT-1 10x v3 sequencing                           
AAACGAACAGGGACTA-1 10x v3 sequencing                           
AAACGAAGTTCTCCAC-1 10x v3 sequencing                           
AAACGAATCATAGCAC-1 10x v3 sequencing                           
AAACGCTAGGAGGCAG-1 10x v3 sequencing                           
...                              ...                        ...
TTTGACTCATGACTGT-6 10x v3 sequencing                           
TTTGACTGTCCGGACT-6 10x v3 sequencing                           
TTTGACTTCTGAGTCA-6 10x v3 sequencing                           
TTTGGAGCAGTTTGGT-6 10x v3 sequencing                           
TTTGGAGTCGTTCTAT-6 10x v3 sequencing                           
                   ethnicity      sex tissue_ontology_term_id
                    <factor> <factor>                <factor>
AAACCCAGTCGGCACT-1   unknown     male          UBERON:0001786
AAACGAACAGGGACTA-1   unknown     male          UBERON:0001786
AAACGAAGTTCTCCAC-1   unknown     male          UBERON:0001786
AAACGAATCATAGCAC-1   unknown     male          UBERON:0001786
AAACGCTAGGAGGCAG-1   unknown     male          UBERON:0001786
...                      ...      ...                     ...
TTTGACTCATGACTGT-6   unknown   female          UBERON:0013682
TTTGACTGTCCGGACT-6   unknown   female          UBERON:0013682
TTTGACTTCTGAGTCA-6   unknown   female          UBERON:0013682
TTTGGAGCAGTTTGGT-6   unknown   female          UBERON:0013682
TTTGGAGTCGTTCTAT-6   unknown   female          UBERON:0013682
                                        tissue
                   "realize"
AAACCCAGTCGGCACT-1             fovea centralis
AAACGAACAGGGACTA-1             fovea centralis
AAACGAAGTTCTCCAC-1             fovea centralis
AAACGAATCATAGCAC-1             fovea centralis
AAACGCTAGGAGGCAG-1             fovea centralis
...                                        ...
TTTGACTCATGACTGT-6 peripheral region of retina
TTTGACTGTCCGGACT-6 peripheral region of retina
TTTGACTTCTGAGTCA-6 peripheral region of retina
TTTGGAGCAGTTTGGT-6 peripheral region of retina
TTTGGAGTCGTTCTAT-6 peripheral region of retina
                   cell_type_ontology_term_id
                                     <factor>
AAACCCAGTCGGCACT-1                           
AAACGAACAGGGACTA-1                 CL:0000125
AAACGAAGTTCTCCAC-1                           
AAACGAATCATAGCAC-1                 CL:0000103
AAACGCTAGGAGGCAG-1                 CL:0000745
...                                       ...
TTTGACTCATGACTGT-6                 CL:0000103
TTTGACTGTCCGGACT-6                 CL:0000103
TTTGACTTCTGAGTCA-6                 CL:0000604
TTTGGAGCAGTTTGGT-6                 CL:0000103
TTTGGAGTCGTTCTAT-6                 CL:0000103
                                   cell_type    ident
                                    <factor> <factor>
AAACCCAGTCGGCACT-1 unknown                      local
AAACGAACAGGGACTA-1 glial cell                   local
AAACGAAGTTCTCCAC-1 foveal cone photoreceptor    local
AAACGAATCATAGCAC-1 bipolar neuron               local
AAACGCTAGGAGGCAG-1 retina horizontal cell       local
...                                      ...      ...
TTTGACTCATGACTGT-6          bipolar neuron      local
TTTGACTGTCCGGACT-6          bipolar neuron      local
TTTGACTTCTGAGTCA-6          retinal rod cell    local
TTTGGAGCAGTTTGGT-6          bipolar neuron      local
TTTGGAGTCGTTCTAT-6          bipolar neuron      local

data.frame( meta )

• How large is it?

dim( meta )

[1] 8217   18

• What categories exist?

colnames( meta )

 [1] "final_cluster_labels"               "libraryID"                         
 [3] "celltype"                           "region"                            
 [5] "development_stage_ontology_term_id" "development_stage"                 
 [7] "disease_ontology_term_id"           "disease"                           
 [9] "assay_ontology_term_id"             "assay"                             
[11] "ethnicity_ontology_term_id"         "ethnicity"                         
[13] "sex"                                "tissue_ontology_term_id"           
[15] "tissue"                             "cell_type_ontology_term_id"        
[17] "cell_type"                          "ident"                             

Subset rows

• colData(), data.frame(), rownames_to_column()

df_meta <- colData( sce ) %>%
  data.frame() %>%
  rownames_to_column()
df_meta

dim( df_meta )

[1] 8217   19

df_meta[ 1:10, 1:3 ]

Subset libraryID

table( df_meta$libraryID )

     fovea_donor_1      fovea_donor_2      fovea_donor_3 
               767               2078                733 
peripheral_donor_1 peripheral_donor_2 peripheral_donor_3 
              1715               1427               1497 

df_fovea2 <- df_meta %>%
  filter( libraryID=="fovea_donor_2" )
df_fovea2

dim( df_fovea2 )

[1] 2078   19

table( df_fovea2$cell_type )

                amacrine cell                bipolar neuron 
                            0                           583 
             endothelial cell     foveal cone photoreceptor 
                           48                           181 
                   glial cell               microglial cell 
                          574                            39 
                pericyte cell peripheral cone photoreceptor 
                           28                             0 
       retina horizontal cell         retinal ganglion cell 
                            8                           103 
             retinal rod cell                       unknown 
                          121                           393 

ggplot gene by cell_type

• full_join(), filter()

df_rhovim_fovea2 <- df_rhovim %>%
  full_join( df_meta, by="rowname" ) %>%
  filter( libraryID=="fovea_donor_2" )
df_rhovim_fovea2

• geom_violin(), facet_wrap(), theme()

ggplot( df_rhovim_fovea2 ) +
  geom_violin( aes( cell_type, value ) ) +
  facet_wrap( vars(name) ) +
  theme( axis.text.x=element_text( angle=90 ) )

3. reducedDim()

Introduce

• How to access the cluster visualization?

help( reducedDim )

• data.frame(), rownames_to_column()

df_umap <- reducedDim( sce, "UMAP" ) %>%
  data.frame() %>%
  rownames_to_column()
df_umap

ggplot cell_type

df_umap_meta <- full_join( df_umap, df_meta, by="rowname" )
df_umap_meta

ggplot( df_umap_meta ) +
  geom_point( aes( UMAP_1, UMAP_2, color=cell_type ) )

ggplot RHO

df_umap_rho <- full_join( df_umap, df_rhovim, by="rowname" ) %>%
  filter( name=="RHO" )
df_umap_rho

ggplot( df_umap_rho ) +
  geom_point( aes( UMAP_1, UMAP_2, color=value ) )

4. findMarkers()

Introduce

• How to access differentially expressed genes?

help( findMarkers )

• Preview of findMarkers() syntax

colLabels( sce ) <- colData( sce )$cell_type
markers <- findMarkers( sce )

• What class of object is this?

class( markers )

[1] "SimpleList"
attr(,"package")
[1] "S4Vectors"

• How large is it?

length( markers )

[1] 12

• Build up to data.frame( markers[[1]] )

markers

List of length 12
names(12): amacrine cell bipolar neuron ... retinal rod cell unknown

hist summary.logFC

markers[1]

List of length 1
names(1): amacrine cell

df_markers_amacrine <- markers[[1]] %>%
  data.frame()
df_markers_amacrine

hist( df_markers_amacrine$summary.logFC )

plotUMAP markers

df_markers_amacrine %>%
  filter( summary.logFC > 0 )

help( plotUMAP )

plotUMAP( sce, colour_by="MEG3", text_by="cell_type" )

Exercise 2

1. Recreate plotUMAP() using ggplot2
   • Plot another gene, cell metadata
2. Recreate plotExpression() using ggplot2
   • Plot another gene by cell metadata
3. Adjust ggplot2 appearance
   • Add titles, customize x/y labels, theme
4. Explore findMarkers() results
   • Use plotUMAP(), plotExpression()

Document software

sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 8

Matrix products: default
BLAS:   /home/idies/R/lib64/R/lib/libRblas.so
LAPACK: /home/idies/R/lib64/R/lib/libRlapack.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils    
[7] datasets  methods   base     

other attached packages:
 [1] forcats_0.5.1               stringr_1.4.0              
 [3] dplyr_1.0.7                 purrr_0.3.4                
 [5] readr_1.4.0                 tidyr_1.1.3                
 [7] tibble_3.1.1                tidyverse_1.3.1            
 [9] scran_1.18.7                scater_1.18.6              
[11] ggplot2_3.3.3               SingleCellExperiment_1.12.0
[13] SummarizedExperiment_1.20.0 Biobase_2.50.0             
[15] GenomicRanges_1.42.0        GenomeInfoDb_1.26.7        
[17] IRanges_2.24.1              S4Vectors_0.28.1           
[19] BiocGenerics_0.36.1         MatrixGenerics_1.2.1       
[21] matrixStats_0.58.0          SeuratObject_4.0.0         
[23] Seurat_4.0.1               

loaded via a namespace (and not attached):
  [1] readxl_1.3.1              backports_1.2.1          
  [3] plyr_1.8.6                igraph_1.2.6             
  [5] lazyeval_0.2.2            splines_4.0.3            
  [7] BiocParallel_1.24.1       listenv_0.8.0            
  [9] scattermore_0.7           digest_0.6.27            
 [11] htmltools_0.5.1.1         viridis_0.6.1            
 [13] fansi_0.4.2               magrittr_2.0.1           
 [15] tensor_1.5                cluster_2.1.0            
 [17] ROCR_1.0-11               limma_3.46.0             
 [19] globals_0.14.0            modelr_0.1.8             
 [21] spatstat.sparse_2.0-0     colorspace_2.0-0         
 [23] rvest_1.0.0               ggrepel_0.9.1            
 [25] haven_2.4.1               xfun_0.22                
 [27] crayon_1.4.1              RCurl_1.98-1.3           
 [29] jsonlite_1.7.2            spatstat.data_2.1-0      
 [31] survival_3.2-7            zoo_1.8-9                
 [33] glue_1.4.2                polyclip_1.10-0          
 [35] gtable_0.3.0              zlibbioc_1.36.0          
 [37] XVector_0.30.0            leiden_0.3.7             
 [39] DelayedArray_0.16.3       BiocSingular_1.6.0       
 [41] future.apply_1.7.0        abind_1.4-5              
 [43] scales_1.1.1              edgeR_3.32.1             
 [45] DBI_1.1.1                 miniUI_0.1.1.1           
 [47] Rcpp_1.0.6                viridisLite_0.4.0        
 [49] xtable_1.8-4              dqrng_0.2.1              
 [51] reticulate_1.19           spatstat.core_2.1-2      
 [53] rsvd_1.0.5                htmlwidgets_1.5.3        
 [55] httr_1.4.2                RColorBrewer_1.1-2       
 [57] ellipsis_0.3.2            ica_1.0-2                
 [59] farver_2.1.0              pkgconfig_2.0.3          
 [61] scuttle_1.0.4             dbplyr_2.1.1             
 [63] sass_0.3.1                uwot_0.1.10              
 [65] deldir_0.2-10             locfit_1.5-9.4           
 [67] utf8_1.2.1                labeling_0.4.2           
 [69] tidyselect_1.1.1          rlang_0.4.10             
 [71] reshape2_1.4.4            later_1.2.0              
 [73] cellranger_1.1.0          munsell_0.5.0            
 [75] tools_4.0.3               cli_2.5.0                
 [77] generics_0.1.0            broom_0.7.6              
 [79] ggridges_0.5.3            evaluate_0.14            
 [81] fastmap_1.1.0             yaml_2.2.1               
 [83] goftest_1.2-2             fs_1.5.0                 
 [85] knitr_1.33                fitdistrplus_1.1-3       
 [87] RANN_2.6.1                pbapply_1.4-3            
 [89] future_1.21.0             nlme_3.1-149             
 [91] sparseMatrixStats_1.2.1   mime_0.10                
 [93] xml2_1.3.2                rstudioapi_0.13          
 [95] compiler_4.0.3            beeswarm_0.3.1           
 [97] plotly_4.9.3              png_0.1-7                
 [99] spatstat.utils_2.1-0      reprex_2.0.0             
[101] statmod_1.4.35            bslib_0.2.4              
[103] stringi_1.5.3             bluster_1.0.0            
[105] lattice_0.20-41           Matrix_1.2-18            
[107] vctrs_0.3.8               pillar_1.6.0             
[109] lifecycle_1.0.0           spatstat.geom_2.1-0      
[111] lmtest_0.9-38             jquerylib_0.1.4          
[113] RcppAnnoy_0.0.18          BiocNeighbors_1.8.2      
[115] data.table_1.14.0         cowplot_1.1.1            
[117] bitops_1.0-7              irlba_2.3.3              
[119] httpuv_1.6.0              patchwork_1.1.0.9000     
[121] R6_2.5.0                  promises_1.2.0.1         
[123] KernSmooth_2.23-17        gridExtra_2.3            
[125] vipor_0.4.5               parallelly_1.25.0        
[127] codetools_0.2-16          MASS_7.3-53              
[129] assertthat_0.2.1          withr_2.4.2              
[131] sctransform_0.3.2         GenomeInfoDbData_1.2.4   
[133] hms_1.0.0                 mgcv_1.8-33              
[135] grid_4.0.3                rpart_4.1-15             
[137] beachmat_2.6.4            rmarkdown_2.7            
[139] DelayedMatrixStats_1.12.3 Rtsne_0.15               
[141] lubridate_1.7.10          shiny_1.6.0              
[143] ggbeeswarm_0.6.0