Questions about identifying differentially accessible peaks in pairwise comparison

[linlejiang]

Great work on the Libra package! But recently, I encountered an issue while using run_de to analyze scATAC-seq data, and I was hoping you could offer some insight. Here’s a brief overview of my workflow:

Objective: I’m using ArchR alongside Libra to identify disease-specific differentially accessible peaks in basal cells.

Subsetting: I first subset the cell type to 'basal' and then refined the ArchR object to include only the disease group and control group for comparison.

Approaches Tried: I’ve experimented with various single-cell and pseudobulk methods, but the resulting peaks show very small effect sizes (avg_logFC between [-0.2, 0.2]). Additionally, the .exp and .pct values for both groups are quite low.

Discrepancy: I noticed significant differences in avg_logFC values between the Wilcoxon test results from Libra and those from ArchR, which has left me uncertain about potential flaws in my analysis.

Below is the relevant code snippet for reference:

idx_copd_condition <- BiocGenerics::which(basal_copd_condition_pseudobulk_peak_matrix$copd_condition %in% c("GOLD_I", "Control"))
cellsSample <- basal_copd_condition_pseudobulk_peak_matrix$cellNames[idx_copd_condition]
basal_peakmatrix_GOLD_I <- basal_copd_condition_pseudobulk_peak_matrix[cellsSample, ]

meta <- basal_peakmatrix_GOLD_I@cellColData[,c('Sample','CellType_epi','copd_condition')]
meta$Sample <- as.vector(meta$Sample)
meta$CellType_epi <- as.vector(meta$CellType_epi)
meta$copd_condition <- relevel(factor(meta$copd_condition), ref = "GOLD_I")

basal_peakmatrix <- getMatrixFromProject(
ArchRProj = basal_peakmatrix_GOLD_I,
useMatrix = "PeakMatrix",
binarize = FALSE,
verbose = TRUE
)

peaks <- rowRanges(basal_peakmatrix)
peak_labels <- paste0(seqnames(peaks), ":", start(peaks), "-", end(peaks))
rownames(basal_peakmatrix) <- peak_labels

DE <- run_de(assay(basal_peakmatrix),
meta = meta,
replicate_col = 'Sample',
cell_type_col = "CellType_epi",
label_col = "copd_condition",
de_family = 'pseudobulk',
de_method = 'edgeR',
de_type = 'LRT',
input_type = 'scATAC',
normalization = 'log_tp10k',
binarization = FALSE,
n_threads = 30
)

Here is the input basal_peakmatrix:

class: RangedSummarizedExperiment
dim: 186457 2755
metadata(0):
assays(1): PeakMatrix
rownames(186457): chr1:819660-820160 chr1:826648-827148 ...
chrX:155881027-155881527 chrX:155887434-155887934
rowData names(1): idx
colnames(2755): LC_18_800014#CTTGTCGTCAATTGGC-1
LC_18_800014#ACAAAGATCCAATCCC-1 ...
LC_18_800024#CCAGATATCCATGACA-1
LC_18_800024#GACTAACGTGCTGTCG-1
colData names(40): BlacklistRatio DoubletEnrichment ... ReadsInPeaks
FRIP

And a brief summary about the distributions of the output generated after applying the pseudobulk methods using Libra (results using Wilcoxon method are not shown in the image, but the distributions are similar to the pseudobulk methods) and the default Wilcoxon method from ArchR:

[Pirhhot]

Hello,

I am in a very similar case as you as I am also trying to use ArchR for the peak calling and Libra for paired Differential accessibility analysis on an experimental condition in a cell specific manner.

I was wondering how did you proceed for the peak calling step as I guess it could have a big impact on your analysis.
Specifically, my problem is that I am not sure of the best approach to define the peaks in this situation.

The two options that I have in mind is to either consider that the peaks specific to each cell type are going to be retrieved either way and that we don't need to specify it to ArchR :

proj_archr <- addGroupCoverages(
  ArchRProj = proj_archr,
  groupBy = "copd_condition", # in your case
  useLabels = TRUE,
  sampleLabels = "Sample",
  minReplicates = number_of_samples # for maximal use of biological reproducibility 
)

The other possibility is to consider that we should subset the ArchR project by condition AND cell type so that it could use the enhanced sensibility for cell type specific peaks :

proj_archr <- addGroupCoverages(
  ArchRProj = proj_archr,
  groupBy = "copd_condition_by_cell.type", # create a combined cellColData column 
  useLabels = TRUE,
  sampleLabels = "Sample",
  minReplicates = number_of_samples
)

I have not made systematic comparison but I would be interested in your opinion.

[linlejiang]

Hello,

I am in a very similar case as you as I am also trying to use ArchR for the peak calling and Libra for paired Differential accessibility analysis on an experimental condition in a cell specific manner.

I was wondering how did you proceed for the peak calling step as I guess it could have a big impact on your analysis. Specifically, my problem is that I am not sure of the best approach to define the peaks in this situation.

The two options that I have in mind is to either consider that the peaks specific to each cell type are going to be retrieved either way and that we don't need to specify it to ArchR :

proj_archr <- addGroupCoverages(
  ArchRProj = proj_archr,
  groupBy = "copd_condition", # in your case
  useLabels = TRUE,
  sampleLabels = "Sample",
  minReplicates = number_of_samples # for maximal use of biological reproducibility 
)

The other possibility is to consider that we should subset the ArchR project by condition AND cell type so that it could use the enhanced sensibility for cell type specific peaks :

proj_archr <- addGroupCoverages(
  ArchRProj = proj_archr,
  groupBy = "copd_condition_by_cell.type", # create a combined cellColData column 
  useLabels = TRUE,
  sampleLabels = "Sample",
  minReplicates = number_of_samples
)

I have not made systematic comparison but I would be interested in your opinion.

The way I approached this problem, as indicated in my description, is similar to your second method. Because I am interested in the differentially accessible peaks specific to certain cell types, this approach seemed more reasonable, and that is why I subset the cell type of interest in advance. But the problem I have is the tiny avg_logFC, exp, and pct values computed by Libra. I am not sure if this pattern was due to the data or the procedures I took.

[Pirhhot]

How many samples do you have ? One of the ideas that I have for your case is that maybe the difference is mainly driven by a subpart of the samples and the single cell approach could enhance this effect ?
Maybe try plotting an heatmap in order to look into inter sample variability ?