[BioC] interaction effects for DEXSeq
Mallon, Eamonn B. (Dr.)
ebm3 at leicester.ac.uk
Wed Jun 12 16:44:37 CEST 2013
Dear Simon,
Thanks I wouldn't have got that in a month of sundays.
Unfortunately no exons come up as significant. Its quite possible that is
true but I wanted to check somethings
> ecs<-testForDEU(ecs, formula0=formula0, formula1=formula1)
Testing for differential exon usage. (Progress report: one dot per 100
genes)
...........................................................................
.......
There were 50 or more warnings (use warnings() to see the first 50)
Warning messages:
1: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
Any idea what this warning is warning me about?
Both
> ecs<-estimatelog2FoldChanges(ecs)
Error in estimatelog2FoldChanges(ecs) :
fitExpToVar parameter is not in the design columns, double check
ecs at designColumns
and
> DEXSeqHTML( ecs, FDR=0.1, fitExpToVar="colony:strain" )
Error in DEXSeqHTML(ecs, FDR = 0.1, fitExpToVar = "colony:strain") :
fitExpToVar parameter is not in the design columns, double check
ecs at designColumns
> DEXSeqHTML( ecs, FDR=0.1, fitExpToVar="colony+strain" )
> ecs at designColumns
[1] "colony" "strain" "type"
I get that in DEXSeqHTML I'm asking for something that doesn't exist. Can
I visualise the results (if I have any) due to interaction? Whats going
wrong with estimatelog2FoldChanges?
Appreciate any light that can be shed
Eamonn
> ecs = read.HTSeqCounts(countfiles =
>file.path("/Users/eamonnmallon/Dropbox/Exon/wholegenome/textfiles",
>paste(rownames(samples), "txt", sep=".")), design = samples,flattenedfile
>= annotationfile)
> sampleNames(ecs) = rownames(samples)
>
> head( counts(ecs),20)
K61 K62 K63 K81 K82 K83 Q61 Q62 Q63 Q81 Q82
XLOC_000001:E001 27 34 23 28 18 33 36 22 49 39 23
XLOC_000001:E002 76 106 106 67 107 159 167 97 127 80 41
XLOC_000001:E003 1 0 2 1 0 0 1 4 0 3 1
XLOC_000001:E004 39 46 63 43 69 106 83 77 73 53 31
XLOC_000001:E005 0 0 0 0 0 0 0 0 0 0 0
XLOC_000002:E001 0 4 14 23 22 1 1 0 0 1 5
XLOC_000003:E001 23 22 10 9 12 32 26 35 24 10 13
XLOC_000004:E001 153 101 121 72 50 198 287 262 201 222 129
XLOC_000005:E001 8 13 6 7 2 11 15 23 20 15 18
XLOC_000006:E001 34 199 301 285 377 102 106 59 54 45 140
XLOC_000007:E001 0 0 2 9 4 0 0 0 0 0 0
XLOC_000008:E001 3 2 3 1 0 4 5 10 2 2 1
XLOC_000009:E001 315 137 87 112 244 185 311 288 179 320 63
XLOC_000010:E001 930 1059 445 317 191 171 485 570 738 865 52
XLOC_000011:E001 19 37 4 8 32 28 9 9 18 12 0
XLOC_000012:E001 1 5 1 2 6 11 27 0 0 4 1
XLOC_000013:E001 5 2 0 10 10 6 7 14 5 12 0
XLOC_000014:E001 1 3 1 0 1 3 9 1 8 0 3
XLOC_000015:E001 15 9 18 8 8 36 12 8 12 7 5
XLOC_000015:E002 15 15 21 12 8 45 30 11 17 13 18
> head(fData(ecs)[,c(1,2,9:12)])
geneID exonID chr start
end strand
XLOC_000001:E001 XLOC_000001 E001 gi|313870964|gb|AELG01010669.1| 2
577 .
XLOC_000001:E002 XLOC_000001 E002 gi|313870964|gb|AELG01010669.1| 656
773 .
XLOC_000001:E003 XLOC_000001 E003 gi|313870964|gb|AELG01010669.1| 870
875 .
XLOC_000001:E004 XLOC_000001 E004 gi|313870964|gb|AELG01010669.1| 1018
1087 .
XLOC_000001:E005 XLOC_000001 E005 gi|313870964|gb|AELG01010669.1| 1088
1088 .
XLOC_000002:E001 XLOC_000002 E001 gi|313870968|gb|AELG01010665.1| 85
662 .
> ecs<-estimateSizeFactors(ecs)
> sizeFactors(ecs)
K61 K62 K63 K81 K82 K83 Q61
Q62 Q63 Q81 Q82
0.7890977 1.4112509 1.4477112 0.8912231 1.0392289 1.0893970 1.2205867
0.8808826 0.9740058 1.0123368 0.7635365
> formuladispersion<-count~sample+(colony*strain)+exon #from anders email
> ecs<-estimateDispersions(ecs, formula=formuladispersion)
Dispersion estimation. (Progress report: one dot per 100 genes)
...........................................................................
.......
Warning messages:
1: In .local(object, ...) :
Exons with less than 10 counts will not be tested. For more details
please see the manual page of 'estimateDispersions', parameter 'minCount'
2: In .local(object, ...) :
Genes with more than 70 testable exons will be omitted from the
analysis. For more details please see the manual page of
'estimateDispersions', parameter 'maxExon'.
> ecs<-fitDispersionFunction(ecs)
> formula0<-count~sample+exon+(colony+strain)*exon
>
>formula1<-count~sample+exon+(colony+strain)*exon+(colony:strain)*I(exon==e
>xonID)
> ecs<-testForDEU(ecs, formula0=formula0, formula1=formula1)
Testing for differential exon usage. (Progress report: one dot per 100
genes)
...........................................................................
.......
There were 50 or more warnings (use warnings() to see the first 50)
> ecs<-estimatelog2FoldChanges(ecs)
Error in estimatelog2FoldChanges(ecs) :
fitExpToVar parameter is not in the design columns, double check
ecs at designColumns
> res1<-DEUresultTable(ecs)
> head(res1)
geneID exonID dispersion pvalue padjust
meanBase
XLOC_000001:E001 XLOC_000001 E001 9.330649e-02 0.1043954 0.8361045
29.695450
XLOC_000001:E002 XLOC_000001 E002 6.637378e-02 0.4553256 0.9404390
98.071027
XLOC_000001:E003 XLOC_000001 E003 9.960544e-01 0.3204769 0.9068146
1.218557
XLOC_000001:E004 XLOC_000001 E004 7.377259e-02 0.3339211 0.9104324
60.072384
XLOC_000001:E005 XLOC_000001 E005 1.000000e+08 NA NA
0.000000
XLOC_000002:E001 XLOC_000002 E001 2.382029e-01 NA NA
6.250462
> table(res1$padjust < 0.1)
FALSE
75341
> table ( tapply( res1$padjust < 0.1, geneIDs(ecs), any ) )
FALSE
4847
> sign <- res1[ which(res1$padjust < 0.1), ]
> DEXSeqHTML( ecs, FDR=0.1, fitExpToVar="colony:strain" )
Error in DEXSeqHTML(ecs, FDR = 0.1, fitExpToVar = "colony:strain") :
fitExpToVar parameter is not in the design columns, double check
ecs at designColumns
> DEXSeqHTML( ecs, FDR=0.1, fitExpToVar="colony+strain" )
Error in DEXSeqHTML(ecs, FDR = 0.1, fitExpToVar = "colony+strain") :
fitExpToVar parameter is not in the design columns, double check
ecs at designColumns
> ecs at designColumns
[1] "colony" "strain" "type"
> warnings()
Warning messages:
1: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
2: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
3: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
4: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
5: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
6: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
7: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
8: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
9: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
10: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
11: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
12: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
13: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
14: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
15: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
16: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
17: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
18: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
19: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
20: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
21: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
22: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
23: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
24: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
25: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
26: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
27: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
28: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
29: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
30: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
31: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
32: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
33: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
34: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
35: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
36: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
37: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
38: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
39: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
40: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
41: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
42: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
43: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
44: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
45: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
46: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
47: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
48: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
49: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
50: In chol.default(XVX + lambda * I, pivot = TRUE) :
the matrix is either rank-deficient or indefinite
Dr Eamonn Mallon
Lecturer in Evolutionary Biology
Adrian 220
Biology Department
University of Leicester
http://www2.le.ac.uk/departments/biology/people/mallon
On 12/06/2013 08:15, "Simon Anders" <anders at embl.de> wrote:
>On 06/06/13 16:22, Mallon, Eamonn B. (Dr.) wrote:
>> Dear all
>> I have 11 samples for a cross infection experiment where there are two
>>colonies (the hosts K or Q) and 2 strains (6 or 8).
>>
>> sample colony strain type
>> K61 k six single-read
>> K62 k six single-read
>> K63 k six single-read
>> K81 k eight single-read
>> K82 k eight single-read
>> K83 k eight single-read
>> Q61 q six single-read
>> Q62 q six single-read
>> Q63 q six single-read
>> Q81 q eight single-read
>> Q82 q eight single-read
>>
>> I am most interested in finding exon usage differences related to the
>>interaction of the colony and strain factors. Following the vignette I
>>put the following code together
>>
>>
>> formuladispersion<-count~sample+(colony:strain)+exon
>> ecs<-estimateDispersions(ecs, formula=formuladispersion)
>> ecs<-fitDispersionFunction(ecs)
>>
>> formula0<-count~sample+exon+(colony:strain)
>> formula1<-count~sample+exon+(colony:strain)*I(exon==exonID)
>> ecs<-testForDEU(ecs, formula0=formula0, formula1=formula1)
>>
>> Does this make sense?
>
>Not quite. This tests for main effects and interaction in one go, so it
>returns as hits all exons whose usage differs between colonies _and/or_
>between strains. You are looking for interactions, i.e., for exons whose
>usage is different between strains _and_ where this difference itself
>differs between colonies. So, you need
>
> formuladispersion<-count~sample+(colony*strain)+exon
>
> formula0<-count~sample+exon+(colony+strain)*exon
>
>formula1<-count~sample+exon+(colony+strain)*exon+(colony:strain)*I(exon==e
>xonID)
>
>Things look a bit simpler if you use the new "TRT" functions:
>
> formuladispersion<-count~sample+(colony*strain)+exon
>
> formula0<-count~sample+exon+(colony+strain)*exon
> formula1<-count~sample+exon+(colony*strain)*exon
>
> Simon
>
>_______________________________________________
>Bioconductor mailing list
>Bioconductor at r-project.org
>https://stat.ethz.ch/mailman/listinfo/bioconductor
>Search the archives:
>http://news.gmane.org/gmane.science.biology.informatics.conductor
More information about the Bioconductor
mailing list