optimal parameters in addGroupCoverages()

[trinevd]

Hi!
Thank you for developing this great tool, it already gave me lots of great biological insights into my datasets.

I am wondering if you could elaborate a bit on how to chose the most optimal parameters in the addGroupCoverages() command? I did read the explanations in the full manual and understand how the replicates are composed in different scenarios. However, I am in doubt of the exact number of pseudobulk replicates I should aim for. And since the downstream peak-based analysis is the core of the scATAC-seq analysis I decided to reach out to you.

Specifically, I am thinking about these 4 parameters:

minReplicates
maxReplicates
minCells
maxCells

For example I have a setup with 8 samples in total (2 conditions with 4 biological replicates). Cells from these 8 samples are clustered into 3 clusters and the cell count per sample and cluster distributes like this:

As you can see from the figure, C1 is a mix of the two conditions whereas C2 and C3 are more enriched for one of the two conditions (WD and chow). Also, when I look at the UMAP, WD-cells and Chow-cells are distributed in slight separate areas of C1 so there might be slight changes in the chromatin landscape of this C1-celltype across conditions.

So my questions would be:

1. Should I optimally aim to get 8 sample-aware pseudobulk replicates and chose minCells and maxCells according to this?
2. In C3, there are quite few cells from 2 samples. Should I, for this reason, set minReplicates = 6?
3. Do I optimally want pseudobulk replicates to be composed from as many (sample-aware) cells as possible and thus set maxCells=1370?

My command would in this case look like this:

MergedProj <- addGroupCoverages(
  ArchRProj = MergedProj,
  groupBy = "Clusters",
  useLabels = TRUE,
  minCells = 40,
  maxCells = 1370, 
  maxFragments = 25 * 10^6,
  minReplicates = 6,
  maxReplicates = 8,
  sampleRatio = 0.8, 
  kmerLength = 6,
  threads = getArchRThreads(),
  returnGroups = FALSE,
  parallelParam = NULL,
  force = TRUE, 
  verbose = TRUE,
  logFile = createLogFile("addGroupCoverages")
)

Any feedback or comments on this approach would be highly appreciated.
Best,
Trine

[rcorces]

What you are asking is extremely specific to your use case and I'm afraid I dont have more to add than what is in the manual. If you would like to ask a general question about one of the 4 parameters you mentioned, perhaps I could answer that.

In general, I think your problem is that you are grouping by cluster which does not suit your use case. You should define a new column in cellColData that looks like condition_cluser (for ex. Chow1_C1) and use this to create your pseudobulks. Otherwise, you will be grouping cells from your two different conditions together because they fall in the same cluster.

[trinevd]

Hi again,
Thanks for your answer. I think it was the grouping parameter I was confused about when I am working with 2 different conditions. It makes very much sense to separate it according to condition_cluster (I assume your mean Chow_C1 in your example above).

In case other people will find it useful I just summarize here what I ended up doing:
I am running the analysis using the commands below. With these parameters, I end up having 3-4 pesudobulk replicates for each condition in each cluster. I have set the minCells and maxCells such that each sample can contribute with a max number of cells to the sample-aware psudobulk replicates. Only in my C3, a pseudobulk replicate is generated from two samples since these have very few cells (Chow2 and Chow4).

#Add info on condition matched with cluster identity to the project
MergedProj <- addCellColData(ArchRProj = MergedProj, data = paste0(MergedProj@cellColData$Clusters,MergedProj@cellColData$Cond), name = "ClusterByCond", cells = getCellNames(MergedProj), force = TRUE)


#Generate pseudobulk replicates
MergedProj <- addGroupCoverages(
  ArchRProj = MergedProj,
  groupBy = "ClusterByCond",
  useLabels = TRUE,
  minCells = 40,
  maxCells = 1400,
  maxFragments = 25 * 10^6,
  minReplicates = 3,
  maxReplicates = 4,
  sampleRatio = 0.8,
  kmerLength = 6,
  threads = getArchRThreads(),
  returnGroups = FALSE,
  parallelParam = NULL,
  force = FALSE, 
  verbose = TRUE,
  logFile = createLogFile("addGroupCoverages")
)

Best,
Trine