edgeR: adjusting prior.n
1
0
Entering edit mode
Guest User ★ 13k
@guest-user-4897
Last seen 10.3 years ago
Hello Bioconductor mailing list, I'm using edgeR to analyze RNA-seq and RIP-seq data. I'm using the moderated gene-wise dispersion method where I choose a prior.n based upon the number of samples I have and how much shrinkage I desire to reduce variance in my data. I've followed the R user guide and a guide a professor made at my university to achieve gene-wise (tag-wise) dispersion as opposed to common dispersion. I can see that the tag- wise dispersion changes for each gene when I change the prior.n value (see below). However I notice that that p-values and logFC values don't change when I change the prior.n (see below). So I'm wondering if I'm doing something wrong or if the p-values and/or logFC are not supposed to change? Or is it just some data is not effected at all by prior.n? The following is an example of my command line code analyzing some RNA-seq data. I illustrated an example where I used a prior.n = 1 or 4 to illustrate how I don't see a change in the final results. Is this correct? Thanks J normFact = calcNormFactors(TotalReads) treatments = substr(colnames(TotalReads), 1, 3) d = DGEList(counts = TotalReads, group = treatments, genes = rownames(TotalReads), lib.size = colSums(TotalReads) * normFact) d = estimateCommonDisp(d) d$common.dispersion [1] 0.1669584 #comparing prior.n 1 v. 4 d1 = estimateTagwiseDisp(d, prior.n = 1) d1 $tagwise.dispersion [1] 0.1607020 0.3674590 0.1058665 0.1080470 0.1545314 20684 more elements ??? d4 = estimateTagwiseDisp(d, prior.n = 4) d4 $tagwise.dispersion [1] 0.1463446 0.4059551 0.1278219 0.1213169 0.2150826 20684 more elements ... #prior.n = 1 edgeHSMvHSFd1 = exactTest(d1, pair = c("HSM", "HSF"), (common.disp = FALSE )) head(edgeHSMvHSFd1) An object of class "DGEExact" $table logFC logCPM PValue ENSG00000000003 -0.07130289 6.290327 6.073496e-01 ENSG00000000005 -2.28060195 -4.439294 1.000000e+00 ENSG00000000419 0.25083103 4.119127 8.780027e-02 ENSG00000000457 -0.19453270 4.557886 8.468506e-02 ENSG00000000460 -0.05015655 1.770160 9.075606e-01 ENSG00000000938 0.92221495 4.526491 1.227442e-11 #prior.n = 4 edgeHSMvHSFd4 = exactTest(d4, pair = c("HSM", "HSF"), (common.disp = FALSE )) head(edgeHSMvHSFd4) An object of class "DGEExact" $table logFC logCPM PValue ENSG00000000003 -0.07130289 6.290327 6.073496e-01 ENSG00000000005 -2.28060195 -4.439294 1.000000e+00 ENSG00000000419 0.25083103 4.119127 8.780027e-02 ENSG00000000457 -0.19453270 4.557886 8.468506e-02 ENSG00000000460 -0.05015655 1.770160 9.075606e-01 ENSG00000000938 0.92221495 4.526491 1.227442e-11 sum(edgeHSMvHSFd1$table[,1]>0) [1] 12770 sum(edgeHSMvHSFd4$table[,1]>0) [1] 12770 #same results for DE genes for both sum(edgeHSMvHSFd1$table[,1]>0 & edgeHSMvHSFd1$table[,3]<0.01) [1] 1763 sum(edgeHSMvHSFd4$table[,1]>0 & edgeHSMvHSFd4$table[,3]<0.01) [1] 1763 -- output of sessionInfo(): > sessionInfo() R version 2.15.0 (2012-03-30) Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods [7] base other attached packages: [1] edgeR_2.6.2 limma_3.12.0 loaded via a namespace (and not attached): [1] tools_2.15.0 -- Sent via the guest posting facility at bioconductor.org.
edgeR edgeR • 1.9k views
ADD COMMENT •
0
Entering edit mode
@gordon-smyth
Last seen 11 hours ago
WEHI, Melbourne, Australia

Dear Jason,

The argument prior.n was removed from the estimateTagwiseDisp() function quite a while ago. The current version of edgeR is 3.4.2:

http://www.bioconductor.org/packages/release/bioc/html/edgeR.html

When prior.n existed, it certainly did affect the tagwise dispersions, logFC and p-values.

