Score Samples

The first step is to score all samples/cells by using specified approach. The score can be composed of 3 terms: TF (term/feature frequency), IDF (inverse document/cell frequency) and IAE (inverse average expression of features). Each term has a couple of available choices with different formats to suit labeled or un-labeled data. Users can use function idf_iae_methods() to see available methods for IDF/IAE term. More details of each term can be seen in help page of each function, e.g. ?idf.

## show available methods
idf_iae_methods()
#>          unlabel HDBSCAN                label IGM              unlabel max 
#>                    "hdb"                    "igm"                      "m" 
#>                     null        label probability label relative frequency 
#>                   "null"                   "prob"                     "rf" 
#>               unlabel SD         unlabel standard 
#>                     "sd"               "standard"

The basic version of TF, IDF and IAE can be termed as:

$\mathbf{TF_{i,j}}=\frac{N_{i,j}}{\sum_j{N_{i,j}}},$ $\mathbf{IDF_i} = \log(1+\frac{n}{n_i+1}),$ $\mathbf{IAE_i} = \log(1+\frac{n}{\sum_j^n\hat N_{i,j}+1})$

Where N_i, j is the counts of feature i in cell j; N̂_i, j is max (0, N_i, j − threshold); n is the total number of documents(cells); n_i is $\sum_{j = 1}^{n} \mathrm{sign}(N_{i,j} > \mathrm{threshold})$.

Here for labeled data, we can choose logTF * IDF_prob * IAE_prob for marker identification: score = log TF * IDF_prob * IAE_prob

The probability version of IDF can be termed as: $$\mathbf{IDF_{i,j}} = \log(1+\frac{\frac{n_{i,j\in D}}{n_{j\in D}}}{\max(\frac{n_{i,j\in \hat D}}{n_{j\in \hat D}})+ e^{-8}}\frac{n_{i,j\in D}}{n_{j\in D}})$$

And the probability version of IAE can be termed as: $$\mathbf{IAE_{i,j}} = \log(1+\frac{\mathrm{mean}(\hat N_{i,j\in D})}{\max(\mathrm{mean}(\hat N_{i,j\in \hat D}))+ e^{-8}}*\mathrm{mean}(\hat N_{i,j\in D}))$$

Where D is the category of cell j; D̂ is the category other than D.

TF here stands for gene frequency, which is similar to CPM, while IDF represents the inverse cell/sample frequency for scRNA-seq data, and IAE is the inverse average expression of each gene across all cells or cells in each labeled group.

Another advantage of smartid is that it can start with raw counts data, with no need for pre-processed data. And the scoring is quite fast.

## compute score
system.time(
  data_sim <- cal_score(
    data_sim,
    tf = "logtf",
    idf = "prob",
    iae = "prob",
    par.idf = list(label = "Group"),
    par.iae = list(label = "Group")
  )
)
#>    user  system elapsed 
#>   0.247   0.004   0.251

## score and tf,idf,iae all saved
assays(data_sim)
#> List of length 7
#> names(7): BatchCellMeans BaseCellMeans BCV CellMeans TrueCounts counts score
names(metadata(data_sim))
#> [1] "Params"     "up_markers" "tf"         "idf"        "iae"

Scale and Transform Score

Scaling is needed to find the markers specific to the group, however, standard scaling might fail due to the rare populations. Here smartid uses a special scaling strategy scale_mgm(), which can scale imbalanced data by given group labels. By doing this, we can avoid the bias towards features with larger numerical ranges during feature selection.

The scale method is depicted as below:

$$z=\frac{x-\frac{\sum_k^{n_D}(\mu_k)}{n_D}}{sd}$$

The score will be transformed using softmax before passing to EM algorithm.

top_m <- top_markers(
  data = data_sim,
  label = "Group",
  n = Inf # set Inf to get all features processed score
)
top_m
#> # A tibble: 4,000 × 3
#> # Groups:   .dot [4]
#>    .dot   Genes    Scores
#>    <chr>  <chr>     <dbl>
#>  1 Group1 Gene352 0.0216 
#>  2 Group1 Gene696 0.0200 
#>  3 Group1 Gene425 0.0171 
#>  4 Group1 Gene391 0.0160 
#>  5 Group1 Gene225 0.0159 
#>  6 Group1 Gene893 0.0147 
#>  7 Group1 Gene913 0.0113 
#>  8 Group1 Gene650 0.0107 
#>  9 Group1 Gene172 0.00998
#> 10 Group1 Gene547 0.00980
#> # ℹ 3,990 more rows

The top n features for each group will be ordered and listed in top_m. smartid provides easy-to-use functions to visualize top feature scores in each group and compare with actual up-regulated DEGs.

It’s clear that the real UP DEGs are popping up to the top n features. And for the negative control “Group 4”, the shape of top feature score is totally different from the ones with DEGs, which can provide more insights to help understand the data.

score_barplot(
  top_markers = top_m,
  column = ".dot",
  f_list = slot(data_sim, "metadata")$up_markers,
  n = 20
)

As we can see, there is an UP DEG ‘Gene76’ not popping up in Group 2, we can check the relative expression of this gene using violin plot. It is clear that this gene is not significantly highly expressed in Group 2 and the average expression is quite low across all cells.

This can also be confirmed in data simulation information, where the scale factor is higher in Group2, but the GeneMean is too small to be confident. Thus this gene won’t be selected by smartid.

