This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Can I integrate two scRNAseq datasets using DE genes?

I am pretty new to the whole bioinformatics world and I got these two datasets that I need to integrate/analyse but they are giving me a hard time.

I am trying to integrate two very similar datasets. The hypothesis is that one is more mature than the other and I would like to see some kind of structure on my UMAP.

The problem is that apparently my data is not as good as I thought it was, I have a very low UMI per cell count and very low captured genes (we loaded 83k cells and only recovered 4k). Also, because the datasets are from very similar cells, the variability in the genes is very low and when I integrate it (with Harmony or Seurat) I get a UMAP that ends up being very diffuse with everything mixed.

I thought about driving my integration with the genes that are differentially expressed between the two datasets instead of the variable features. So I ran it like this:

merged_obj <- merge(
  x = saliva_83k,
  y = sg_40k_2,
  add.cell.ids = c("saliva", "sg")
)
merged_obj <- JoinLayers(merged_obj)

# Find differentially expressed genes between saliva and salivary_gland
de_genes <- FindMarkers(
  merged_obj,
  ident.1 = "saliva",
  ident.2 = "salivary_gland",
  min.pct = 0.1,           
  logfc.threshold = 0.25,  
  test.use = "wilcox",
  verbose = TRUE
)

# Filter for significant genes
sig_de_genes <- rownames(de_genes[de_genes$p_val_adj < 0.05, ])

# =============================================
# SCALE DATA USING ONLY DE GENES
# =============================================
merged_obj <- ScaleData(
  merged_obj,
  features = sig_de_genes,
  verbose = FALSE
)

# =============================================
# RUN PCA USING ONLY DE GENES AS FEATURES
# =============================================
merged_obj <- RunPCA(
  merged_obj,
  features = sig_de_genes,
  npcs = 30,
  verbose = FALSE
)

# Check elbow plot
ElbowPlot(merged_obj, ndims = 30)

# =============================================
# RUN HARMONY
# =============================================
integrated_obj <- RunHarmony(
  merged_obj,
  group.by.vars = "dataset",
  dims.use = 1:10,  
  theta = 0.5,
  lambda = 1,
  sigma = 0.1,
  max.iter.harmony = 10,
  verbose = FALSE
)

If I do this I get a beautiful UMAP. But I know that because I am using DEG between them to drive the integration, I can't really do any broad affirmations about anything I might find down the road. Everything would have to be viewed in a light of 'saliva vs gland'. At least that is my understanding of this.

Can anyone tell me if what I did makes sense and if I can actually proceed with this analysis or if this whole thing is just delusional of me?

seurat single cell integration
• 1,573 views • link
updated by antonioggsousa antonioggsousa • written by tais.bv • tais.bv

Hi,

Actually there is a software that does integration based on DEG genes, but not in the way you've described. The R package is called SCIBER (see the paper at https://academic.oup.com/bioinformatics/article/39/1/btac819/6957084?searchresult=1). First, it finds the DEGs in each cluster, within each dataset. Then, it matches the clusters between datasets that share the DEGs.

Integration assumes that the cells are shared across datasets. Therefore it attempts and "forces" the alignment between datasets. Thus, it is normal that after integration the cells overlap in the UMAP space. If you don't expect them to align, i.e. if you know that the cells across datasets are not the same, you shouldn't perform integration. In that case, you can merge the samples and analyze them together, which will probably give you a projection where the cells coming from the two datasets will not overlap in the UMAP space.

Best,

António

• 0 views • link
written by antonioggsousa antonioggsousa

1 answer

Use highly variable genes for the integration and try other methods CCA, RPCA to compare with Harmony. Choose the one that reflects the biological differences.

• 0 views • link
written by Arup Ghosh Arup Ghosh

I did that. I tried Harmony, CCA and RPCA both v4 and v5. Everything just gives me a difuse object with no structure. Because the datasets are too similar, the HVG doesn't really do anything for me. That's why I wanted to try with DEG, to force some kind of separation

• 0 views • link
written by tais.bv • tais.bv

What do you mean by difuse object with no structure? No cell clusters or the two datasets don't show separate cell clusters, share the umap plot.

• 0 views • link
written by Arup Ghosh Arup Ghosh

I played around with Harmony a little bit more today and noticed that if I change theta to 1 I get the structure I'm looking for. But I'm afraid I'm doing something wrong because I've played around with this data so much, I don't trust anything anymore. This is what I'm running and the UMAPs I'm getting

# =============================================================================
# LOAD DATA
# =============================================================================

saliva <- readRDS("RDS files/S.O.spz_saliva_83k.rds")
saliva$dataset <- 'saliva'

gland <- readRDS("RDS files/S.O.spz_sg_40k_2.rds")
gland$dataset <- 'gland'

# ============================================================================
# HARMONY INTEGRATION 
# ============================================================================

