Labeling PBMC in the Seurat tutorial

One result of the tutorial analysis of the 3000 cell subset is a table of cell types and expression-based markers of cell identity. The first three columns of the table below are from concluding material in the Seurat tutorial; the remaining columns are created by “manual” matching between the Seurat terms and terms found in Cell Ontology.

kable(stab <- seur3kTab())

Relationships asserted in the Cell Ontology

Given the informally selected tags in the table above, we can sketch the Cell Ontology graph connecting the associated cell types. The ontoProc package adds functionality to ontologyPlot with make_graphNEL_from_ontology_plot. This allows use of all Rgraphviz and igraph visualization facilities for graphs derived from ontology structures.

Here we display the PBMC cell sets reported in the Seurat tutorial.

library(ontoProc)
cl = getOnto("cellOnto", "2023")
onto_plot2(cl, stab$tag)

Molecular features asserted in the Cell Ontology

The CLfeats function traces relationships and properties from a given Cell Ontology class. Briefly, each class can assert that it is the intersection_of other classes, and has_part, lacks_part, has_plasma_membrane_part, lacks_plasma_membrane_part can be asserted as relationships holding between cell type instances and cell components. The components are often cross-referenced to Protein Ontology or Gene Ontology. When the Protein Ontology component has a synonym for which an HGNC symbol is provided, that symbol is retrieved by CLfeats. Here we obtain the listing for a mature CD1a-positive dermal dendritic cell.

suppressMessages({
kable(CLfeats(cl, "CL:0002531", pr=pr, go=go))
})

The ctmarks function starts a shiny app that generates tables of this sort for selected cell types.

ctmarks snapshot

Mapping from gene ‘presence/role’ to cell type

The sym2CellOnto function helps find mention of given gene symbols in properties or parts of cell types.

kable(sdf <- as.data.frame(sym2CellOnto("ITGAM", cl, pr)))

table(sdf$cond)

## 
## lacksPMP 
##       22

kable(as.data.frame(sym2CellOnto("FOXP3", cl, pr)))

Use case: a set of cell types defined by “diagonal expression”

For a demonstration, we consider the discussion in Bakken et al. (2017), of a ‘diagonal’ expression pattern defining a group of novel cell types. A set of genes is identified and cells are distinguised by expressing exactly one gene from the set.

Diagonal expression pattern.

The necessary information is collected in a vector. The vector is the set of genes, the name of element i is the tag to be associated with the type of cell that expresses gene i and does not express any other gene in the set.

sigels = c("CL:X01"="GRIK3", "CL:X02"="NTNG1", "CL:X03"="BAGE2",
             "CL:X04"="MC4R", "CL:X05"="PAX6", "CL:X06"="TSPAN12", 
             "CL:X07"="hSHISA8", "CL:X08"="SNCG", "CL:X09"="ARHGEF28", 
             "CL:X10"="EGF")

A data.frame defining the cell types and their properties

The cyclicSigset function produces a data.frame instance connecting cell types with the genes expressed or unexpressed.

cs = cyclicSigset(sigels)
dim(cs)

## [1] 100   3

cs[c(1:5,9:13),]

##        gene   type     cond
## 1  ARHGEF28 CL:X09   hasExp
## 2     GRIK3 CL:X09 lacksExp
## 3     NTNG1 CL:X09 lacksExp
## 4     BAGE2 CL:X09 lacksExp
## 5      MC4R CL:X09 lacksExp
## 9      SNCG CL:X09 lacksExp
## 10      EGF CL:X09 lacksExp
## 11    BAGE2 CL:X03   hasExp
## 12    GRIK3 CL:X03 lacksExp
## 13    NTNG1 CL:X03 lacksExp

table(cs$cond)

## 
##   hasExp lacksExp 
##       10       90

It is expected that a tabular layout like this will suffice to handle general situations of cell type definition.

Translating the data.frame elements to OBO Term instances

The most complicated aspect of novel OBO term construction is the proper specifications of relationships with existing ontology components. A prolog that is mostly shared by all terms is generated programmatically for the diagonal pattern task.

 makeIntnProlog = function(id, ...) {
 # make type-specific prologs as key-value pairs
     c(  
       sprintf("id: %s", id),
       sprintf("name: %s-expressing cortical layer 1 interneuron, human", ...),
       sprintf("def: '%s-expressing cortical layer 1 interneuron, human described via RNA-seq observations' [PMID 29322913]", ...),
       "is_a: CL:0000099 ! interneuron",
       "intersection_of: CL:0000099 ! interneuron")
 }

The ldfToTerms API uses this to create a set of strings that can be parsed as a term.

pmap = c("hasExp"="has_expression_of", lacksExp="lacks_expression_of")
head(unlist(tms <- ldfToTerms(cs, pmap, sigels, makeIntnProlog)), 20)

##  [1] "[Term]"                                                                                                        
##  [2] "id: CL:X01"                                                                                                    
##  [3] "name: GRIK3-expressing cortical layer 1 interneuron, human"                                                    
##  [4] "def: 'GRIK3-expressing cortical layer 1 interneuron, human described via RNA-seq observations' [PMID 29322913]"
##  [5] "is_a: CL:0000099 ! interneuron"                                                                                
##  [6] "intersection_of: CL:0000099 ! interneuron"                                                                     
##  [7] "has_expression_of: PR:000008242 ! GRIK3"                                                                       
##  [8] "lacks_expression_of: PR:000011467 ! NTNG1"                                                                     
##  [9] "lacks_expression_of: PR:000004625 ! BAGE2"                                                                     
## [10] "lacks_expression_of: PR:000001237 ! MC4R"                                                                      
## [11] "lacks_expression_of: PR:000012318 ! PAX6"                                                                      
## [12] "lacks_expression_of: PR:000016738 ! TSPAN12"                                                                   
## [13] "lacks_expression_of: PR:B8ZZ34 ! hSHISA8"                                                                      
## [14] "lacks_expression_of: PR:000015325 ! SNCG"                                                                      
## [15] "lacks_expression_of: PR:000013942 ! ARHGEF28"                                                                  
## [16] "lacks_expression_of: PR:000006928 ! EGF"                                                                       
## [17] "[Term]"                                                                                                        
## [18] "id: CL:X02"                                                                                                    
## [19] "name: NTNG1-expressing cortical layer 1 interneuron, human"                                                    
## [20] "def: 'NTNG1-expressing cortical layer 1 interneuron, human described via RNA-seq observations' [PMID 29322913]"