sin_score_boxplot(
  metadata(data_sim)$tf,
  features = "Gene76",
  ref.group = "Group2",
  label = data_sim$Group
)
#> Warning: `add_rownames()` was deprecated in dplyr 1.0.0.
#> ℹ Please use `tibble::rownames_to_column()` instead.
#> ℹ The deprecated feature was likely used in the smartid package.
#>   Please report the issue at
#>   <https://github.com/DavisLaboratory/smartid/issues>.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

## sim gene info
SummarizedExperiment::elementMetadata(data_sim)[76, ]
#> DataFrame with 1 row and 8 columns
#>          Gene BaseGeneMean OutlierFactor   GeneMean DEFacGroup1 DEFacGroup2
#>   <character>    <numeric>     <numeric>  <numeric>   <numeric>   <numeric>
#> 1      Gene76   0.00628427             1 0.00628427           1     1.73743
#>   DEFacGroup3 DEFacGroup4
#>     <numeric>   <numeric>
#> 1           1           1

Marker Selection

As we can see from above, there is a distinctly different distribution of feature score between group with DEGs and without DEGs. And there is a clear separation (break point) between the real DEGs and non-DEGs.

To help automatically select real markers for each group, smartid used an expectation maximization (EM) approach to identify which genes fall into the real DEGs distribution and which are not.

Regarding the distribution of scores as a mixture model, here we can choose function markers_mixmdl() in smartid to separate features. There are 2 available mixture model to choose: normal (Gaussian) or gamma. We choose “norm” here as it runs faster.

smartid also allows to plot the mixture distribution plot after EM. It’s obvious that the top 2 components of Group 4 share quite similar distribution, thus no markers will be selected for this group.

set.seed(123)
marker_ls <- markers_mixmdl(
  top_markers = top_m,
  column = ".dot",
  ratio = 2,
  dist = "norm",
  plot = TRUE
)
#> number of iterations= 88
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#> number of iterations= 81
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#> number of iterations= 36
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#> number of iterations= 507
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

marker_ls
#> $Group1
#>  [1] "Gene352" "Gene696" "Gene425" "Gene391" "Gene225" "Gene893" "Gene913"
#>  [8] "Gene650" "Gene172" "Gene547" "Gene436" "Gene904" "Gene312" "Gene308"
#> 
#> $Group2
#> [1] "Gene419" "Gene332" "Gene141" "Gene37"  "Gene682" "Gene778" "Gene818"
#> 
#> $Group3
#> [1] "Gene833" "Gene357" "Gene518" "Gene26"  "Gene405" "Gene462" "Gene28" 
#> 
#> $Group4
#> NULL

We can also compare our selected markers with real DEGs. As there is no markers or DEG in group 4, only show overlap from Group1-3. It’s clear that all the markers identified by smartid are real DEGs, with a couple of missing genes in Group 1 and 2. But as what we showed above, those genes only exhibit low mean expression across all cells, thus not confident enough to be selected as markers.

library(UpSetR)

upset(fromList(c(slot(data_sim, "metadata")$up_markers, marker_ls)), nsets = 6)

smartid also provides some other implementation of marker selection. Here is another example using mclust. Different from markers_mixmdl(), markers_mclust() doesn’t need a pre-defined number of components (which is default 3 in markers_mixmdl()), instead, it will select the number of components by searching a series of potential numbers. This method is sometimes more robust than markers_mixmdl() but will also take longer running time.

Similarly, this method also allows to plot the mixture distribution for each component, but separately.

set.seed(123)
marker_ls_new <- markers_mclust(
  top_markers = top_m,
  column = ".dot",
  method = "max.one",
  plot = TRUE
)

names(marker_ls_new) <- paste(names(marker_ls_new), "new")

We can compare the marker list with the previous one. The overlap result shows 2 methods can achieve almost the same marker list for each group.

upset(fromList(c(marker_ls, marker_ls_new)), nsets = 6)

Score Samples

Similarly, the first step is to score samples/cells using the specified method. This step also starts with raw counts data, without need for data pre-processing, which is quite convenient and fast.

## compute score without label
system.time(
  data_sim <- cal_score(
    data_sim,
    tf = "logtf",
    idf = "sd",
    iae = "sd",
    new.slot = "score_unlabel"
  )
)
#>    user  system elapsed 
#>   0.182   0.017   0.206

## new score is saved and tf,idf,iae all updated
assays(data_sim)
#> List of length 8
#> names(8): BatchCellMeans BaseCellMeans BCV ... counts score score_unlabel
names(metadata(data_sim))
#> [1] "Params"     "up_markers" "tf"         "idf"        "iae"

Compute Overall Score for Gene-set

To compare overall score of the given gene-set, we don’t need to scale and transform score this time. Using gs_score() can easily compute the overall score for each cell based on the given gene-set list.

## compute score for each group marker list
data_sim <- gs_score(
  data = data_sim,
  features = marker_ls[1:3], # group 4 has no markers
  slot = "score_unlabel",
  suffix = "score.unlabel" # specify the suffix of names to save
)

## saved score
colnames(colData(data_sim))
#> [1] "Cell"                 "Batch"                "Group"               
#> [4] "ExpLibSize"           "Group1.score.unlabel" "Group2.score.unlabel"
#> [7] "Group3.score.unlabel"

Now we get 3 columns of score for each group markers. We can then visualize the score across groups, see how well it can discern the target group.

It’s evident that the score can sufficiently separate the target group from all others.

as.data.frame(colData(data_sim)) |>
  tidyr::pivot_longer("Group1.score.unlabel":"Group3.score.unlabel",
    names_to = "group markers",
    values_to = "score"
  ) |>
  ggplot(aes(x = Group, y = score, fill = Group)) +
  geom_boxplot() +
  facet_wrap(~`group markers`, scales = "free") +
  theme_bw()