Can any one help with the below? Dr Eamonn Mallon Lecturer in Evolutionary Biology Adrian 220 Biology Department University of Leicester http://www2.le.ac.uk/departments/biology/people/mallon On 07/06/2013 14:30, "Mallon, Eamonn B. (Dr.)" <ebm3 at="" mail.cfs.le.ac.uk=""> wrote: >My initial question still stands, but I have also tried a simpler model by >just looking at the effect of strain on one colony. > >So there are two levels to condition (strain = 6 or 8) > >It runs fine but I get >> table(res1$padjust < 0.1) > >FALSE TRUE >65836 5987 > > >This seems excessive > >And the warning (50 times) > >> DEXSeqHTML( ecs, FDR=0.1, ) >There were 50 or more warnings (use warnings() to see the first 50) >> warnings() >Warning messages: >1: In chol.default(XVX + lambda * I, pivot = TRUE) : > the matrix is either rank-deficient or indefinite > > > >Any ideas what is happening > >Thanks in advance > >Eamonn > >Dr Eamonn Mallon >Lecturer in Evolutionary Biology >Adrian 220 >Biology Department >University of Leicester > >http://www2.le.ac.uk/departments/biology/people/mallon > > > > > >On 06/06/2013 15:22, "Mallon, Eamonn B. (Dr.)" <ebm3 at="" leicester.ac.uk=""> >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? >> >>Thanks in advance >>Eamonn >> >>Dr Eamonn Mallon >>Lecturer in Evolutionary Biology >>Adrian 220 >>Biology Department >>University of Leicester >> >>http://www2.le.ac.uk/departments/biology/people/mallon >> >> >> [[alternative HTML version deleted]] >> >>_______________________________________________ >>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 >

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(ex on==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(ex on==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

Just looking over my own plot and noticed a typo: On 12/06/13 09:15, Simon Anders wrote: > 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 ^ should be "*exon" > > formula0<-count~sample+exon+(colony+strain)*exon > > formula1<-count~sample+exon+(colony+strain)*exon+(colony:strain)*I(e xon==exonID) > > > Things look a bit simpler if you use the new "TRT" functions: > > formuladispersion<-count~sample+(colony*strain)+exon ^ here too > > formula0<-count~sample+exon+(colony+strain)*exon > formula1<-count~sample+exon+(colony*strain)*exon