The content in tms can then be appended to the content of the Cell Ontology cl.obo as text for import with ontologyIndex::get_OBO.

A data.frame mapping from informal to formal terms

Aaron Lun has produced a mapping from informal terms used in the Human Primary Cell Atlas to Cell Ontology tags. We provisionally include a copy of this mapping in ontoProc:

hpca_map = read.csv(system.file("extdata/hpca.csv", package="ontoProc"), strings=FALSE)
head(hpca_map)

##                          uncontrolled controlled
## 1        DC:monocyte-derived:immature CL:0000840
## 2      DC:monocyte-derived:Galectin-1 CL:0000451
## 3             DC:monocyte-derived:LPS CL:0000451
## 4                 DC:monocyte-derived CL:0000451
## 5 Smooth_muscle_cells:bronchial:vit_D CL:0002598
## 6       Smooth_muscle_cells:bronchial CL:0002598

We will rename columns of this map for convenience of our bind_formal_tags method.

names(hpca_map) = c("informal", "formal")  # obligatory for now

Binding formal tags to the HPCA data

I am turning this code off for now because there is no standard approach to getting the mapping from the SummarizedExperment yet. When SingleR merges the ‘standardized’ branch, this will come back.

Let’s retrieve the HPCA data from SingleR:

library(SummarizedExperiment)
library(SingleCellExperiment)
library(celldex)
hpca_sce = HumanPrimaryCellAtlasData()

Now bind the formal tags:

hpca_sce = bind_formal_tags(hpca_sce, "label.fine", hpca_map)
length(unique(hpca_sce$label.ont))

## [1] 66

We don’t check for failed mappings:

length(xx <- which(is.na(hpca_sce$label.ont)))

## [1] 2

if (length(xx)>0) print(colData(hpca_sce)[xx,])

## DataFrame with 2 rows and 3 columns
##            label.main             label.fine   label.ont
##           <character>            <character> <character>
## GSM154081     T_cells T_cell:CCR10+CLA+1,2..          NA
## GSM154084     T_cells T_cell:CCR10-CLA+1,2..          NA

sum(hpca_sce$label.ont == "", na.rm=TRUE) # iPS and BM

## [1] 50

Subsetting using the class hierarchy of Cell Ontology

cell_onto = ontoProc::getOnto("cellOnto", "2023")

## loading from cache

hpca_mono = subset_descendants( hpca_sce, cell_onto, "^monocyte$" )
table(hpca_mono$label.fine)

## 
##                          Monocyte                    Monocyte:CD14+ 
##                                27                                 3 
##                    Monocyte:CD16+                    Monocyte:CD16- 
##                                 6                                 7 
##                    Monocyte:CXCL4   Monocyte:F._tularensis_novicida 
##                                 2                                 6 
##                     Monocyte:MCSF Monocyte:S._typhimurium_flagellin 
##                                 2                                 1 
##             Monocyte:anti-FcgRIIB           Monocyte:leukotriene_D4 
##                                 2                                 4

table(hpca_mono$label.ont) # not much diversity

## 
## CL:0000576 CL:0001054 
##         57          3

hpca_tcell = subset_descendants( hpca_sce, cell_onto, "^T cell$" )
table(hpca_tcell$label.fine)

## 
##                    T_cell:CD4+              T_cell:CD4+_Naive 
##                             12                              6 
##     T_cell:CD4+_central_memory    T_cell:CD4+_effector_memory 
##                              5                              4 
##                    T_cell:CD8+     T_cell:CD8+_Central_memory 
##                             16                              3 
##    T_cell:CD8+_effector_memory T_cell:CD8+_effector_memory_RA 
##                              4                              4 
##              T_cell:CD8+_naive              T_cell:Treg:Naive 
##                              4                              2 
##                T_cell:effector             T_cell:gamma-delta 
##                              4                              2

table(hpca_tcell$label.ont) # 

## 
## CL:0000624 CL:0000625 CL:0000798 CL:0000895 CL:0000900 CL:0000904 CL:0000905 
##         12         16          2          6          4          5          4 
## CL:0000907 CL:0000911 CL:0000913 CL:0001062 CL:0002677 
##          3          4          4          4          2

uu = unique(hpca_tcell$label.ont)
onto_plot2(cell_onto, uu)

Visually identifying differences between ontology versions for a set of given terms

Using the ontoDiff function, you can visually identify which related terms exist in the newest provided version of an ontology but not in the old version.

cl = getOnto("diseaseOnto")

## downloading 1 resources

## retrieving 1 resource

## loading from cache

cl2 = getOnto(ontoname = "diseaseOnto", year_added = "2021")

## downloading 1 resources

## retrieving 1 resource

## loading from cache

cl3k = c("DOID:0040064","DOID:0040076","DOID:0081127","DOID:0081126","DOID:0081131","DOID:0060034")
ontoDiff(cl,cl2,cl3k)