Dear all,
I am trying edgeR to detect DE gene for my time course data. I have one experimental condition with 4 time point, But while using design matrix its give me a error.
My code is: (in short:)
> label
[1] WT WT WT WT WT WT WT WT WT WT WT WT WT WT WT WT WT
Levels: WT
> Time
[1] T1 T1 T1 T1 T1 T1 T2 T2 T2 T2 T3 T3 T3 T3 T4 T4 T4
Levels: T1 T2 T3 T4
# Before we fit GMLs, we need to define our design matrix
design <- model.matrix (~Time + label)
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
Any advice? Thanks in advance guys
•
written
by
bioinforupesh2009.au •
Since label serves no purpose, try design <- model.matrix(~Time)and let us know if that doesn't fix things.
•
written
by
Devon Ryan
Thank you very much Ryan....you are fabulous. its work :)))
•
written
by
bioinforupesh2009.au •
Dear Ryan:
I am wondering where I have Time 1 as I have 4 time in my series ??
> design
(Intercept) TimeT2 TimeT3 TimeT4
WT.T1.1 1 0 0 0
WT.T1.2 1 0 0 0
WT.T1.3 1 0 0 0
WT.T1.4 1 0 0 0
WT.T1.5 1 0 0 0
WT.T1.6 1 0 0 0
WT.T2.1 1 1 0 0
WT.T2.2 1 1 0 0
WT.T2.3 1 1 0 0
WT.T2.4 1 1 0 0
WT.T3.1 1 0 1 0
WT.T3.2 1 0 1 0
WT.T3.3 1 0 1 0
WT.T3.4 1 0 1 0
WT.T4.1 1 0 0 1
WT.T4.2 1 0 0 1
WT.T4.3 1 0 0 1
attr(,"assign")
[1] 0 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Time
[1] "contr.treatment"
You specified to keep the intercept, so it's not an explicit coefficient (everything is being compared to it). If you want to use contrasts instead, then use model.matrix(~0+Time).
Dear Devon,
Sorry for late comments for the same question: Advised model.matrix(~0+Time) may I use coef =2:4?
lrt <- glmLRT(fit,coef=2:4) is it ok?
If yes, then I got about 700 genes (miRNAs) with hugest logFC out of ~1200. is it right? I have doubt? Something is wrong!
#Estimate Dispersion
design <- model.matrix (~0 +Time)
rownames(design) <- colnames(y)
design
y <- estimateGLMCommonDisp(y, design, verbose =T)
y <- estimateGLMTrendedDisp(y,design)
y <- estimateGLMTagwiseDisp(y,design)
fit <- glmFit(y,design)
lrt <- glmLRT(fit,coef=2:4)
topTags(lrt)
Coefficient: TimeT2 TimeT3 TimeT4
logFC.TimeT2 logFC.TimeT3 logFC.TimeT4 logCPM LR PValue FDR
mmu-miR-5619-5p ** -20.45565** -20.22821 -19.54368 1.159783 43579141 0 0
mmu-miR-28c ** -20.45565 ** -19.69322 -19.92847 1.161204 43576050 0 0
mmu-miR-6395 **-20.45565** -19.93606 -20.08351 1.163118 43334784 0 0
.....
If you're going to use model.matrix(~0+Time) then you really need to use contrasts. Otherwise, you end up testing whether you have signal from an miRNA at T2, T3 and/or T4...which is probably not what you wanted to test. If you want to use glmLRT(fit,coef=2:4), then use model.matrix(~Time) instead. This would then ask the question, "Is there an effect of time on a given miRNA." This is likely what you actually want to know.
Thank Devon for you generous advice. I tried the same and got 0 significant miRNA of FDR < 0.05 while having about 51 significant raw pvalue, could you please tell me how it happened.
# noemalization
y <- DGEList(counts=wt,group=label)
#Estimate Dispersion
design <- model.matrix (~Time)
rownames(design) <- colnames(y)
design
y <- estimateGLMCommonDisp(y, design, verbose =T)
y <- estimateGLMTrendedDisp(y,design)
y <- estimateGLMTagwiseDisp(y,design)
fit <- glmFit(y,design)
lrt <- glmLRT(fit,coef=2:4)
topTags(lrt)
> sum(lrt$table$PValue < 0.05)
[1] 51
> sum(lrt$table$FDR < 0.05)
[1] 0
If there were only 51 genes with raw p-values <0.05, then it's entirely expected that you'd have no significant findings. It seems that time has no effect in your dataset.
Do you think, gene filtering step can affect to find the significant genes ??? May i try DESeq2 packages as he has gene filtering steps ??? Do you think its worth to try ??? what edgeR does not have it.
Sure, give it a try. In your case I doubt it'll make a difference, but it won't hurt to try.
•
written
by
Devon Ryan
Dear Dev,
Could you help me once again?
I tried but got error while using DESeqDataSeq
library(DESeq2)
x= read.table("counts_filt.txt", header=T)
rownames(x) <- x$miRNA
x <- x [, -1]
design=read.table("design_deseq.txt", header=T)
>design
strain time.m replicates id
WT.T1.1 wt 1 r1 wk
WT.T1.2 wt 1 r1 wk
WT.T1.3 wt 1 r1 wk
WT.T2.1 wt 12 r2 year
WT.T2.2 wt 12 r2 year
WT.T2.3 wt 12 r2 year
WT.T2.4 wt 12 r2 year
WT.T3.1 wt 21 r3 two_years
.......
.................
WT.T4.6 wt 36 r4 old
ddsTc <- DESeqDataSet(x, design= ~design$strain + design$time.m + design$strain:design$time.m)
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function 'assays' for signature '"data.frame"'
Please read help(DESeqDataSet), noting in particular the DESeqDataSetFromMatrix() function.
•
written
by
Devon Ryan
##countData; I called edgeR's DGEList bcz got error in counts function of DESEQ2
countData <- DGEList(counts=x,genes=x[,1],)
dds <- DESeqDataSetFromMatrix(countData= countData, colData = colData, design = ~ Time + label + Time:label)
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more level
However, I used only Time
dds <- DESeqDataSetFromMatrix(countData= countData, colData = colData, design = ~ Time)
ddsLRT= DESeq (dds, test="LRT", reduced = ~Time)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
**NOTE: fitType='parametric', but the dispersion trend was not well captured by the
function: y = a/x + b, and a local regression fit was automatically substituted.
specify fitType='local' or 'mean' to avoid this message next time.
final dispersion estimates
fitting model and testing
Error in nbinomLRT(object, full = full, reduced = reduced, betaPrior = betaPrior, :
less than one degree of freedom, perhaps full and reduced models are not in the correct order
Why it happened?
•
written
by
bioinforupesh2009.au •
I seriously doubt that using a DGElist as the input matrix will work properly. Also, you need to show what colData is.
•
written
by
Devon Ryan
wait I will update it soon, I got it. Thanks for advice.
•
written
by
bioinforupesh2009.au •
It appears you're trying to get through an analysis without reading any of the documentation for the software. When you encounter an error there is a bit of troubleshooting you should do on your part. For one thing read the help for the function that gave an error, here that means type ?DESeq into the R console and read the meaning of the full and reduced formula.
The help files tell you what kind of object each argument is asking for. By typing class(x) you can see if 'x' is the appropriate object to be giving for each argument.
•
written
by
Michael Love
Log in to answer this question.
As Devon say below label serves no purpose, when trying to generate the design.matrix you need more than one level, other wise R cannot create the model representation.