This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Can't seem to integrate similar datasets

Hi everyone!

I've been having a hard time integrating some scRNAseq datasets and wanted to know if anyone can shed a ligth on my problem.

I have 3 datasets from the same experiment that have been previously published by another group and 2 datasets from my lab. They are all from the same condition/tissue, with the only difference being that the published dataset was prepared using v3 Chromium Single Cell and sequenced on a HiSeq X Ten and our datasets were prepared using HIVE and sequenced on a NovaSeqXSeries.

I got the raw fastq files from the published datasets and processed it using our pipeline. The paper says they loaded around 5k - 8k cells and our pipeline recovered ~4k cells. For our dataset, we loaded 83k/40k cells and recovered around the same 4k. Maybe because of that, the published data had a much higher UMI count avg than our data.

For the integration, I tried everything! I did Harmony, Seurat CCA v4 and v5 with NormalizeData and SCTransform. I even tried downsampling the published data to have a similar UMI avg.

No matter what I did, they won't integrate. And I am sure of that because every time I get two different clusters that were supposed to be the same stage. One coming from the published data and the other from our data. Either that or they become a diffuse blob with zero structure.

Can anyone give me a tip on how to integrate these guys? Anything at all is appreciated!

seurat integration
• 1,319 views • link
written by tais.bv • tais.bv

What are the genes driving each of your clusters ? Could they be related to your experimental design ? PolyA capture differences ?

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

That having said, can you show plots and relevant code please? Sounds like a comparably trivial integration, so there is probably a strong confounder here that needs to be identified. Could be a set of genes not captured at all in one of the two experimental protocols, but well captured in the other. Might make sense to subset to the "common universe" first, so genes detected in at least x % of cells in BOTH datasets, and removing genes nested by protocol.

• 0 views • link
written by ATpoint ATpoint

This is a great idea, thank you!

I created a list of genes that were expressed in at least 5% of cells and in at least 4 of my 5 datasets and it was only 312 genes long (this very low number is probably the reason to all of my problems!). Then I tried Harmony integration but doing scaling and PCA using features = common_genes. Unfortunately that didn't work, my data was still clustering away from the published data. This is how I ran that:

seurat_list <- list(m38, m39, m40, saliva, gland)
joint_integration <- merge(
  x = seurat_list[[1]], 
  y = seurat_list[2:length(seurat_list)],
  add.cell.ids = c("m38", "m39", "m40", "saliva", "gland")
)

# Normalize all genes
joint_integration <- NormalizeData(joint_integration, verbose = FALSE)

# Scale only the common genes
joint_integration <- ScaleData(
  joint_integration,
  features = common_genes
)

# Run PCA on common genes only
joint_integration <- RunPCA(
  joint_integration,
  features = common_genes,
  verbose = FALSE
)

joint_integration <- RunHarmony(
  joint_integration,
  group.by.vars = 'orig.ident',
  dims.use = 1:20,  
  theta = 2,
  verbose = FALSE
)

joint_integration <- RunUMAP(
  joint_integration,
  reduction = "harmony",
  dims = 1:10,
  n.neighbors = 30,
  min.dist = 0.5,
  verbose = FALSE,

)

DimPlot(joint_integration, group.by = "orig.ident")

So I tried doing Seurat's CCA v4 integration. This is what I ran:

combined <- list(m38, m39, m40, 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:15)

combined <- IntegrateData(anchorset = combined_anchors, dims = 1:15, 
                          new.assay.name = "CCA")
combined <- ScaleData(combined, features = common_genes, verbose = FALSE, assay = "CCA")
combined <- RunPCA(combined, features = common_genes, npcs = 30, verbose = FALSE, assay = "CCA",
                   reduction.name = "PCA_on_CCA")
combined <- RunUMAP(combined, reduction = "PCA_on_CCA", dims = 1:15, 
                    reduction.name = "UMAP_on_CCA",  
                    n.neighbors = 20, min.dist = 0.5)

DimPlot(combined, group.by = 'orig.ident')

This actually gave me something that was mixed. Not exaclty what I was expecting to see, but I can work with that! I just wanted to make sure what I did is correct and defensible in a paper and won't affect my downstream analysis...

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

By the way, this is the UMAP for the Harmony integration enter image description here And this is for the Seurat v4 enter image description here

I was expecting the blue dataset to cluster separately from the pink, but oh well!

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

0 answers

No answers yet.

Log in to answer this question.

More posts like this