Integrating joint snATAC/RNA-seq experiments

[cristanchoa]

Hello,

First off, thanks for the great tools in Signac/Seurat. I have a conceptual question that I think combines a few different of the vignettes and tutorials, but want to see if I'm thinking about it correctly before proceeding.

I have 16 10x runs from the joint ATAC/RNA platform that are for 4 conditions (2 treatments and by each sex) that each have 4 biological replicates. I am interested in various combination comparisons between treatment and sex for differential expression and accessibility.

I'm trying to understand when along the analysis it is appropriate to use NormalizeData, ScalaData, and Integrate with Harmony and if there are any special considerations when using the data from the joint chemistries or should I think of them as separate assays for the purpose of normalizing the data across the conditions/replicates? I had seen this but it is only for the one samples (https://satijalab.org/signac/articles/pbmc_multiomic.html)

Trying to decide on a workflow. Here are a couple I came up with:

1. Merge the 4 biological replicates from each condition into one Seurat object, then normalize and scale data, and the integrate the 4 conditions? It seems like I'd still want to go back to the scaled data for the DE and DA analysis so this doesn't quite seem right.
2. Make one large Seurat object with all the data and normalize and scale the data? It seems like I might end up normalizing for some of the interesting biological data this way, but perhaps could get around that with vars.to.regress for the nCount for the RNA seq? And then for the ATAC seq integrate the peaks using something like one of these vignettes (https://satijalab.org/signac/articles/integration.html or https://satijalab.org/signac/articles/merging.html)?

Thanks in advance for any guidance.
Ana

[kmwinkley]

I'm interested in the same topic here as well. I'm in a somewhat similar situation, where we have multiple captures of the 10x multiome (snRNA/ATAC) for several different groups of data that we eventually want to integrate. My approach so far is below:

• Assuming the difference between the captures for each group are only technical variation, I simply "merge" the Seurat objects for each capture of a group. I normalize the RNA using SCTransform with a batch variable of "capture". I the ATAC is recounted using a combined peak set that is identical for all groups and captures, and then normalized with TF-IDF. This produces a single Seurat object for each group that has been normalized to reduce sources of technical variability.
• Then the Seurat objects for each group are joined using the "integration" workflow. The RNA and ATAC are integrated separately and then joined into a single Seurat object which will have 5 assays (ATAC, RNA, SCT, integratedRNA, and integratedATAC). The integratedRNA and integratedATAC assays will be used for WNN dimensional reduction and clustering.

The only issue I'm encountering so far is some errors on the IntegrateData step of the ATAC integration. I'd be more than open to input on this workflow and it's shortcomings.

[cristanchoa]

Thanks for the breakdown of what you have done! What kind of error are you getting at the IntegrateData?

[kmwinkley]

The traceback seems to show that I'm ultimately failing during a merge.ChromatinAssay step. I was using Signac 1.1.0.9009, but there was a bug fix in the most recent version that fixed issue #355 that I think might be related to the error I'm getting. I've updated to Signac 1.1.1 and I'm trying again currently.

[coopershawna]

@kmwinkley Did the new version of Signac fix it?

[timoast]

Hi, this is an interesting question, and we're still working on methods to handle this properly (integration of multimodal datasets). In the meantime, I would suggest finding integration anchors using the RNA modality (using whatever options you want: CCA, rPCA). You can then use the same set of anchors to correct the cell embeddings for both the RNA data (run IntegrateData(), followed by ScaleData and RunPCA on the integrated assay), and for the LSI (run IntegrateEmbeddings() with reductions = "lsi"). You can then use both integrated dimension reductions for multimodal clustering (run FindMultiModalNeighbors() with the integrated RNA data and the corrected LSI as inputs).

[timoast]

@kmwinkley reductions needs to be a dimension reduction object containing all the cells to be integrated, so you also need to first merge the objects and run LSI on the merged object, then supply that LSI to the reductions parameter in IntegrateEmbeddings()

[kmwinkley]

Thanks for the clarification. That makes sense. The object that I'm integrating contains ~190k cells, so I'd like to avoid creating too many extra copies of these cells if possible. Would running lsi on the merged seurat object that results from the RNA integration step, and then supplying that DimReduc object to the IntegrateEmbeddings function be a reasonable alternative? I just want to make sure I'm not overlooking some modifications that may happen to the ATAC modality of the RNA-integrated object that may alter the results.

[timoast]

Yes, it just needs to be a merged object. If you have run integration on the RNA assay, you should have an object with an integrated RNA assay, the original RNA assay, and the original ATAC assay. You can run LSI on the ATAC assay to create the LSI with all cells present

[coopershawna]

Thank you! Not only are your vignettes amazing, REALLY appreciate willingness to answer questions too!

[kmwinkley]

Great, thanks for the help!

[JonathanNoonan]

@kmwinkley - Hello! Did you manage to get this to work? I'm trying it now and getting the same issue...but I believe I have ensured there is an LSI reduction for the combined objects. LSI is present for all cells, but it seems either i'm doing something wrong or its not happy with what I've done.

Would you mind sharing the script that worked for you if you were successful?

Cheers!

[hchintalapudi]

Hello all,
This is a very interesting discussion!
I have a similar use case where I have two multiome datasets (ATAC + RNA).
I followed this method and used the RNA integration anchors for integrating ATAC data but I got poor results or my integrated ATAC clusters.
How can I use a common/intersecting peak set in this scenario? I used the merge vignette's logic and that improved my ATAC integration marginally.
In my case, I start from two multiome objects and integrate them with each assay, but how can I contain the peaks to only the common peaks?

I did it separately by creating separate RNA integrated object and sepratae ATAC object (containg common peaks) but that had its own set of problems and is against the philosophy of this workflow.

If anyone could give some tips/advice on this, it would be appreciated.

Thanks.

[cristanchoa]

Do you mean between datasets it seems like the ATAC isn't overlapping well (ie not detecting the same cell types well?) In that case you might have perform an integration step with Harmony (https://satijalab.org/signac/0.2/articles/integration.html).

[hchintalapudi]

Thanks for the reply, @cristanchoa.
I have two datasets which belong to the same tumor/cell line but one is transplanted and one is cultured, sequenced separately and hence have a batch effect.
Ideally, they should integrate well and the RNA part integrated very well.
When it comes to ATAC, the ideal way to integrate is to choose or quantify the common peakset in both the samples, so I followed that.
Although I did not get overlapping clusters like RNA, there was a slight improvement.
So, I want to try that in the method described above. It looks easy but I'm missing some simple logic.
I followed the approach described above, used RNA integration anchors for
zmel_WNN.pdf
ATAC integration and then did WNN graph, see the plot attached
zmel_integrtaed_commonPeakset.pdf
.

[cristanchoa]

Running the above didn’t quite work for me for the exact same reason that there was a significant effect on my case of two treatments in developing brain (I had 16 samples so it was clearly not just a batch effect but a true treatment effect). I had to integrate the data with the runHarmony wrapper in Signac. It was in the vignette from Signac 1.2 but on my phone I don’t see it in the updated version. It worked. I created a “harmony” reduction. I then used that along with the rna pca into the parameters of wnn pipeline and it all worked quite well to create on integrated object I could use for downstream analysis.

Certainly open to feedback if this wasn’t quite right but the clusters of cell types I got in the brain were exactly as expected so hopefully pretty close.

[hchintalapudi]

Aah! Interesting.
I have not used Harmony with ATAC data.
It looks like there's no runHarmony function with Signac's latest version 1.4.0 I'm using.
Do you mind sharing that chunk of code you used for Harmony integration?
I'm looking to see if it is still doable with Harmony package rather than Signac's version of Harmony.
Right now, I only see this old vignette

Thanks.

[cristanchoa]

I think below is what you are looking for. Hope it helps. Good luck!
Quick edit: this is after you have made it all the data into one object with combined peaks using the initial vignettes.

`library(Signac)
library(Seurat)
library(SeuratWrappers)
library(ggplot2)
library(patchwork)
library(GenomeInfoDb)
library(EnsDb.Mmusculus.v79)
library(tidyverse)
library(dplyr)
library(grid)
library(readr)
library(harmony)

DefaultAssay(all_runs)<-"peaks"
all_runs<- RunTFIDF(all_runs)

all_runs <- FindTopFeatures(all_runs, min.cutoff = 20)
all_runs <- RunSVD(all_runs)

all_runs <- RunUMAP(all_runs, dims = 1:50, reduction = 'lsi', reduction.name = "umap.ATAC")

all_runs <- RunHarmony(object = all_runs, group.by.vars = 'condition',
reduction = 'lsi',
assay.use = 'peaks',
project.dim = FALSE)

all_runs<- RunUMAP(all_runs, dims = 1:50, reduction = 'harmony', reduction.name = "umap.ATAC_harm")

DefaultAssay(all_runs)<-"RNA"

all_runs<-SCTransform(all_runs, verbose = TRUE)
all_runs<- RunPCA(all_runs, verbose = TRUE)
all_runs <- RunUMAP(all_runs, dims = 1:50, verbose = TRUE)

all_runs <- FindMultiModalNeighbors(all_runs, reduction.list = list("pca", "harmony"), dims.list = list(1:50, 1:50))
all_runs <- RunUMAP(all_runs, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")

`

[hchintalapudi]

Thanks, @cristanchoa !
I tried this and this is what I noticed.
Running Harmony only improved things marginally when I used a separate merged ATAC object with a common peakset.
It was pretty much the same when I used the integration approach described above.
I think my data is such but I'm back to the same problem again:
How do I create a merged multiome object where I have merged RNA data (which I will integrate later) and only common set of peaks in the ATAC data? In the merging vignette, the fragment objects and the seurat objects are constructed from scratch but I'm not able to add my RNA counts since the number of cells differ in the samples.
Probably I'm missing some trick with subsetting the objects

[micosacak]

Hi @cristanchoa ,
It is not clear for me, how to get all_runs. Could you please share the part of code to get all_runs?

Basically, I have 2 multiome datasets and trying to integrate them.
thanks,
ilyas.

[cristanchoa]

I basically followed the vignette for the ATAC (https://stuartlab.org/signac/articles/merging) and it was the older version of the SCT transform integration with the RNA (there is now an updated one with the seurat V5 object that I think is cleaner). Then, as suggested by Tim above used the weighted nearest neighbors (https://satijalab.org/seurat/articles/weighted_nearest_neighbor_analysis) analysis to bring it all together.