On my wishlist but completely optional :3
So far unfortunately only in python.
Original publication
Github matlab implementation (original author), they mention R implementation but not done, author did not respond
Github python implementation 1 (original author)
Github python implementation 2("fast" version)

Just to keep track.
Will be needed for exploratory data analysis (e.g. unsupervised clustering) later.

After merge.HITObject there are only two samples in the matrix for x@aggregated_profile$pseudobulk$layer1$all
Should be all samples.

Just to keep track

In get.aggregated.profile there is a part where names are assigned to columns (for cell types). (see reference to code section below)
If I understand correctly:
The first part handles the case if there are multiple cell types. (?)
The second part if there is only one cell type.

Can't this be handled by a general loop, valid for ncol >=1 ?
I.e. like this

for(av in names(avg.exp)){ colnames(avg.exp[[i]]) <- unique([email protected][!is.na([email protected][[group.by.aggregated[[i]]]]), group.by.aggregated[[i]]]) }

Alternatively, shouldn't the first line read:

if(ncol(avg.exp[[i]]) > 1){
?

Just to keep track. As an alternative to silhouette. Other cluster scoring methods can be added in the future too.

Hello,

I am very interested in using your tools which seems great for cell classification. However I am currently using R version 4.1.3 and HiTME requires R >= 4.3.1. Would it be possible to make it compatible with a lower version of R?

Thanks

Noemie

When using get.cluster.score with batching, batch size should be considered. For example, if one batch has 30 observations and the other one only 3-5, then the latter will never produce a meaningful p-value, whereas the first prior one should be very meaningful. Fisher's method does not account for that and weights them equally. A better option would be the Weighted Z-Score Method.

Fisher's Method
Use Cases: Fisher's method is suitable when the p-values being combined are independent and come from studies testing the same null hypothesis.
Advantages: It is straightforward to implement and does not require assumptions about the distribution of the underlying test statistics.
Disadvantages: Fisher's method tends to lose power when the p-values being combined are not uniformly distributed under the null hypothesis.

Weighted Z-Score Method
Use Cases: Weighted Z-score method is useful when some studies are deemed more reliable or informative than others.
Advantages: Allows researchers to assign weights to individual studies based on their sample sizes, effect sizes, or other relevant factors, which can improve power and accuracy.
Disadvantages: Requires prior knowledge or assumptions about the relative reliability of individual studies, and incorrect weighting can bias results.

Hello,

Thanks again for creating this package. When using scGate alone, we have the option to choose whether the cell called as Multi should be considered as NA or not. This option is not included in Run.HiTME. Could you please include this parameter in the Run.HiTME function to be able to get the scGate multi cell.

Thanks

Hey Josep. I just noticed that when I Run.HiTME on a Seurat object, it messes up the data structure but only if I run ref.maps. Here is a minimal working example:

scGate_models_DB <- get_scGateDB(branch = "master", verbose = T, force_update = TRUE)
scGate_models <- scGate_models_DB$human$HiTME

layer1.links <- list("CD8" = "CL:0000625",
                     "CD4" = "CL:0000624",
                     "DC" = "CL:0000451",
                     "MoMac" = "CL:0000576_CL:0000235")

for(a in names(ref.maps)){
  ref.maps[[a]]@misc$layer1_link <- layer1.links[[a]]
}

# Take as input any small single sample that you have for testing
a <- subset(obj, subset = sample == obj$sample[1])

# This gives an error on running Run.HiTME a second time
b <- Run.HiTME(a, scGate.model = scGate_models, ref.maps = ref.maps)
c <- Run.HiTME(b, scGate.model = scGate_models, ref.maps = ref.maps)

# This gives no error on running Run.HiTME a second time
a <- subset(obj, subset = sample == obj$sample[1])
b <- Run.HiTME(a, scGate.model = scGate_models, ref.maps = NULL)
c <- Run.HiTME(b, scGate.model = scGate_models, ref.maps = NULL)

The error I get is:

Error in validObject(object = x) : 
  invalid class "Seurat" object: 1: all cells in assays must be in the same order as the Seurat object
invalid class "Seurat" object: 2: 'active.idents' must be named with cell names

I assume it might be due to subsetting the cells.

Hi Josep,

at one point during "Run.HiTME", the ref maps used get saved to obj@misc (see reference to specific line of code below).

However, the current behaviour only saves ref maps that were actually run, not all that were tried. For example:
I Run.HiTME on a list of seurat objects ("obj.list") with the CD4 and CD8 ref maps.
obj.list[[1]] has CD4 and CD8 cells, so obj.list[[1]]@misc$Projectils_param = c("CD4", "CD8")
obj.list[[2]] has CD4 but no CD8 cells, so obj.list[[2]]@misc$Projectils_param = c("CD4")

Wouldn't it be nicer to store all ref maps run on the object for traceability?
Alternatively, we could save both, e.g:

obj.list[[2]]@misc$Projectils_param$all <- c("CD4", "CD8")
obj.list[[2]]@misc$Projectils_param$executed <- c("CD4")

I point this out because it might be important to know if a cell type was not present, so the ref map was not executed.
PS: It would also be nice to store all cell types within the ref maps, either in @misc or in $functional.cluster as factor. This is in order to identify, which cell subtypes are NA (because layer 1 cell type did not map to layer 2) and which are 0 (layer 1 cell type was mapped to layer 2 but not detected in this sample) for the calculation of the cell type composition if some or all subtypes are equal to zero.

What do you think?

To drop a celltype from a HiTObject
Alternatively, the annotation can be removed from the Seurat object metadata column before running get.HiTObject

plot.celltype.freq currently operates on the HiTObject_list. I suggest that we do this from the HiT_summary resulting from merge.HiTObjects as this object uses less memory (in my case HiTObject_list (56 elements, 1GB) vs HiT_summary (391 MB)) and because downstream analyses are carried out with the HiT_summary.

It would also make it easier (simplifies the code) to create grouped boxplots, e.g. comparing responders vs non-responder celltype composition.

Would be nice to have a paramter freq.type to switch between absolute relative abundance and clr (to better visualize low abundant cell types).

NA values can be interesting. However, I would exclude NAs from the relative abundances as they add noise and make it difficult to interpret the ratios in the plots visually (i.e. keep according to what user specified in get.HiTObject useNA parameter, so by default exclude). I suggest we show number of NAs per sample in a separate plot, e.g. like this, right in the beginning (after Run.HiTME on the obj.list) as a quality control:

plot.nas <- function (obj.list = NULL,
                      annot.col = c("scGate_multi", "functional.cluster"),
                      bottom.mar = 10.2) {
  if (is.null(obj.list) &
      !is.list(obj.list) &
      !all(lapply(obj.list, inherits, "Seurat"))) {
    stop("Please provide a list of seurat objects")
  }
  
  for (col in annot.col) {
    na_perc_per_sample <- c()
    for (i in names(obj.list)) {
      if (col %in% names(obj.list[[i]]@meta.data)) {
        percs <- prop.table(table(obj.list[[i]]@meta.data[[col]], useNA = "ifany"))*100
        nas_perc <- unname(percs[is.na(names(percs))])
        na_perc_per_sample <- c(na_perc_per_sample, setNames(nas_perc, i))
      } else {
        stop(paste(col, " not found in obj.list item ", i))
      }
    }
    par(mar = c(bottom.mar, 4.1, 4.1, 2.1))
    barplot(na_perc_per_sample,
            ylab = "Percentage of NA values",
            las=2)
  }
}

plot.nas(obj.list)

Major issue. Make Run.HiTME preferably one line of code with default settings (scGate, ProjecTILs, etc.).

If there are NAs in the metadata column used for grouping, get.cluster.score fails. Add parameter to drop samples with NA (default). Alternatively, make NAs as.character, then as.factor the whole column.

Using get.cluster.score: with batching results in inconsistent data structure

Just to keep track