This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Can't reproduce harmony results in Seurat

Hello,

I integrated some scRNA-seq data using Harmony on Seurat V4 back in 2022, but I didn’t save the Seurat object at that time. When I needed to perform further analysis, I reran Harmony on the same data, but the clustering results were quite different. Specifically, the number of clusters decreased from 20 to 18, and most importantly, I can’t find cluster 17, which is crucial for our analysis.

enter image description here enter image description here

Is there any way to recover the clustering results from 2022? If reproducing the previous results isn't possible, I plan to increase the resolution parameter until I see the cluster of interest and annotate the new clusters. How can I justify that the new annotation closely matches the previous one?

Thank you for your help!

seurat scrna-seq harmony
• 5,149 views • link
written by fifty_fifty • fifty_fifty

Are you using the exact same computer (same processor and OS) ?

• 0 views • link
written by Bastien Hervé Bastien Hervé

No, I have a diferent MacBook now. Do you think that can affect harmony embeddings?

• 0 views • link
written by fifty_fifty • fifty_fifty

looks like my it's related to M2 processor. I ran the script on Linux and it was able to find the cluster of interest. Still, the number of clusters is 19 but I used Seurat V5, maybe if I go back to V4 I'll be able to reproduce my initial result. Thank you!

enter image description here

• 0 views • link
written by fifty_fifty • fifty_fifty

2 answers

Hi, I have in the past been able to reproduce old results using the most recent version of Seurat (v5). Can you share your code that you were (are?) using?

• 0 views • link
written by yura.grabovska • yura.grabovska

I am using the same code and data that I had 2 years ago

Kim_brain_data <- Kim_data[,grepl("NS_",colnames(Kim_data))]
Kim_brain <- CreateSeuratObject(Kim_brain_data, project="Kim_brain", min.cells=3, min.features = 200)
Kim_brain[["percent.mt"]]<-PercentageFeatureSet(Kim_brain, pattern="^MT-")
Kim_brain[["percent.rb"]]<-PercentageFeatureSet(Kim_brain, pattern="RP[SL]")
VlnPlot(Kim_brain, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), ncol = 4, pt.size = 0.3)
Kim_brain_sub <- subset(Kim_brain, subset = nFeature_RNA < 7500 & nFeature_RNA > 500 & nCount_RNA < 50000 & percent.mt < 15 & percent.rb < 50)
Kim_annot_brain_sub <- Kim_annot_brain[rownames(Kim_annot_brain) %in% colnames(Kim_brain_sub),]


Kim_brain_sub <- AddMetaData(Kim_brain_sub, Kim_annot_brain_sub)
Kim_brain_sub <- Seurat::NormalizeData(Kim_brain_sub, verbose = FALSE) %>%
  FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>% 
  ScaleData(verbose = FALSE) %>% 
  RunPCA(pc.genes = Kim_brain_sub@var.genes, npcs = 20, verbose = FALSE)

Kim_brain_sub <- Kim_brain_sub %>% 
  RunHarmony("Sample", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(Kim_brain_sub, 'harmony')

Kim_brain_sub <- Kim_brain_sub %>% 
  RunUMAP(reduction = "harmony", dims = 1:20) %>% 
  FindNeighbors(reduction = "harmony", dims = 1:20) %>% 
  FindClusters(resolution = 0.5) %>% 
  identity()
DimPlot(object = Kim_brain_sub, reduction = "umap", pt.size = .1, label = T)
• 0 views • link
written by fifty_fifty • fifty_fifty

When you run UMAP these days you get this:

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session

This is probably why your resulting UMAP looks different as well as the already mentioned difference in how the runs are seeded

• 0 views • link
written by yura.grabovska • yura.grabovska

thank you, I didn't notice that warning.

So, I ran it on Seurat 4.2.0 and got umap viz that doesn't look like the initial clusters either. Again, it found 18 clusters instead of 20.

enter image description here

here, I set.seed(1) in the very beginning before CreateSeuratObject. the visualization looks different but the number and structure of clusters did not change.

enter image description here

• 0 views • link
written by fifty_fifty • fifty_fifty

Can I just double check. You're basically saying you're getting a different number of clusters. But the general stucture looks similar. Do you know what the markers of cluster 17 were? Can you check for presense absense based on gene expression rather than visually looking for it?

• 0 views • link
written by yura.grabovska • yura.grabovska

The markers are expressed in two cell populations. One can be found with no problem but the second one is smaller. In my initial run from 2022, it was detected as a separate cluster 17

enter image description here

I still see those two cell populations in the new clustering visualization but the smaller cluster can't be detected as a separate cluster. It is a part of cluster 3 here

enter image description here

I technically can increase resolution until I see it but I was hoping that I would be able to reproduce my 2022 results. We've already published the 2022 results and need cell-cell interaction analysis with the same data.

• 0 views • link
written by fifty_fifty • fifty_fifty

Have you tried iterating over the seed with set.seed()? I doubt there is a way to directly infer the seed but you could probably identify your UMAP in a lineup.

Set.seed in Seurat

• 0 views • link
written by Jeremy Leipzig Jeremy Leipzig

IIRC the older versions of Seurat use the seed 42. But because methodology has changed between versions, especially how multiple samples are stored and normalised by default, it might not just be a seed issue. For example, Harmony as ran in a Seurat wrapper, and Harmony ran from the package have different input parameters.

• 0 views • link
written by yura.grabovska • yura.grabovska

I tried to change seed but it didn't affect the number or structure of clusters, it changed umap visualization. I am using the same script, same resolution parameter and getting less clusters. That's why I thought that maybe harmony creates different embeddings

• 0 views • link
written by fifty_fifty • fifty_fifty

Log in to answer this question.

More posts like this