You have made a programming error by enclosing the argument common.dispersion=FALSE in parentheses, which prevents correct argument matching.

Best wishes
Gordon
ADD COMMENT •
0
Entering edit mode
Hi, Dr. Gordon and all: I used limma before, but not much on limma-voom package. I run into some issue of warnings and then errors, not sure what is. Below are my commands, only show essential parts: >show(head(targets)) ID RasType Subject Type 1 T4626_01A KRas TCGA_38_4626 Tumor 2 N4626_11A KRas TCGA_38_4626 Normal 3 T2668_01A KRas TCGA_44_2668 Tumor 4 N2668_11A KRas TCGA_44_2668 Normal 5 T6145_01A KRas TCGA_44_6145 Tumor > RasTum<-paste(targets$RasType, targets$Type,sep=".") > RasTum<-factor(RasTum, levels=c("KRas.Tumor","KRas.Normal","KRasP.Tumor","KRasP.Normal")); > RasTum [1] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor [6] KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal [11] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor [16] KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal [21] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRasP.Tumor [26] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal [31] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor [36] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal [41] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor [46] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal [51] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor [56] KRasP.Normal KRasP.Tumor KRasP.Normal Levels: KRas.Tumor KRas.Normal KRasP.Tumor KRasP.Normal >Patient<-factor(targets$Subject); >rawdata<-rawdata[,c("symbols","entrezID" ,"gene_id",targets$ID)]; >y<-DGEList(count=rawdata[,4:ncol(rawdata)],genes=rawdata[,1:3]); > show(dim(y)); [1] 20531 58 > y3<-y[rowSums(y$counts)>=50,]; > show(dim(y3)); [1] 18215 58 > y3$samples$lib.size <- colSums(y3$counts) > rownames(y3$counts) <- rownames(y3$genes) <- y3$genes$entrezID > y3 <- calcNormFactors(y3,method=c("TMM")); > show(y3$samples[1:6,]) group lib.size norm.factors T4626_01A 1 97435132 0.9834587 N4626_11A 1 62583672 0.9154837 T2668_01A 1 77034746 0.9687860 N2668_11A 1 58631311 0.8644352 T6145_01A 1 31634419 1.0303445 N6145_11A 1 29872761 0.9615021 > design<-model.matrix(~Patient+RasTum); > show(design[1:6,c(1:10,30:32)]) (Intercept) PatientTCGA_38_4626 PatientTCGA_38_4627 PatientTCGA_38_4632 1 1 1 0 0 2 1 1 0 0 3 1 0 0 0 4 1 0 0 0 5 1 0 0 0 6 1 0 0 0 PatientTCGA_44_2655 PatientTCGA_44_2657 PatientTCGA_44_2661 1 0 0 0 2 0 0 0 3 0 0 0 4 0 0 0 5 0 0 0 6 0 0 0 PatientTCGA_44_2662 PatientTCGA_44_2665 PatientTCGA_44_2668 RasTumKRas.Normal 1 0 0 0 0 2 0 0 0 1 3 0 0 1 0 4 0 0 1 1 5 0 0 0 0 6 0 0 0 1 RasTumKRasP.Tumor RasTumKRasP.Normal 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 > colnames(design)<-sub("Patient","",colnames(design)); > colnames(design)<-sub("RasTum","",colnames(design)); > colnames(design)[1] <- "Intercept" > show(design[1:6,c(1:10,30:32)]) Intercept TCGA_38_4626 TCGA_38_4627 TCGA_38_4632 TCGA_44_2655 TCGA_44_2657 1 1 1 0 0 0 0 2 1 1 0 0 0 0 3 1 0 0 0 0 0 4 1 0 0 0 0 0 5 1 0 0 0 0 0 6 1 0 0 0 0 0 TCGA_44_2661 TCGA_44_2662 TCGA_44_2665 TCGA_44_2668 KRas.Normal KRasP.Tumor 1 0 0 0 0 0 0 2 0 0 0 0 1 0 3 0 0 0 1 0 0 4 0 0 0 1 1 0 5 0 0 0 0 0 0 6 0 0 0 0 1 0 KRasP.Normal 1 0 2 0 3 0 4 0 5 0 6 0 > dim(design) [1] 58 32 > con.matrix<-makeContrasts( + KRas.Tumor_KRas.Normal = -KRas.Normal, + KRasP.Tumor_KRasP.Normal = KRasP.Tumor-KRasP.Normal, + KRas.Tumor_KRasP.Tumor=-KRasP.Tumor, + levels=design) > v <- voom(y3,design); Coefficients not estimable: KRasP.Normal Warning message: Partial NA coefficients for 18215 probe(s) > fit<- lmFit(v,design) Coefficients not estimable: KRasP.Normal Warning message: Partial NA coefficients for 18215 probe(s) > fit1<-contrasts.fit(fit, con.matrix) Error in contrasts.fit(fit, con.matrix) : trying to take contrast of non-estimable coefficient I used a similar way before for array data, but not with limma-voom for RNAseq data. Any suggestions here? Thanks a lot in advance! Ming ATRF/NCI-Frederick Maryland, USA Here is the sessionInfo: > show(sessionInfo()) R version 3.0.1 (2013-05-16) Platform: x86_64-unknown-linux-gnu (64-bit) locale: [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 [7] LC_PAPER=C LC_NAME=C [9] LC_ADDRESS=C LC_TELEPHONE=C [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] edgeR_3.4.2 limma_3.16.8 [[alternative HTML version deleted]]
ADD REPLY •
0
Entering edit mode

Hi Ming,

This has nothing to do with voom. It is the same message that you will always receive from limma when you try to fit an over-parametrized model. (By comparison, edgeR would have simply told you the design matrix was not of full rank and stopped.)

You can only used a paired analysis (+Patient etc) when the treatments to be compared are applied to the same patient. You can't compare different tumour types (KRas and KRasP) which affect different patients.

Gordon

ADD REPLY •
0
Entering edit mode
Hi, Dr. Gordon: Thx a lot for pointing that out. I just new to voom and so keep thinking about issue might come from voom. I shall be aware of that issue in limma. It works now. By the way, a small quick question: in the topTable result of limma, the logFC is base 2 or is natural base (2.71828 etc)? I though was base 2, just want to confirm. Thx again and best Ming > Date: Wed, 19 Feb 2014 17:32:06 +1100 > From: smyth@wehi.EDU.AU > To: yi02@hotmail.com > CC: bioconductor@r-project.org > Subject: Re: Problem in limma-voom > > Hi Ming, > > This has nothing to do with voom. It is the same message that you will > always receive from limma when you try to fit an over-parametrized model. > (By comparison, edgeR would have simply told you the design matrix was not > of full rank and stopped.) > > You can only used a paired analysis (+Patient etc) when the treatments to > be compared are applied to the same patient. You can't compare different > tumour types (KRas and KRasP) which affect different patients. > > Gordon > > > On Wed, 19 Feb 2014, Ming Yi wrote: > > > Hi, Dr. Gordon and all: > > > > > > > > I used limma before, but not much on limma-voom package. I run into some > > issue of warnings and then errors, not sure what is. Below are my > > commands, only show essential parts: > > > > > > > >> show(head(targets)) > > > > ID RasType Subject Type > > 1 T4626_01A KRas TCGA_38_4626 Tumor > > 2 N4626_11A KRas TCGA_38_4626 Normal > > 3 T2668_01A KRas TCGA_44_2668 Tumor > > 4 N2668_11A KRas TCGA_44_2668 Normal > > 5 T6145_01A KRas TCGA_44_6145 Tumor > > > > > RasTum<-paste(targets$RasType, targets$Type,sep=".") > >> RasTum<-factor(RasTum, levels=c("KRas.Tumor","KRas.Normal","KRasP.Tumor","KRasP.Normal")); > >> RasTum > > [1] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor > > [6] KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal > > [11] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor > > [16] KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal > > [21] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRasP.Tumor > > [26] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal > > [31] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor > > [36] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal > > [41] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor > > [46] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal > > [51] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor > > [56] KRasP.Normal KRasP.Tumor KRasP.Normal > > Levels: KRas.Tumor KRas.Normal KRasP.Tumor KRasP.Normal > > > >> Patient<-factor(targets$Subject); > > > >> rawdata<-rawdata[,c("symbols","entrezID" ,"gene_id",targets$ID)]; > > > >> y<-DGEList(count=rawdata[,4:ncol(rawdata)],genes=rawdata[,1:3]); > > > >> show(dim(y)); > > [1] 20531 58 > > > >> y3<-y[rowSums(y$counts)>=50,]; > > > >> show(dim(y3)); > > [1] 18215 58 > >> y3$samples$lib.size <- colSums(y3$counts) > >> rownames(y3$counts) <- rownames(y3$genes) <- y3$genes$entrezID > >> y3 <- calcNormFactors(y3,method=c("TMM")); > >> show(y3$samples[1:6,]) > > group lib.size norm.factors > > T4626_01A 1 97435132 0.9834587 > > N4626_11A 1 62583672 0.9154837 > > T2668_01A 1 77034746 0.9687860 > > N2668_11A 1 58631311 0.8644352 > > T6145_01A 1 31634419 1.0303445 > > N6145_11A 1 29872761 0.9615021 > >> design<-model.matrix(~Patient+RasTum); > >> show(design[1:6,c(1:10,30:32)]) > > (Intercept) PatientTCGA_38_4626 PatientTCGA_38_4627 PatientTCGA_38_4632 > > 1 1 1 0 0 > > 2 1 1 0 0 > > 3 1 0 0 0 > > 4 1 0 0 0 > > 5 1 0 0 0 > > 6 1 0 0 0 > > PatientTCGA_44_2655 PatientTCGA_44_2657 PatientTCGA_44_2661 > > 1 0 0 0 > > 2 0 0 0 > > 3 0 0 0 > > 4 0 0 0 > > 5 0 0 0 > > 6 0 0 0 > > PatientTCGA_44_2662 PatientTCGA_44_2665 PatientTCGA_44_2668 RasTumKRas.Normal > > 1 0 0 0 0 > > 2 0 0 0 1 > > 3 0 0 1 0 > > 4 0 0 1 1 > > 5 0 0 0 0 > > 6 0 0 0 1 > > RasTumKRasP.Tumor RasTumKRasP.Normal > > 1 0 0 > > 2 0 0 > > 3 0 0 > > 4 0 0 > > 5 0 0 > > 6 0 0 > >> colnames(design)<-sub("Patient","",colnames(design)); > >> colnames(design)<-sub("RasTum","",colnames(design)); > >> colnames(design)[1] <- "Intercept" > >> show(design[1:6,c(1:10,30:32)]) > > Intercept TCGA_38_4626 TCGA_38_4627 TCGA_38_4632 TCGA_44_2655 TCGA_44_2657 > > 1 1 1 0 0 0 0 > > 2 1 1 0 0 0 0 > > 3 1 0 0 0 0 0 > > 4 1 0 0 0 0 0 > > 5 1 0 0 0 0 0 > > 6 1 0 0 0 0 0 > > TCGA_44_2661 TCGA_44_2662 TCGA_44_2665 TCGA_44_2668 KRas.Normal KRasP.Tumor > > 1 0 0 0 0 0 0 > > 2 0 0 0 0 1 0 > > 3 0 0 0 1 0 0 > > 4 0 0 0 1 1 0 > > 5 0 0 0 0 0 0 > > 6 0 0 0 0 1 0 > > KRasP.Normal > > 1 0 > > 2 0 > > 3 0 > > 4 0 > > 5 0 > > 6 0 > >> dim(design) > > [1] 58 32 > >> con.matrix<-makeContrasts( > > + KRas.Tumor_KRas.Normal = -KRas.Normal, > > + KRasP.Tumor_KRasP.Normal = KRasP.Tumor-KRasP.Normal, > > + KRas.Tumor_KRasP.Tumor=-KRasP.Tumor, > > + levels=design) > >> v <- voom(y3,design); > > Coefficients not estimable: KRasP.Normal > > Warning message: > > Partial NA coefficients for 18215 probe(s) > > > >> fit<- lmFit(v,design) > > Coefficients not estimable: KRasP.Normal > > Warning message: > > Partial NA coefficients for 18215 probe(s) > >> fit1<-contrasts.fit(fit, con.matrix) > > Error in contrasts.fit(fit, con.matrix) : > > trying to take contrast of non-estimable coefficient > > > > > > I used a similar way before for array data, but not with limma- voom for RNAseq data. Any suggestions here? > > > > > > > > Thanks a lot in advance! > > > > > > > > Ming > > > > > > > > ATRF/NCI-Frederick > > > > Maryland, USA > > > > > > > > Here is the sessionInfo: > > > >> show(sessionInfo()) > > R version 3.0.1 (2013-05-16) > > Platform: x86_64-unknown-linux-gnu (64-bit) > > > > locale: > > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > > [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 > > [7] LC_PAPER=C LC_NAME=C > > [9] LC_ADDRESS=C LC_TELEPHONE=C > > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > > > > attached base packages: > > [1] stats graphics grDevices utils datasets methods base > > > > other attached packages: > > [1] edgeR_3.4.2 limma_3.16.8 > > > > > > > > > > > > > > > > > > > > > > ______________________________________________________________________ > The information in this email is confidential and inte...{{dropped:9}}
ADD REPLY •
0
Entering edit mode
Hi, On Wed, Feb 19, 2014 at 7:13 AM, Ming Yi <yi02 at="" hotmail.com=""> wrote: > By the way, a small quick question: in the topTable result of limma, the logFC is base 2 or is natural base (2.71828 etc)? I though was base 2, just want to confirm. Take a look at the help file from ?topTable under the Value section, which is where the output is described. Specifically, for the logFC output: """ estimate of the log2-fold-change corresponding to the effect or contrast (for topTableF there may be several columns of log-fold-changes) """ -steve -- Steve Lianoglou Computational Biologist Genentech
ADD REPLY •
0
Entering edit mode
Thx a lot, Steve! I thought I saw this somwhere before. Best Ming > Date: Wed, 19 Feb 2014 09:01:11 -0800 > Subject: Re: [BioC] Problem in limma-voom > From: lianoglou.steve@gene.com > To: yi02@hotmail.com > CC: smyth@wehi.edu.au; bioconductor@r-project.org > > Hi, > > On Wed, Feb 19, 2014 at 7:13 AM, Ming Yi <yi02@hotmail.com> wrote: > > By the way, a small quick question: in the topTable result of limma, the logFC is base 2 or is natural base (2.71828 etc)? I though was base 2, just want to confirm. > > Take a look at the help file from ?topTable under the Value section, > which is where the output is described. > > Specifically, for the logFC output: > > """ > estimate of the log2-fold-change corresponding to the effect or > contrast (for topTableF there may be several columns of > log-fold-changes) > """ > > -steve > > -- > Steve Lianoglou > Computational Biologist > Genentech [[alternative HTML version deleted]]
ADD REPLY •
0
Entering edit mode
Hi, Dr. Gordon and all: I run into some limma design issue not sure what they are: I had three tumor types want to compare with each other: RasP, RasP1Hit, RasPNot here are the critical commands: >targets<-targets[targets$RasType!="RasOnly" & targets$Type=="Tumor",] > show(head(targets)); Subject SampleName Type RasPType RasType 1 TCGA_38_4625 T4625_01A Tumor RasP RasP 5 TCGA_38_4627 T4627_01A Tumor RasP RasP 7 TCGA_38_4632 T4632_01A Tumor RasP RasP 9 TCGA_44_2655 T2655_01A Tumor RasP RasP 11 TCGA_44_2657 T2657_01A Tumor RasP RasP 13 TCGA_44_2661 T2661_01A Tumor RasP RasP > RasTum<-paste(targets$RasType, targets$Type,sep=".") > RasTum<-factor(RasTum, levels=c("RasP1Hit.Tumor","RasP.Tumor","RasPNot.Tumor")); > show(RasTum); [1] RasP.Tumor RasP.Tumor RasP.Tumor RasP.Tumor RasP.Tumor [6] RasP.Tumor RasP.Tumor RasP.Tumor RasP1Hit.Tumor RasP1Hit.Tumor [11] RasP.Tumor RasPNot.Tumor RasP1Hit.Tumor RasP.Tumor RasP1Hit.Tumor [16] RasP1Hit.Tumor RasP.Tumor RasP.Tumor RasPNot.Tumor RasP1Hit.Tumor [21] RasP1Hit.Tumor RasP1Hit.Tumor RasPNot.Tumor RasP.Tumor RasPNot.Tumor [26] RasPNot.Tumor RasPNot.Tumor RasPNot.Tumor RasPNot.Tumor RasPNot.Tumor [31] RasPNot.Tumor RasPNot.Tumor RasPNot.Tumor RasPNot.Tumor RasPNot.Tumor [36] RasPNot.Tumor RasPNot.Tumor RasPNot.Tumor RasPNot.Tumor RasP.Tumor [41] RasP.Tumor RasP.Tumor RasP.Tumor RasPNot.Tumor RasPNot.Tumor Levels: RasP1Hit.Tumor RasP.Tumor RasPNot.Tumor > design<-model.matrix(~RasTum); > colnames(design)<-sub("RasTum","",colnames(design)); > colnames(design)[1] <- "Intercept" > rownames(design)<-targets$Sample > show(design[1:5,]) Intercept RasP.Tumor RasPNot.Tumor T4625_01A 1 1 0 T4627_01A 1 1 0 T4632_01A 1 1 0 T2655_01A 1 1 0 T2657_01A 1 1 0 ... > con.matrix<-makeContrasts( + RasP.Tumor_RasPNot.Tumor=RasP.Tumor- RasPNot.Tumor,RasPRasP1Hit.Tumor_RasPNot.Tumor=RasP.Tumor- RasPNot.Tumor, + RasP1Hit.Tumor_RasPNot.Tumor=-RasPNot.Tumor, + levels=design) > show(con.matrix); Contrasts Levels RasP.Tumor_RasPNot.Tumor RasPRasP1Hit.Tumor_RasPNot.Tumor Intercept 0 0 RasP.Tumor 1 1 RasPNot.Tumor -1 -1 Contrasts Levels RasP1Hit.Tumor_RasPNot.Tumor Intercept 0 RasP.Tumor 0 RasPNot.Tumor -1 My first question is: as shown in show(RasTum), there are three subtype of tumors I want to compare between: Levels: RasP1Hit.Tumor, RasP.Tumor and RasPNot.Tumor shown in show(design[1:5,]), everything is relative to RasP1Hit.Tumor. So in my con.matrix<-makeContrasts command: RasP1Hit.Tumor_RasPNot.Tumor=-RasPNot.Tumor RasP.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor in particular, RasP1Hit.Tumor_RasPNot.Tumor=-RasPNot.Tumor since everything is relative to RasP1Hit.Tumor. However, for contrast: RasPRasP1Hit.Tumor_RasPNot.Tumor, I intended to pool both RasP class and RasP1Hit class together to compare to RasPNot class, I end up with: RasPRasP1Hit.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor the right side is the same as RasP.Tumor_RasPNot.Tumor=RasP.Tumor-RasPNot.Tumor because everything is relative to RasP1Hit.Tumor. So what is the right way to set up contrast RasPRasP1Hit.Tumor_RasPNot.Tumor in makeContrasts command? kind of confused here? Because of the setting, the results are as below: > show(summary(de <- decideTests(fit))); RasP.Tumor_RasPNot.Tumor RasPRasP1Hit.Tumor_RasPNot.Tumor -1 0 0 0 17764 17764 1 0 0 RasP1Hit.Tumor_RasPNot.Tumor -1 0 0 17764 1 0 My 2nd question is: I did a similar design with including all normal samples, since at the same model I want to get tumor vs normal contrasts as well. > con.matrix<-makeContrasts( + RasP.Tumor_RasP.Normal = RasP.Tumor- RasP.Normal,RasPRasP1Hit.Tumor_RasPRasP1Hit.Normal = RasP.Tumor-RasP .Normal-RasP1Hit.Normal, + RasPNot.Tumor_RasPNot.Normal=RasPNot.Tumor-RasPNot.Normal, + RasP.Tumor_RasPNot.Tumor=RasP.Tumor- RasPNot.Tumor,RasPRasP1Hit.Tumor_RasPNot.Tumor=RasP.Tumor- RasPNot.Tumor, + RasP.Normal_RasPNot.Normal=RasP.Normal- RasPNot.Normal,RasPRasP1Hit.Normal_RasPNot.Normal=RasP1Hit.Normal+RasP .Normal-RasPNot.Normal, + Tumor_Normal=RasP.Tumor+ RasPNot.Tumor+-RasP.Normal- RasPNot.Normal, + Tumor_Normal_2=RasP.Tumor+ RasPNot.Tumor-RasP.Normal-RasPNot .Normal-RasP1Hit.Normal, + levels=design) here everything is relative to RasP1Hit.Tumor. and the result as below: > show(summary(de <- decideTests(fit))); RasP.Tumor_RasP.Normal RasPRasP1Hit.Tumor_RasPRasP1Hit.Normal -1 3947 4522 0 9912 8830 1 4274 4781 RasPNot.Tumor_RasPNot.Normal RasP.Tumor_RasPNot.Tumor -1 4809 24 0 8461 18102 1 4863 7 RasPRasP1Hit.Tumor_RasPNot.Tumor RasP.Normal_RasPNot.Normal -1 24 608 0 18102 17158 1 7 367 RasPRasP1Hit.Normal_RasPNot.Normal Tumor_Normal Tumor_Normal_2 -1 2001 5611 5580 0 13934 6728 6651 1 2198 5794 5902 RasP.Tumor_RasPNot.Tumor contrast has small number of DEGs (24+7) at FDR 5% compared to using tumor data alone for the model shown earlier above, this might be cuased by filtering as well and I am fine with that small difference. What is bothering me is RasP.Normal_RasPNot.Normal contrast has about 608+367=~1k DEGs, which is way more than RasP.Tumor_RasPNot.Tumor contrast. I did plot plotMDS and see those normals are mixed together (I did see two subgroups seems having batch effects, but for each batch (if there is), there are RasP.Normal and RasPNot.Normal mixed there), all seem not supposed to having 1k DEGs either. Plus, the tumor contrast RasP.Tumor_RasPNot.Tumor has such small set of DEGs. Biologically, it might be true (population, ethics group, batch effect etc), but I just want to make sure what I did is reasonable. Any advice would be highly appreciated! Thanks a lot for your help! Best Ming > Date: Wed, 19 Feb 2014 17:32:06 +1100 > From: smyth@wehi.EDU.AU > To: yi02@hotmail.com > CC: bioconductor@r-project.org > Subject: Re: Problem in limma-voom > > Hi Ming, > > This has nothing to do with voom. It is the same message that you will > always receive from limma when you try to fit an over-parametrized model. > (By comparison, edgeR would have simply told you the design matrix was not > of full rank and stopped.) > > You can only used a paired analysis (+Patient etc) when the treatments to > be compared are applied to the same patient. You can't compare different > tumour types (KRas and KRasP) which affect different patients. > > Gordon > > > On Wed, 19 Feb 2014, Ming Yi wrote: > > > Hi, Dr. Gordon and all: > > > > > > > > I used limma before, but not much on limma-voom package. I run into some > > issue of warnings and then errors, not sure what is. Below are my > > commands, only show essential parts: > > > > > > > >> show(head(targets)) > > > > ID RasType Subject Type > > 1 T4626_01A KRas TCGA_38_4626 Tumor > > 2 N4626_11A KRas TCGA_38_4626 Normal > > 3 T2668_01A KRas TCGA_44_2668 Tumor > > 4 N2668_11A KRas TCGA_44_2668 Normal > > 5 T6145_01A KRas TCGA_44_6145 Tumor > > > > > RasTum<-paste(targets$RasType, targets$Type,sep=".") > >> RasTum<-factor(RasTum, levels=c("KRas.Tumor","KRas.Normal","KRasP.Tumor","KRasP.Normal")); > >> RasTum > > [1] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor > > [6] KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal > > [11] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor > > [16] KRas.Normal KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal > > [21] KRas.Tumor KRas.Normal KRas.Tumor KRas.Normal KRasP.Tumor > > [26] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal > > [31] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor > > [36] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal > > [41] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor > > [46] KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal > > [51] KRasP.Tumor KRasP.Normal KRasP.Tumor KRasP.Normal KRasP.Tumor > > [56] KRasP.Normal KRasP.Tumor KRasP.Normal > > Levels: KRas.Tumor KRas.Normal KRasP.Tumor KRasP.Normal > > > >> Patient<-factor(targets$Subject); > > > >> rawdata<-rawdata[,c("symbols","entrezID" ,"gene_id",targets$ID)]; > > > >> y<-DGEList(count=rawdata[,4:ncol(rawdata)],genes=rawdata[,1:3]); > > > >> show(dim(y)); > > [1] 20531 58 > > > >> y3<-y[rowSums(y$counts)>=50,]; > > > >> show(dim(y3)); > > [1] 18215 58 > >> y3$samples$lib.size <- colSums(y3$counts) > >> rownames(y3$counts) <- rownames(y3$genes) <- y3$genes$entrezID > >> y3 <- calcNormFactors(y3,method=c("TMM")); > >> show(y3$samples[1:6,]) > > group lib.size norm.factors > > T4626_01A 1 97435132 0.9834587 > > N4626_11A 1 62583672 0.9154837 > > T2668_01A 1 77034746 0.9687860 > > N2668_11A 1 58631311 0.8644352 > > T6145_01A 1 31634419 1.0303445 > > N6145_11A 1 29872761 0.9615021 > >> design<-model.matrix(~Patient+RasTum); > >> show(design[1:6,c(1:10,30:32)]) > > (Intercept) PatientTCGA_38_4626 PatientTCGA_38_4627 PatientTCGA_38_4632 > > 1 1 1 0 0 > > 2 1 1 0 0 > > 3 1 0 0 0 > > 4 1 0 0 0 > > 5 1 0 0 0 > > 6 1 0 0 0 > > PatientTCGA_44_2655 PatientTCGA_44_2657 PatientTCGA_44_2661 > > 1 0 0 0 > > 2 0 0 0 > > 3 0 0 0 > > 4 0 0 0 > > 5 0 0 0 > > 6 0 0 0 > > PatientTCGA_44_2662 PatientTCGA_44_2665 PatientTCGA_44_2668 RasTumKRas.Normal > > 1 0 0 0 0 > > 2 0 0 0 1 > > 3 0 0 1 0 > > 4 0 0 1 1 > > 5 0 0 0 0 > > 6 0 0 0 1 > > RasTumKRasP.Tumor RasTumKRasP.Normal > > 1 0 0 > > 2 0 0 > > 3 0 0 > > 4 0 0 > > 5 0 0 > > 6 0 0 > >> colnames(design)<-sub("Patient","",colnames(design)); > >> colnames(design)<-sub("RasTum","",colnames(design)); > >> colnames(design)[1] <- "Intercept" > >> show(design[1:6,c(1:10,30:32)]) > > Intercept TCGA_38_4626 TCGA_38_4627 TCGA_38_4632 TCGA_44_2655 TCGA_44_2657 > > 1 1 1 0 0 0 0 > > 2 1 1 0 0 0 0 > > 3 1 0 0 0 0 0 > > 4 1 0 0 0 0 0 > > 5 1 0 0 0 0 0 > > 6 1 0 0 0 0 0 > > TCGA_44_2661 TCGA_44_2662 TCGA_44_2665 TCGA_44_2668 KRas.Normal KRasP.Tumor > > 1 0 0 0 0 0 0 > > 2 0 0 0 0 1 0 > > 3 0 0 0 1 0 0 > > 4 0 0 0 1 1 0 > > 5 0 0 0 0 0 0 > > 6 0 0 0 0 1 0 > > KRasP.Normal > > 1 0 > > 2 0 > > 3 0 > > 4 0 > > 5 0 > > 6 0 > >> dim(design) > > [1] 58 32 > >> con.matrix<-makeContrasts( > > + KRas.Tumor_KRas.Normal = -KRas.Normal, > > + KRasP.Tumor_KRasP.Normal = KRasP.Tumor-KRasP.Normal, > > + KRas.Tumor_KRasP.Tumor=-KRasP.Tumor, > > + levels=design) > >> v <- voom(y3,design); > > Coefficients not estimable: KRasP.Normal > > Warning message: > > Partial NA coefficients for 18215 probe(s) > > > >> fit<- lmFit(v,design) > > Coefficients not estimable: KRasP.Normal > > Warning message: > > Partial NA coefficients for 18215 probe(s) > >> fit1<-contrasts.fit(fit, con.matrix) > > Error in contrasts.fit(fit, con.matrix) : > > trying to take contrast of non-estimable coefficient > > > > > > I used a similar way before for array data, but not with limma- voom for RNAseq data. Any suggestions here? > > > > > > > > Thanks a lot in advance! > > > > > > > > Ming > > > > > > > > ATRF/NCI-Frederick > > > > Maryland, USA > > > > > > > > Here is the sessionInfo: > > > >> show(sessionInfo()) > > R version 3.0.1 (2013-05-16) > > Platform: x86_64-unknown-linux-gnu (64-bit) > > > > locale: > > [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C > > [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8 > > [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 > > [7] LC_PAPER=C LC_NAME=C > > [9] LC_ADDRESS=C LC_TELEPHONE=C > > [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C > > > > attached base packages: > > [1] stats graphics grDevices utils datasets methods base > > > > other attached packages: > > [1] edgeR_3.4.2 limma_3.16.8 > > > > > > > > > > > > > > > > > > > > > > ______________________________________________________________________ > The information in this email is confidential and inte...{{dropped:9}}
ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 765 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