merged_obj <- merge(
  x = saliva, 
  y = gland,
  add.cell.ids = c("saliva", "gland")
)
merged_obj <- JoinLayers(merged_obj)
merged_obj <- NormalizeData(merged_obj, verbose = FALSE)
merged_obj <- FindVariableFeatures(merged_obj, selection.method = "vst")
merged_obj <- ScaleData(merged_obj)
merged_obj <- RunPCA(merged_obj, verbose = FALSE)
harmony_theta1 <- RunHarmony(
  merged_obj,
  group.by.vars = 'dataset',
  dims.use = 1:15,  
  theta = 1,
  verbose = FALSE
)
harmony_theta1 <- RunUMAP(
  harmony_theta1,
  reduction = "harmony",
  dims = 1:7,
  n.components = 2,
  n.neighbors = 30,
  min.dist = 0.5,
  verbose = FALSE
)

harmony_theta2 <- RunHarmony(
  merged_obj,
  group.by.vars = 'dataset',
  dims.use = 1:15,  
  theta = 2,
  verbose = FALSE
)
harmony_theta2 <- RunUMAP(
  harmony_theta2,
  reduction = "harmony",
  dims = 1:7,
  n.components = 2,
  n.neighbors = 30,
  min.dist = 0.5,
  verbose = FALSE
)

# ============================================================================
# SEURAT V5
# ============================================================================

merged_obj <- merge(
  x = saliva, 
  y = gland,
  add.cell.ids = c("saliva", "gland")
)

merged_obj <- JoinLayers(merged_obj)
merged_obj[["RNA"]] <- split(merged_obj[["RNA"]], f = merged_obj$orig.ident)
merged_obj <- NormalizeData(merged_obj, verbose = FALSE)
merged_obj <- FindVariableFeatures(merged_obj, selection.method = "vst")
merged_obj <- ScaleData(merged_obj)
merged_obj <- RunPCA(merged_obj, verbose = FALSE)

seurat_v5 <- IntegrateLayers(
  object = merged_obj, method = CCAIntegration,
  orig.reduction = "pca", new.reduction = "integrated.cca",
  verbose = FALSE
)

seurat_v5 <- RunUMAP(
  seurat_v5,
  reduction = "integrated.cca",
  dims = 1:7,
  n.components = 2,
  n.neighbors = 30,
  min.dist = 0.5,
  reduction.name = "umap_cca",
  verbose = FALSE
)

# ============================================================================
# SEURAT V4
# ============================================================================

combined <- list(saliva, gland)
for (i in 1:length(combined)) {
  combined[[i]] <- FindVariableFeatures(combined[[i]], selection.method = "vst", nfeatures = nrow(combined[[i]])*.3, verbose = FALSE)
}
combined_anchors <- FindIntegrationAnchors(
  object.list = combined, 
  dims = 1:30
)
combined <- IntegrateData(
  anchorset = combined_anchors, 
  dims = 1:30, 
  new.assay.name = "CCA"
)
combined <- ScaleData(
  combined,
  verbose = FALSE, 
  assay = "CCA"
)
combined <- RunPCA(
  combined, 
  npcs = 30, 
  verbose = FALSE, 
  assay = "CCA",
  reduction.name = "PCA_on_CCA"
)
combined <- RunUMAP(
  combined, 
  reduction = "PCA_on_CCA", 
  dims = 1:7, 
  reduction.name = "UMAP_on_CCA", 
  n.neighbors = 30, 
  min.dist = 0.5,
  verbose = FALSE
)

# =============================================================================
# PLOT UMAPS
# =============================================================================

wrap_plots(
  DimPlot(harmony_theta1, 
          reduction = "umap", 
          group.by = "orig.ident") +
    NoAxes() +
    ggtitle("Harmony Theta 1"),

  DimPlot(harmony_theta2, 
          reduction = "umap", 
          group.by = "orig.ident") +
    NoAxes() +
    ggtitle("Harmony Theta 2"),

  DimPlot(seurat_v5, 
          reduction = "umap_cca", 
          group.by = "orig.ident") +
    NoAxes() +
    ggtitle("Seurat v5"),

  DimPlot(combined, 
          reduction = "UMAP_on_CCA", 
          group.by = "orig.ident") +
    NoAxes() +
    ggtitle("Seurat v4"),

  ncol = 2
) + plot_layout(guides = "collect")

enter image description here

• 0 views • link
written by tais.bv • tais.bv

The default theta suggested by Harmony is 2, and both sets look well integrated if we assume the source of variation is from the batch. If you believe the variation in expression is only from biological differences, stick with the lower theta values, either 0 or 1. A higher theta value will make the clusters more homogeneous (cells from multiple batches in a cluster).

Ref: https://davetang.github.io/muse/harmony.html

• 0 views • link
written by Arup Ghosh Arup Ghosh

I was indeed expecting to see saliva_83k cluster away from sg_40K_2 even if they shared some cells, so theta = 1 did it for me. If the coding is sound then I guess I'll proceed with that.

Thank you for you help and for the reference!

• 0 views • link
written by tais.bv • tais.bv

Log in to answer this question.

More posts like this