how to make makeContrast in edgeR?
2
0
Entering edit mode
wang peter ★ 2.0k
@wang-peter-4647
Last seen 10.3 years ago
I meet such problem, this is the coding raw.data <- read.table("expression-table.txt",row.names=1) lib_size <- read.table("lib_size.txt"); lib_size <- unlist(lib_size) d <- raw.data[, 2:dim(raw.data)[2]] length<-raw.data[, 1] d <- DGEList(counts = d, lib.size = lib_size) #normalization dge <- calcNormFactors(dge) treatment=factor(c(rep('control',6),rep('treated',24),rep('control',5) )) time=factor(c('0h','0h','0h','24h','24h','24h','0h','0h','0h','6h','6h ','6h','6h','12h','12h','12h','12h','18h','18h','18h','18h', '24h','24h','24h','36h','36h','36h','48h','48h','48h','6h ','12h','18h','36h','48h')) design <- model.matrix(~treatment*time) dge <- estimateGLMCommonDisp(dge, design) dge <- estimateGLMTagwiseDisp(dge, design) glmfit.dge <- glmFit(dge, design, dispersion=dge$tagwise.dispersion) i want to make contrast between 6h_treated vs 6h_contol, 12h_treated vs 12h_contol, and so on.... how to use makeContrast or any other function to compare them thank you very much -- shan gao Room 231(Dr.Fei lab) Boyce Thompson Institute Cornell University Tower Road, Ithaca, NY 14853-1801 Office phone: 1-607-254-1267(day) Official email:sg839 at cornell.edu Facebook:http://www.facebook.com/profile.php?id=100001986532253
• 5.1k views
ADD COMMENT •
0
Entering edit mode
KJ Lim ▴ 420
@kj-lim-5288
Last seen 4.3 years ago
Finland
Dear Shan Gao, You could consider import your experiment to edgeR from a target file(which contains information of your experiment). It is much easier for you to track back anything wrong in between if there was one or two. The Section 3.3 of the latest edgeR User's Guide is worth for you to have a look. http://bioconductor.org/packages/2.11/bioc/vignettes/edgeR/inst/doc/ed geRUsersGuide.pdf You can find the guide how to utilize the makeContrast function. Hopefully this help. Best regards, KJ Lim On 25 August 2012 07:15, wang peter <wng.peter@gmail.com> wrote: > I meet such problem, this is the coding > > raw.data <- read.table("expression-table.txt",row.names=1) > lib_size <- read.table("lib_size.txt"); > lib_size <- unlist(lib_size) > d <- raw.data[, 2:dim(raw.data)[2]] > length<-raw.data[, 1] > > d <- DGEList(counts = d, lib.size = lib_size) > #normalization > dge <- calcNormFactors(dge) > treatment=factor(c(rep('control',6),rep('treated',24),rep('control', 5))) > > time=factor(c('0h','0h','0h','24h','24h','24h','0h','0h','0h','6h',' 6h','6h','6h','12h','12h','12h','12h','18h','18h','18h','18h', > > '24h','24h','24h','36h','36h','36h','48h','48h','48h','6h','12h','1 8h','36h','48h')) > design <- model.matrix(~treatment*time) > > dge <- estimateGLMCommonDisp(dge, design) > dge <- estimateGLMTagwiseDisp(dge, design) > glmfit.dge <- glmFit(dge, design, dispersion=dge$tagwise.dispersion) > > i want to make contrast between 6h_treated vs 6h_contol, 12h_treated > vs 12h_contol, > and so on.... > > how to use makeContrast or any other function to compare them > > thank you very much > > -- > shan gao > Room 231(Dr.Fei lab) > Boyce Thompson Institute > Cornell University > Tower Road, Ithaca, NY 14853-1801 > Office phone: 1-607-254-1267(day) > Official email:sg839@cornell.edu > Facebook:http://www.facebook.com/profile.php?id=100001986532253 > > _______________________________________________ > Bioconductor mailing list > Bioconductor@r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > [[alternative HTML version deleted]]
ADD COMMENT •
0
Entering edit mode
thank you for your suggestion but it is not relevant to my question if you donot use a target file, and write some coding following my current coding. can i do those comparitions? shan
ADD REPLY •
0
Entering edit mode
@sean-davis-490
Last seen 4 months ago
United States
On Sat, Aug 25, 2012 at 12:15 AM, wang peter <wng.peter@gmail.com> wrote: > I meet such problem, this is the coding > > raw.data <- read.table("expression-table.txt",row.names=1) > lib_size <- read.table("lib_size.txt"); > lib_size <- unlist(lib_size) > d <- raw.data[, 2:dim(raw.data)[2]] > length<-raw.data[, 1] > > d <- DGEList(counts = d, lib.size = lib_size) > #normalization > dge <- calcNormFactors(dge) > treatment=factor(c(rep('control',6),rep('treated',24),rep('control', 5))) > > time=factor(c('0h','0h','0h','24h','24h','24h','0h','0h','0h','6h',' 6h','6h','6h','12h','12h','12h','12h','18h','18h','18h','18h', > > '24h','24h','24h','36h','36h','36h','48h','48h','48h','6h','12h','1 8h','36h','48h')) > design <- model.matrix(~treatment*time) > > dge <- estimateGLMCommonDisp(dge, design) > dge <- estimateGLMTagwiseDisp(dge, design) > glmfit.dge <- glmFit(dge, design, dispersion=dge$tagwise.dispersion) > > i want to make contrast between 6h_treated vs 6h_contol, 12h_treated > vs 12h_contol, > and so on.... > > how to use makeContrast or any other function to compare them > > I'd suggest making a single factor, call it fac, derived by combining the time and treatment variables so that your factor looks like: control_0h, control_0h, ..., treated_0h, ..., treated_48h, etc. design = model.matrix( ~ fac) contmat = makeContrasts(treated_6h-control_6h,treated_12h-control_12h) I didn't test this, but this is a simple approach to the problem that often works nicely for two-factor designs. Sean [[alternative HTML version deleted]]
ADD COMMENT •
0
Entering edit mode
thank you very much that is really a simpler way to do so but if i want to follow the way i used, how to use the funtion ? i think there is must a way to do so by two factor design shan
ADD REPLY •
0
Entering edit mode
it also generate another problem if i use one factor, as you suggested : control_0h, control_0h, ..., treated_0h, ..., treated_48h, etc. design = model.matrix( ~ fac) contmat = makeContrasts(treated_6h-control_6h,treated_12h-control_12h) will the glm use all the samples to do regression for single comparision? i hope glm use all to do regression, even on single test like:treated_6h-control_6h
ADD REPLY •
0
Entering edit mode
On Sat, Aug 25, 2012 at 4:31 PM, wang peter <wng.peter@gmail.com> wrote: > it also generate another problem > > if i use one factor, as you suggested : control_0h, control_0h, ..., > treated_0h, ..., treated_48h, etc. > > design = model.matrix( ~ fac) > contmat = makeContrasts(treated_6h-control_6h,treated_12h- control_12h) > > will the glm use all the samples to do regression for single comparision? > > I believe it will "do the right thing", yes. Sean > i hope glm use all to do regression, even on single test > like:treated_6h-control_6h > > _______________________________________________ > Bioconductor mailing list > Bioconductor@r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > [[alternative HTML version deleted]]
ADD REPLY •
0
Entering edit mode
contrasts.fit cannnot work on glmFit such is my coding: d <- DGEList(counts = d, lib.size = lib_size) #normalization dge <- calcNormFactors(d) rm(d) fac=factor(c('c0h','c0h','c0h','c24h','c24h','c24h','t0h','t0h','t0h', 't6h','t6h','t6h','t6h','t12h','t12h','t12h','t12h','t18h','t18h','t18 h','t18h', 't24h','t24h','t24h','t36h','t36h','t36h','t48h','t48h',' t48h','c6h','c12h','c18h','c36h','c48h')) design = model.matrix(~fac) colnames(design) <- levels(fac) contrast.matrix = makeContrasts("t0h-c0h","t6h-c6h","t12h-c12h","t18h-c18h","t24h-c24h ","t36h-c36h","t48h-c48h","t48h-c0h",levels=design) dge <- estimateGLMCommonDisp(dge, design)# dge <- estimateGLMTagwiseDisp(dge, design) glmfit.dge <- glmFit(dge, design, dispersion=dge$tagwise.dispersion)# fit1 <- contrasts.fit(glmfit.dge, contrast.matrix) Error in cov2cor(fit$cov.coefficients) : 'V' is not a square numeric matrix fit2 <- eBayes(fit1) topTable(fit2, coef=1, adjust="BH")
ADD REPLY •
0
Entering edit mode
Dear Gao Shan, > fit1 <- contrasts.fit(glmfit.dge, contrast.matrix) > fit2 <- eBayes(fit1) > topTable(fit2, coef=1, adjust="BH") It seemed that you are not going to carry out your analysis with the GLM likelihood ratio test of edgeR package. If I'm not mistaken, the *eBayes*function is belongs to the Limma package. If you would like to analyse your RNA-Seq data with Limma package, you need to do data conversion with voom() as mentioned in the Limma User's Guide. The Limma User's Guide is worth for a read in this case. > makeContrasts("t0h-c0h","t6h-c6h","t12h-c12h","t18h-c18h","t24h-c24h ","t36h-c36h","t48h-c48h","t48h-c0h",levels=design) The section 3.3.1 of Chapter 3 in the latest edgeR User's Guide has a clear guide how you could build your contrast matrix if you still prefer to carry your analysis with the edgeR. That section was added recently because several users have been asking similar question as you did via mailing list . Have a nice Sunday. Best regards, KJ Lim On 26 August 2012 05:35, wang peter <wng.peter@gmail.com> wrote: > > contrasts.fit cannnot work on glmFit > such is my coding: > > d <- DGEList(counts = d, lib.size = lib_size) > #normalization > dge <- calcNormFactors(d) > rm(d) > > fac=factor(c('c0h','c0h','c0h','c24h','c24h','c24h','t0h','t0h','t0h', 't6h','t6h','t6h','t6h','t12h','t12h','t12h','t12h','t18h','t18h','t18 h','t18h', > 't24h','t24h','t24h','t36h','t36h','t36h','t48h','t48h','t48h','c6h', 'c12h','c18h','c36h','c48h')) > design = model.matrix(~fac) > colnames(design) <- levels(fac) > contrast.matrix = > makeContrasts("t0h-c0h","t6h-c6h","t12h-c12h","t18h-c18h","t24h-c24h ","t36h-c36h","t48h-c48h","t48h-c0h",levels=design) > > dge <- estimateGLMCommonDisp(dge, design)# > dge <- estimateGLMTagwiseDisp(dge, design) > glmfit.dge <- glmFit(dge, design, dispersion=dge$tagwise.dispersion)# > > fit1 <- contrasts.fit(glmfit.dge, contrast.matrix) > > Error in cov2cor(fit$cov.coefficients) : > 'V' is not a square numeric matrix > > fit2 <- eBayes(fit1) > topTable(fit2, coef=1, adjust="BH") > > _______________________________________________ > Bioconductor mailing list > Bioconductor@r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor [[alternative HTML version deleted]]
ADD REPLY •
0
Entering edit mode
Dear lim, i read the manual and did such work > d <- DGEList(counts = d, lib.size = lib_size) > #normalization > dge <- calcNormFactors(d) > rm(d) > > fac=factor(c('c0h','c0h','c0h','c24h','c24h','c24h','t0h','t0h','t0h ','t6h','t6h','t6h','t6h','t12h','t12h','t12h','t12h','t18h','t18h','t 18h','t18h', + 't24h','t24h','t24h','t36h','t36h','t36h','t48h','t48h' ,'t48h','c6h','c12h','c18h','c36h','c48h')) > design = model.matrix(~fac) > colnames(design) <- levels(fac) > > dge <- estimateGLMCommonDisp(dge, design) > dge <- estimateGLMTagwiseDisp(dge, design) > glmfit.dge <- glmFit(dge, design, dispersion=dge$tagwise.dispersion) > aa = makeContrasts("t0h-c0h",levels=design) > lrt.dgh <- glmLRT(dge, glmfit.dge, aa) Error in coef.names[coef] : only 0's may be mixed with negative subscripts > sessionInfo() R version 2.15.1 (2012-06-22) Platform: x86_64-pc-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=Chinese (Simplified)_People's Republic of China.936 [2] LC_CTYPE=Chinese (Simplified)_People's Republic of China.936 [3] LC_MONETARY=Chinese (Simplified)_People's Republic of China.936 [4] LC_NUMERIC=C [5] LC_TIME=Chinese (Simplified)_People's Republic of China.936 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] edgeR_2.6.10 limma_3.12.1 loaded via a namespace (and not attached): [1] tools_2.15.1
ADD REPLY •
0
Entering edit mode
thank you lim, i read the manual and write such coding but meet an error > d <- DGEList(counts = d, lib.size = lib_size) > #normalization > dge <- calcNormFactors(d) > rm(d) > > fac=factor(c('c0h','c0h','c0h','c24h','c24h','c24h','t0h','t0h','t0h ','t6h','t6h','t6h','t6h','t12h','t12h','t12h','t12h','t18h','t18h','t 18h','t18h', + 't24h','t24h','t24h','t36h','t36h','t36h','t48h','t48h' ,'t48h','c6h','c12h','c18h','c36h','c48h')) > design = model.matrix(~fac) > colnames(design) <- levels(fac) > dge <- estimateGLMCommonDisp(dge, design)# > dge <- estimateGLMTagwiseDisp(dge, design) > glmfit.dge <- glmFit(dge, design, dispersion=dge$tagwise.dispersion)# > t0_c0 = makeContrasts("t0h-c0h",levels=design) > lrt.dge <- glmLRT(dge, glmfit.dge, contrast = t0_c0) > result <- topTags(lrt.dge, n=dim(dge)[1], adjust.method="BH", sort.by="p.value") Error in abs(object$table$logFC) : Non-numeric argument to mathematical function R version 2.15.1 (2012-06-22) Platform: x86_64-pc-mingw32/x64 (64-bit) locale: [1] LC_COLLATE=Chinese (Simplified)_People's Republic of China.936 [2] LC_CTYPE=Chinese (Simplified)_People's Republic of China.936 [3] LC_MONETARY=Chinese (Simplified)_People's Republic of China.936 [4] LC_NUMERIC=C [5] LC_TIME=Chinese (Simplified)_People's Republic of China.936 attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] edgeR_2.6.10 limma_3.12.1
ADD REPLY •
0
Entering edit mode
Dear Gao Shan, Kindly please spare time to read through the Chapter 3 from beginning till the end in the latest edgeR User's guide. The guides in that Chapter are useful to your analysis. > design = model.matrix(~fac) This design matrix will include an intercept column in your design matrix. If you would like to utilize the makeContrasts function for glmLRT analysis and to avoid confusion, you could modify your design matrix: > design = model.matrix(~0+fac) The contrast vector could be contracted: > t0_c0 = makeContrasts(t0h-c0h,levels=design) ## assuming label of design matrix column is right > lrt.dge <- glmLRT(dge, glmfit.dge, contrast = t0_c0) There is a poster worth for a read either: https://www.stat.math.ethz.ch/pipermail/bioconductor/2012-May/045700.h tml Best regards, KJ Lim On 26 August 2012 22:53, wang peter <wng.peter at="" gmail.com=""> wrote: > > thank you lim, > i read the manual and write such coding > but meet an error > > > d <- DGEList(counts = d, lib.size = lib_size) > > #normalization > > dge <- calcNormFactors(d) > > rm(d) > > > > fac=factor(c('c0h','c0h','c0h','c24h','c24h','c24h','t0h','t0h','t 0h','t6h','t6h','t6h','t6h','t12h','t12h','t12h','t12h','t18h','t18h', 't18h','t18h', > + 't24h','t24h','t24h','t36h','t36h','t36h','t48h','t48 h','t48h','c6h','c12h','c18h','c36h','c48h')) > > design = model.matrix(~fac) > > colnames(design) <- levels(fac) > > dge <- estimateGLMCommonDisp(dge, design)# > > dge <- estimateGLMTagwiseDisp(dge, design) > > glmfit.dge <- glmFit(dge, design, dispersion=dge$tagwise.dispersion)# > > t0_c0 = makeContrasts("t0h-c0h",levels=design) > > lrt.dge <- glmLRT(dge, glmfit.dge, contrast = t0_c0) > > result <- topTags(lrt.dge, n=dim(dge)[1], adjust.method="BH", sort.by="p.value") > > Error in abs(object$table$logFC) : > Non-numeric argument to mathematical function > > R version 2.15.1 (2012-06-22) > Platform: x86_64-pc-mingw32/x64 (64-bit) > > locale: > [1] LC_COLLATE=Chinese (Simplified)_People's Republic of China.936 > [2] LC_CTYPE=Chinese (Simplified)_People's Republic of China.936 > [3] LC_MONETARY=Chinese (Simplified)_People's Republic of China.936 > [4] LC_NUMERIC=C > [5] LC_TIME=Chinese (Simplified)_People's Republic of China.936 > > attached base packages: > [1] stats graphics grDevices utils datasets methods base > > other attached packages: > [1] edgeR_2.6.10 limma_3.12.1
ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 406 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6