I have a differential expression design with multiple conditions, say A, B, C, D.

Using edgeR I'd like to find genes that are upregulated in condition A vs each of the others conditions, i.e. genes upregulated in: A vs B and in A vs C and in A vs D.

My current strategy is to do all the pairwise comparisons and than intersect the output to find genes with e.g. LogFC > 0 and FDR < 0.05 in each comparisons. Like (partially pseudo-code):

design <- model.matrix(~0 + group)
contrasts <- makeContrasts(
    A_vs_B= A - B,
    A_vs_C= A - C,
    A_vs_D= A - D,
    levels= design
)
fit <- glmFit(y, design)

foreach contrast in contrasts:
    lfc <- glmTreat(fit, contrast= contrast)
    detable <- append(detable, topTags(lfc, n= Inf))

# Now query detable

Is there a better strategy? I think a contrast as A_vs_all= A - (B + C + D)/3 won't work because it compares A vs the average of B, C, D.

(For the record, I don't have genes, I have ATACseq regions and I'm trying to answer the question which ATAC sites are specific to one particular condition?)

Kevin has attempted a stricter interpretation of your statement about finding regions specific to one condition, but I will assume that you are simply trying to find regions that are significantly more open (or less open) in A than in any of the other conditions. That's the question posed by the title of your question and is a perfectly sensible way to proceeed. In that case, doing the three comparisons of A with B,C,D as you have is exactly what you need to do.

Collating the three results is quite easy, for example:

AvsB <- glmTreat(fit, contrast=contrast[,1])
AvsC <- glmTreat(fit, contrast=contrast[,2])
AvsD <- glmTreat(fit, contrast=contrast[,3])
Results <- cbind(decideTests(AvsB), decideTests(AvsC), decideTests(AvsD))
RegionsUp <- which( rowSums(Results)== 3 )
RegionsDown <- which( rowSums(Results)== -3 )

Defining a reference set of peaks/regions

IMO there are four major ways to define the regions for differential accessibility analyses of ATAC-seq.

1. Focus on genes, i.e., preset putative promoter regions or gene bodies (e.g., PMID: 23375371).

2. Tile windows across the genome, possibly merging adjacent windows that appear concordant (e.g., csaw package).

3. Pool all the libraries together and use a peak-caller on the the pool. Do a differential accessibility anaysis for the resulting consensus peaks.

4. Call peaks on each library separately and combine the peaks using some ad hoc intersection/union method.

I use Methods 1-3, all of which control the FDR correctly and have no problem with differing sequencing depths. Method 4 is the most commonly used in the literature, but also the only one that I consider unsatisfactory.

Hey Dario, good to see you here.

which ATAC sites are specific to one particular condition?

Is this what a differential expression comparison (any comparison) is going to show? A differential expression comparison will tell you what is statistically significantly different between, e.g., A and the other groups, but not necessarily what is unique / specific to A.

In my mind, what you need to do here are the pairwise comparisons, take common peaks / regions from these, and also apply some further metric, such as Z-score or read cut-offs. That is, you need p-values but also some other threshold to define what is truly open [chromatin].

Kevin