Inconsistency between CHOIR clusters with Harmony Integration

[alifarhat40]

Hello,

I am trying to apply CHOIR on the RNA assay of a multiomic dataset. I did this in two ways. First, I applied Harmony integration using the harmony package followed by seurat clustering with resolution of 2 to be conservative. The second way I used the following command in CHOIR:

cellranger <- readRDS("integrated_cellranger_peaks_RNA.rds")

cellranger

DefaultAssay(cellranger) <- "RNA"
cellranger

cellranger <- CHOIR(cellranger, 
                batch_correction_method = "Harmony",
                batch_labels = "sample",
                use_assay = "RNA",
                n_cores= 50)

However, the clusters do not look good when compared to seurat's clustering. See comparison below:

1. How is pca used when CHOIR says computing pca? Does it make sense to start the tree with pca? Shouldn't CHOIR use the dimensionality reduction obtained from Harmony instead to build the tree and the graphs?

2. Should I add use_slot = counts or = data?

3. I would like to understand more how the dimensionality reduction works when setting subtree_reductions = TRUE .

Warning message in .validInput(subtree_reductions, "subtree_reductions", list(reduction, :
“Supplied dimensionality reduction matrix for parameter 'subtree_reductions' will only be used for the root tree. Thereafter, the dimensionality reductions for each subtree will be calculated according to the specified 'reduction_method'. To use only the supplied dimensionality reduction matrix, set parameter 'subtree_reductions' to FALSE.”

How are subsequent dimensionality reductions computed? It is my understanding that Harmony only provides a reduced dimensionality reduciton, and thats it. It does not provide a corrected counts or corrected data matrix. So do you compute PCA every time on the counts matrix for subsequent trees? If so, where does Harmony integration come to play here? I suspect the discrepancy is coming from the buildTree() function. These integration methods do not provide corrected counts or corrected scaled data matrices, so how does integration and batch correction come into play here.

Similarly, I get errors and issues with ATAC Harmony and Multiomic Harmony integrations. Also, should I be using counts or data below in the use_slot parameter?

cellranger <- CHOIR(cellranger,
                    use_assay = c("peaks", "RNA"),
                    use_slot = c("data", "data"),
                    atac = c(TRUE, FALSE),
                    batch_correction_method = c("Harmony", "Harmony"),
                    batch_labels = c("sample","sample"),
                    distance_approx = FALSE,
                    n_cores= 50)

Note that I have 8 different batches.

[cathrinesant]

How is pca used when CHOIR says computing pca? Does it make sense to start the tree with pca? Shouldn't CHOIR use the dimensionality reduction obtained from Harmony instead to build the tree and the graphs?

CHOIR does use the dimensionality reduction obtained from Harmony to build the tree/graphs. PCA is computed, followed by Harmony, and the Harmony-corrected dimensionality reduction is used for every downstream step.

Should I add use_slot = counts or = data?

CHOIR has been most thoroughly tested with log normalized data, which is typically stored in slot/layer "data" for Seurat objects after Seurat::NormalizeData() has been run.

I would like to understand more how the dimensionality reduction works when setting subtree_reductions = TRUE

So do you compute PCA every time on the counts matrix for subsequent trees? If so, where does Harmony integration come to play here?

Rundown of how subtrees play a role in building/pruning the tree:

First, a set of variable features P0_var_features and a dimensionality reduction P0_reduction are generated using all cells. The P0_reduction is used to generate the "root tree" that starts with a single cluster composed of all cells, and branches until the mean silhouette score is maximized. At this point, we'll call the clusters at the bottom of the root tree, with the maximum silhouette score, "parent clusters". These will be numbered P1, P2, P3, and so on. Often, parent clusters will reflect the major cell types present in the dataset.

Now, the cells belonging to each parent cluster will be subsetted out in turn. When subtree_reductions = TRUE, a new set of variable features P1_var_features and a new dimensionality reduction P1_reduction are generated using only the cells belonging to parent cluster P1. PCA is used here by default, again followed by Harmony. The P1_reduction is then used to generate a subtree, which extends until the clusters are demonstrably overclustered.

This is repeated for each parent cluster, and the resulting subtrees are then all stitched together along with the root tree to form the final clustering tree that function CHOIR::buildTree() outputs. When pruning this tree using CHOIR::pruneTree(), the root/subtree information is used as follows:

• If two clusters being compared belong to the same subtree (e.g., both composed of cells from parent cluster P1), then the P1_reduction will be used to determine the distance between the clusters (and whether distance filters should be imposed via parameter distance_awareness), the P1_graph_nn will be used to determine the number of nearest neighbors shared between the two clusters (and whether adjacency filters should be imposed via parameter min_connections), and the features used to train the random forest classifier will come from P1_var_features (unless parameter feature_set = "all").
• If two clusters being compared belong to different subtrees (e.g., one is composed of cells from parent cluster P2 and the other composed of cells from parent cluster P3) or either/both are composed of a mix of cells from various parent clusters, then P0_reduction will be used to determine the distance between the clusters, the P0_graph_nn will be used to determine the number of nearest neighbors shared between the two clusters, and the features used to train the random forest classifier will come from P0_var_features.

This allows CHOIR to use more nuanced information when comparing similar clusters, helping to distinguish cell states or cell subtypes. For example, in the brain, a set of variable features generated from the entire set of cells might not include the genes that best distinguish different interneuron subtypes or microglial inflammatory states, if these genes are relatively lowly expressed overall. But in a set of variable features generated only from the cells belonging to a certain parent cluster, these genes are more likely to make the cut.

So, by setting subtree_reductions = FALSE, you forfeit this more nuanced information, and the P0_var_features and P0_reduction are used for all comparisons. This can be totally fine for certain datasets, but for others, you may miss valuable differences.

CHOIR also uses batch information within each permutation test comparison. See #28 (comment) for a recent explainer of some of the ways this can get complicated as the number of batches increases, particularly when each batch represents a single sample. Note that CHOIR's batch correction was designed to work best when each batch contains more than 1 sample. Trying out the new parameter in the dev branch, max_n_batch, as mentioned in #28 (comment), might help with your odd clustering results here.

Similarly, I get errors and issues with ATAC Harmony and Multiomic Harmony integrations. Also, should I be using counts or data below in the use_slot parameter?

CHOIR for datasets containing ATAC-seq has been most thoroughly tested using ArchR objects. I'm not sure whether counts/data makes the most sense for a peak matrix. You could also try the dev branch though, because I think I recently resolved the issue preventing distance_approx = TRUE from working with Signac object (no promises).

[alifarhat40]

Thank you @catpetersen for your thoughtful responses. I have been using this method ever since it came out and it never disappoints.

1. So to make sure I understand you correctly: If I have an experiment where I do multiome (atac + rna) sequencing on cells at different days: example day1, day 30, day 60, day 90, then using batch correction on these cells should be okay with the batch being assigned to the day collected assuming they are the same type of cells, right?

2. I am unclear what max_n_batch is supposed to represent here. If Hamrony does batch correction, then why do you say:
   "Note that CHOIR's batch correction was designed to work best when each batch contains more than 1 sample."
   

I thought a sample is a batch. In my data, the batches are under "sample" hence why I do batch_labels = "sample". Each sample is data collected at a different timepoint. To me, I thought seurat clustering on harmony should at least match the clustering obtained from CHOIR (using the Harmony option). Sure, CHOIR might fail and merge additional clusters than it should, but at least they should match, but here it seems all mushy and not using the neighbors or communities the right way?

3. In the main branch, I tried with ArchR on this dataset I described in (1) and it seemed to have worked, but I had to remove the vector of batch_correction_method and just mention it as below since it threw out error when I followed the guildeline on the main CHOIR page. Here is my command:

proj <- CHOIR(
    proj,
    ArchR_matrix = c("GeneScoreMatrix", "GeneExpressionMatrix"),
    ArchR_depthcol = c("nFrags", "Gex_nUMI"),
    atac = c(TRUE, FALSE),
    batch_correction_method = "Harmony",
    batch_labels = "sample",
    n_cores = 50
)

But I do have a concern. In my logfile, it said it performed IterativeLSI with 2000 variable features, and I am assuming this is for the RNA gene expression. shouldn't it use PCA for the gene expression? Why is it using IterativeLSI?

[cathrinesant]

Thanks, that's kind of you to say!

So to make sure I understand you correctly: If I have an experiment where I do multiome (atac + rna) sequencing on cells at different days: example day1, day 30, day 60, day 90, then using batch correction on these cells should be okay with the batch being assigned to the day collected assuming they are the same type of cells, right?

If the day here represents days in culture, or some such experimental variable, then I would worry that batch correction using "day collected" as the batch will confound true biological differences due to time in culture with possible batch effects. CHOIR assumes that cell groups of interest, such as unique cell types or states, are not batch-confounded. This sounds like it may be an impractical assumption for your experiment, in which case, CHOIR's batch correction may erroneously merge some clusters. If days is not confounded with any relevant biological or experimental variable, and just represents technical replicates collected on different days, then the above does not apply.

Regardless, you can always look through the records object@misc$CHOIR$records$comparison_records for more details about each pairwise cluster comparison. Your data looks like it has a pretty continuous gradient, rather than several clearly distinct cell populations, which also makes things more difficult.

I am unclear what max_n_batch is supposed to represent here.

See #28 (comment) for a detailed description of parameter max_n_batch. It's still in active development.

In my logfile, it said it performed IterativeLSI with 2000 variable features, and I am assuming this is for the RNA gene expression. shouldn't it use PCA for the gene expression? Why is it using IterativeLSI?

Since the next step is to generate a joint dimensionality reduction by combining the ATAC reduction and the RNA reduction, we felt it was more reasonable to use the same dimensionality reduction method by default. This is also the suggested procedure when using ArchR without CHOIR for multi-omic analysis.

[alifarhat40]

@catpetersen I see. My dataset is just one mouse organ at different ages, we are trying to determine the genetic drivers of aging. So its a few subtypes of the same celltype. One sample is at 1 month, another at 2 months, then 8 months, then 18 months. Nothing else confounds them, just the age of development. In that case, can I use CHOIR with Harmony by setting the batches to my "age sample"? Or does it break down as you explain above?

[alifarhat40]

So I tried running the following on my multiome (scRNa + scATAC) ArchR object:

proj <- CHOIR(proj,
              ArchR_matrix = c("GeneScoreMatrix", "GeneExpressionMatrix"),
              ArchR_depthcol = c("nFrags", "Gex_nUMI"),
              atac = c(TRUE, FALSE),
              batch_correction_method = "Harmony",
              batch_labels = "sample",
              n_cores= 50)

But the results do not look good

Compare it to the WNN_Harmony integration from here with resolution 2:

This issue seems consistent with what others are facing with integrations and batches. However, I will try it once more without batch correction and see if it gives any better results.