
Hi, Mike and All:

 

I am testing DESeq2 for multi-factor contrast setting for my own data with more complex meta data but currently use the simpler dataset from the user guide for testing purpose, and run into some issues that need your input and advice.  Here are the commands (only show some more relevant outputs):

 

library("DESeq2")

library("Biobase")
library("pasilla")
data("pasillaGenes")
countData <- counts(pasillaGenes)
colData <- pData(pasillaGenes)[,c("condition","type")]
colData<-data.frame(colData,paste(colData$condition,colData$type,sep="."))
colnames(colData)[3]<-"condition_type";
> dds <- DESeqDataSetFromMatrix(countData = countData,colData = colData, design = ~ condition_type)
> colData(dds)
DataFrame with 7 rows and 4 columns
             condition        type        condition_type sizeFactor
              <factor>    <factor>              <factor>  <numeric>
treated1fb     treated single-read   treated.single-read  1.5116926
treated2fb     treated  paired-end    treated.paired-end  0.7843521
treated3fb     treated  paired-end    treated.paired-end  0.8958321
untreated1fb untreated single-read untreated.single-read  1.0499961
untreated2fb untreated single-read untreated.single-read  1.6585559
untreated3fb untreated  paired-end  untreated.paired-end  0.7117763
untreated4fb untreated  paired-end  untreated.paired-end  0.783745

>dds <- DESeq(dds)

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting generalized linear model
> resultsNames(dds)
[1] "Intercept"
[2] "condition_type_treated.single.read_vs_treated.paired.end"
[3] "condition_type_untreated.paired.end_vs_treated.paired.end"
[4] "condition_type_untreated.single.read_vs_treated.paired.end"

> colData(dds)
DataFrame with 7 rows and 4 columns
             condition        type        condition_type sizeFactor
              <factor>    <factor>              <factor>  <numeric>
treated1fb     treated single-read   treated.single-read  1.5116926
treated2fb     treated  paired-end    treated.paired-end  0.7843521
treated3fb     treated  paired-end    treated.paired-end  0.8958321
untreated1fb untreated single-read untreated.single-read  1.0499961
untreated2fb untreated single-read untreated.single-read  1.6585559
untreated3fb untreated  paired-end  untreated.paired-end  0.7117763
untreated4fb untreated  paired-end  untreated.paired-end  0.7837458


Then I tried  a slight diiffermry setting of the design:

> dds <- DESeqDataSetFromMatrix(countData = countData,colData = colData, design = ~0+ condition_type)

> dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting generalized linear model
580 rows did not converge in beta, labelled in mcols(object)$betaConv. Use larger maxit argument with nbinomWaldTest
> resultsNames(dds)
[1] "condition_typetreated.paired.end"
[2] "condition_type_treated.single.read_vs_treated.paired.end"
[3] "condition_type_untreated.paired.end_vs_treated.paired.end"
[4] "condition_type_untreated.single.read_vs_treated.paired.end"
> colData(dds)
DataFrame with 7 rows and 4 columns
             condition        type        condition_type sizeFactor
              <factor>    <factor>              <factor>  <numeric>
treated1fb     treated single-read   treated.single-read  1.5116926
treated2fb     treated  paired-end    treated.paired-end  0.7843521
treated3fb     treated  paired-end    treated.paired-end  0.8958321
untreated1fb untreated single-read untreated.single-read  1.0499961
untreated2fb untreated single-read untreated.single-read  1.6585559
untreated3fb untreated  paired-end  untreated.paired-end  0.7117763
untreated4fb untreated  paired-end  untreated.paired-end  0.7837458


Then supposedly, I can use results(dds, "condition_type_treated.single.read_vs_treated.paired.end") to get DEGs for each contrast shown in resultsNames(dds).

 

here are my questions:

1. I used design = ~0+ condition_type instead of design = ~ condition_type in 2nd case, try to skip the intercept so that I can easily get all possible contrasts,  but seem not working the way I want.  

2. I tried to get all possible contrasts:  but besides the contrasts shown in resultsNames(dds) in both cases, the contrasts like  untreated.single.read vs treated.single.read, untreated.paired.end vs untreated.single.read not even exists in the resultsNames(dds). also I like the contrast generally like: treated (including both treated.single.read and treated.paired-end) vs untreated (including both untreated.single-read and untreated  paired-end). I know for this case, we can just to design = ~condition, but I wish to do this in the same roof of one single design model although I can do a separate design. In limma and edgeR, there is a function like: con.matrix<-makeContrasts() where one can set up any contrasts under the design at will. Is there anythign like that in DESeq2? I understand we can do  design(dds) <- formula(~ condition_type), but no contrast setting can be made at will. Anything in DESeq2 can get around that?

3. Also for simple contrast, I understand one can use relevel(colData(dds)$condition,"control") kind of command to set base level, but for multiple-factors contrasts as I am after, I almost need some kind of makeContrasts() mechanism to set up contrasts at will or have to do that individually one by one, which obvioulsy would be tedious and also these contrasts won't be in a single model roof. Anything can get around like that as well? if question 2 is addressed, this one shall be no problem.

 

Thanks in advance for your help! Appreciated very much!

 

Best

 

Ming 

ATRF/NCI-Frederick,

Maryland, USA.

 

    

 



 
 		 	   		  
	[[alternative HTML version deleted]]

