Hi,
I am new to DESeq2 and I am finding some difficulties when building my
contrasts.
In my experimental design I have 2 subjects (defined by depth from which I sampled the marine microbial commmunity, Depth), to which I performed 3 different treatments (Treatment) and I took samples 3 times (Time). I
have therefore the following conditions:
-Treatment: 3 leves (C, DOM, DOMNP)
-Time: 3 levels (00d, 02d, 25d)
-Depth: 2 levels (25m, 125m)
I start with a phyloseq object (PHYLUM) that I convert to DESeq (Phylum)
easily, and to make it easy I start with a simple contrast in one
factor:
> Phylum=phyloseq_to_deseq2(PHYLUM,~Treatment)
> Phylum
class: DESeqDataSet
dim: 57 30
exptData(0):
assays(1): counts
rownames(57): OTU1 OTU2 ... OTU56 OTU57
rowData metadata column names(0):
colnames(30): T2_SF5C T25_SF1B ... T2_DP1A T2_DP1C
colData names(3): Treatment Time Depth
#I first Set base level to treatment Control and time 0
> colData(Phylum)$Treatment<-relevel(colData(Phylum)$Treatment,"C")
> colData(Phylum)$Time<-relevel(colData(Phylum)$Time,"00d")
#I then apply deseq:
> Phylum<-estimateSizeFactors(Phylum)
> Phylum<-estimateDispersions(Phylum,fitType="local")
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
> plotDispEsts(Phylum)
> Phylum<-nbinomWaldTest(Phylum)
> plotMA(Phylum)
> meanSdPlot(log2(counts(Phylum,normalized=TRUE)))
#I then call resultsNames()
> resultsNames(Phylum)
[1] "Intercept" "TreatmentC" "TreatmentDOM" "TreatmentDOMNP"
Here is where I am comfused. Attending to the manual I should get:
> resultsNames(Phylum)
[1] "Intercept" "Treatment_DOM_vs_C" "Treatment_DOMNP_vs_C"
Since the resultsNames() has not provided contrasts but just the
different levels of the factor Treatment I find this error when I cal
results():
> res<-results(Phylum, "TreatmentDOMNP")
Error in results(Phylum, "TreatmentDOMNP") :
'contrast', as a character vector of length 3, should have the form:
contrast = c('factorName','numeratorLevel','denominatorLevel'),
see the manual page of ?results for more information
> head(res,2)
Error in head(res, 2) :
error in evaluating the argument 'x' in selecting a method for
function 'head': Error: object 'res' not found
I have 2 options to do the contrast (below) but I would like to know
what I am doing wrong so I do not get the contrast in the resultsNames()
from the begining.
OPTION A
> res<-results(Phylum)
> res<-res[order(res$padj),]
> head(res)
log2 fold change (MAP): Treatment DOMNP vs C
Wald test p-value: Treatment DOMNP vs C
DataFrame with 6 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue
padj
<numeric> <numeric> <numeric> <numeric> <numeric>
<numeric>
OTU39 11.53402 1.0638491 0.2560363 4.155071 3.251859e-05
0.001268225
OTU14 47.93680 -0.4802443 0.1285006 -3.737291 1.860137e-04
0.001846475
OTU21 635.45233 1.3622669 0.3590578 3.794005 1.482367e-04
0.001846475
OTU48 1097.05931 1.3071717 0.3501877 3.732775 1.893821e-04
0.001846475
OTU10 2979.93247 1.4262973 0.4060197 3.512877 4.432830e-04
0.003457607
OTU2 53.32792 0.9679247 0.3266231 2.963430 3.042308e-03
0.017682230
OPTION B
> DOMvC<-results(Phylum,contrast=c("Treatment", "DOMNP",
"C"),independentFiltering=FALSE,cooksCutoff=FALSE)
> DOMvC<-DOMvC[order(DOMvC$padj,na.last=NA),]
> DOMvC
log2 fold change (MAP): Treatment DOMNP vs C
Wald test p-value: Treatment DOMNP vs C
DataFrame with 57 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue
padj
<numeric> <numeric> <numeric> <numeric> <numeric>
<numeric>
OTU39 11.53402 1.0638491 0.2560363 4.155071 3.251859e-05
0.001853560
OTU14 47.93680 -0.4802443 0.1285006 -3.737291 1.860137e-04
0.002698695
OTU21 635.45233 1.3622669 0.3590578 3.794005 1.482367e-04
0.002698695
OTU48 1097.05931 1.3071717 0.3501877 3.732775 1.893821e-04
0.002698695
OTU10 2979.93247 1.4262973 0.4060197 3.512877 4.432830e-04
0.005053426
... ... ... ... ... ...
...
OTU9 2.583912 -0.05975886 0.3911841 -0.15276402 0.8785844
0.9361963
OTU11 96.893172 0.01630965 0.1146960 0.14219896 0.8869229
0.9361963
OTU1 5.168939 -0.03418351 0.3801815 -0.08991366 0.9283558
0.9521476
OTU25 5.163909 0.01949842 0.2752228 0.07084593 0.9435204
0.9521476
OTU29 9.956546 -0.01226835 0.2044383 -0.06001004 0.9521476
0.9521476
One quick note:
Unfortunately I do not have replicates for everything and that may be
the problem:
For each depth, at Time 00d I just have 1 replicate fo the control
(C_00d). At time 2 days I have 3 replicates for the control (C_02d), 2
replicates for the DOM treatment (DOM_02d) and 3 replicates for the
DOMNP treatment (DOMNP_02). At time 25 days I have 3 replicates for the
control (C_25d), 2 replicates for the DOM treatment (DOM_02d) and 2
replicates for the DOMNP treatment.
THANKS SO SO MUCH!!!
