DESeq2 design and stats question
Entering edit mode
rbronste ▴ 60
Last seen 4.0 years ago

So will just list current design first and then get into the question:

rangedCounts <- dba.peakset(Adult_BN_count, bRetrieve=TRUE)

nrows <- 1054182 
ncols <- 12
counts <- matrix(runif(nrows * ncols, 1, 1e4), nrows)

sampleNames<-c("MBV1",    "MBV2",    "MBV3",    "FBV1",    "FBV2",    "FBV3", "MBE7", "MBE8", "MBE9", "FBE1", "FBE2", "FBE3")
sampleSex<-c("male", "male", "male", "female", "female", "female", "male", "male", "male", "female", "female", "female")
sampleTreatment<-c("vehicle", "vehicle", "vehicle", "vehicle", "vehicle", "vehicle", "BB", "BB", "BB", "BB", "BB", "BB")
colData<-data.frame(sampleName=sampleNames, sex=sampleSex, treatment=sampleTreatment)

# Retrieve counts from dba.count without the site interval ranges

counts <- as.matrix(mcols(rangedCounts))

# Construct a SummarizedExperiment

se<-SummarizedExperiment(assays=list(counts=counts),rowRanges=rowRanges, colData=colData)

# DESeq2 analysis

table(sampleSex, sampleTreatment)

dds<-DESeqDataSet(se, design= ~sex + treatment)

dds$group <- factor(paste0(dds$sex, dds$treatment))

design(dds) <- ~0 + group

dds <- dds[ rowSums(counts(dds)) > 1, ]

dds <- DESeq(dds, betaPrior = FALSE)

So my question is about two types of results I am aiming for, group comparisons and pairwise comparisons, results method listed below:

#pairwise comparison

BN_group_male<-results(dds, format = c("GRanges"), contrast=c("group", "maleBB", "malevehicle"))

#group comparison

resFemVeh_BN<-results(dds, lfcThreshold = 1, altHypothesis = "greater", format = c("GRanges"), contrast=c(-1/3, 1, -1/3, -1/3))

My overall goal here is to see:

1. Example: What are the male vs female diff irrespective of treatment etc.

2. Example: What are differential peaks specific to only the Male BB condition and no other.

I am wondering if this is the proper way to get the results for these particular comparisons given the above design? 

deseq2 differential binding analysis deseq atac-seq • 1.3k views
Entering edit mode
Last seen 11 hours ago
United States

The first is straightforward: a comparison of the BB vs vehicle in males. The second is a comparison of the second group against an average of the others. Use plotCounts to get a sense of the top genes that the results tables are pulling out.

Entering edit mode


My central question here is whether the design makes sense for the two ways in which I am retrieving results? 

Also, in addition, I was noticing that when I lfcThreshold the #1 result type (a comparison of the BB vs vehicle in males) at the same cutoff as #2 I get pretty much no diff peaks (hence no lfc thresholding in the above #1 results formula), and wondering why that might be? Thanks again.

Entering edit mode

"Male vs female irrespective of treatment" you could get by averaging male vs female in each group. So putting 1/2 and -1/2 in the contrast. You don't seem to have anything like that in your code.

One vs the rest you can get with your second contrast.

Getting no diff peaks with cutoff: do some exploration and see if you find the answer. I would look at e.g. an MA plot of the results table.

Entering edit mode

Ok great thanks will put in that additional comparison, but in terms of the actual contrast I guess it would be the following, for male vs female only:

contrast=c(1/2, -1/3, -1/2, -1/3)

I guess just a little confused in this format for how to specify the comparison, not as much in the vignette. 

Thanks again!

Entering edit mode

Put 1/2 for coefficients associated with M and -1/2 for the coefficients associated with F. This is a bit of a complex contrast because it isn’t represented by any coefficient in the design. There isn’t a M vs F “irrespective” when you are fitting a group for each combination of sex and treatment. You may want to discuss this in more detail with a statistical collaborator. I mostly provide software feedback here, but have limited time for further statistical instruction.

Entering edit mode

Thanks again!


Login before adding your answer.

Traffic: 871 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6