Another limma factorial that needs controlling and pairing
2
0
Entering edit mode
sabrina ▴ 10
@sabrina-4210
Last seen 9.9 years ago
Hi, Karl: I have a similar question as yours, did you figure it out since your last post? Other question I have is about using topTable. In your previous design, did you use topTable for each contrast, meaning for each particular contrast, how you find the top list of genes? If that is the case, what is the best way to sort the genes? Did you use B value or adj.p? Or you just used decideTests? Thanks Sabrina
• 1.1k views
ADD COMMENT
0
Entering edit mode
k. brand ▴ 420
@k-brand-1874
Last seen 9.9 years ago
Hi Sabrina, Yes, i believe i did figure out a contrast matrix which allowed me to use limma to indentify DE genes between two factors of interest ('Tissue' and 'Pperiod' in my case) whilst controlling for a third factor ('Time' in my case where the samples were collected at different times of the day as part of another question incorporated in the same experiment but not discussed here). So to be clear, my goal was to find DE genes between all pairwise comparisons of 6 data sets, namely, R.S, R.E, R.L, C.S, C.E and C.L and also for any genes significant for interaction between 'Tissue' and 'Pperiod'. In my original post "[BioC] Another limma factorial that needs controlling and pairing" i'd outlined an approach based closly on the first example discussed in Gorodn Smyth's chapter 23*, pg412. In the end this approach didn't work for me. Instead i succeeded using the "treatment-contrast" parametrization (pg414) example. You'll find a representation of the code i used for my analysis below. I haven't run it, but for the purpose of illustration, together with my targets file should suffice in explaining exactly how i used limma. The most diffciult part for me was understanding the way R sets up this "treatment-contrast" parametrization. To be sure about how i was setting up the contrasts, i manually calculated the Log2FC for each contrast and compared this with the limma output. This confirmation was critical for a dilatante like myself to be sure my costrasts were indeed what i thought they were. I would advise similar such corss checking of your contrasts. As you'll see from my code, i used 'decide tests' since overall significance of genes across all contrasts was of greater help to me in understanding the biology. As to which correction method you use, there is plenty of online info (BioC archives, limma user guide) to clarify this. In short, "global" seemed the best balanced method for my question(s). Final note on controlling for the paring of tissues (being obtained from the same animal): i omitted this for simplicty in the code below. In fact Prof. Smyth clarified out to deal with this using duplicateCorrelation in a later BioC post if this is relevant for you; and i can also repost if needed. hth, Karl *http://www.statsci.org/smyth/pubs/limma-biocbook-reprint.pdf ###not run library(limma) targets <- read.delim("RNA_Targets.txt") Tissue <- factor(targets$Tissue, levels = c("R", "C")) Pperiod <- factor(targets$Pperiod, levels = c("E", "L", "S")) Time <- factor(targets$Time) design <- model.matrix(~Tissue * Pperiod + Time) colnames(design) fit <- lmFit(rma.pp, design) cont.matrix <- cbind( R.E_R.S = c(0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), R.E_R.L = c(0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), R.S_R.L = c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), C.E_C.S = c(0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1), C.E_C.L = c(0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0), C.S_C.L = c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1), R.S_C.S = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1), R.E_C.E = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), R.L_C.L = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0), R.E_R.S__C.E_C.S = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1),#int R.E_R.L__C.E_C.L = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0))#int.2 fit2 <- contrasts.fit(fit, cont.matrix) fit3 <- eBayes(fit2) #write resutls of contrasts for 0.05 results <- (decideTests(fit3, method="global", adjust.method="BH", p.value=0.05)) write.fit(fit3, results, file="w.int.glob.BH.0.05.csv", sep=",") #write resutls of contrasts for 0.01 results <- (decideTests(fit3, method="global", adjust.method="BH", p.value=0.01)) write.fit(fit3, results, digits=10, file="w.int.glob.BH.0.01.csv", sep=",") ###end "RNA_Targets.txt": FileName Tissue Pperiod Time Animal 01file.CEL R S 1 1 02file.CEL C S 1 1 03file.CEL R S 2 2 04file.CEL C S 2 2 05file.CEL R S 3 3 06file.CEL C S 3 3 07file.CEL R S 4 4 08file.CEL C S 4 4 09file.CEL R S 5 5 10file.CEL C S 5 5 11file.CEL R S 6 6 12file.CEL C S 6 6 13file.CEL R S 7 7 14file.CEL C S 7 7 15file.CEL R S 8 8 16file.CEL C S 8 8 17file.CEL R S 9 9 18file.CEL C S 9 9 19file.CEL R S 10 10 20file.CEL C S 10 10 21file.CEL R S 11 11 22file.CEL C S 11 11 23file.CEL R S 12 12 24file.CEL C S 12 12 25file.CEL R S 13 13 26file.CEL C S 13 13 27file.CEL R S 14 14 28file.CEL C S 14 14 29file.CEL R S 15 15 30file.CEL C S 15 15 31file.CEL R S 16 16 32file.CEL C S 16 16 33file.CEL R E 1 17 34file.CEL C E 1 17 35file.CEL R E 2 18 36file.CEL C E 2 18 37file.CEL R E 3 19 38file.CEL C E 3 19 39file.CEL R E 4 20 40file.CEL C E 4 20 41file.CEL R E 5 21 42file.CEL C E 5 21 43file.CEL R E 6 22 44file.CEL C E 6 22 45file.CEL R E 7 23 46file.CEL C E 7 23 47file.CEL R E 8 24 48file.CEL C E 8 24 49file.CEL R E 9 25 50file.CEL C E 9 25 51file.CEL R E 10 26 52file.CEL C E 10 26 53file.CEL R E 11 27 54file.CEL C E 11 27 55file.CEL R E 12 28 56file.CEL C E 12 28 57file.CEL R E 13 29 58file.CEL C E 13 29 59file.CEL R E 14 30 60file.CEL C E 14 30 61file.CEL R E 15 31 62file.CEL C E 15 31 63file.CEL R E 16 32 64file.CEL C E 16 32 65file.CEL R L 1 33 66file.CEL C L 1 33 67file.CEL R L 2 34 68file.CEL C L 2 34 69file.CEL R L 3 35 70file.CEL C L 3 35 71file.CEL R L 4 36 72file.CEL C L 4 36 73file.CEL R L 5 37 74file.CEL C L 5 37 75file.CEL R L 6 38 76file.CEL C L 6 38 77file.CEL R L 7 39 78file.CEL C L 7 39 79file.CEL R L 8 40 80file.CEL C L 8 40 81file.CEL R L 9 41 82file.CEL C L 9 41 83file.CEL R L 10 42 84file.CEL C L 10 42 85file.CEL R L 11 43 86file.CEL C L 11 43 87file.CEL R L 12 44 88file.CEL C L 12 44 89file.CEL R L 13 45 90file.CEL C L 13 45 91file.CEL R L 14 46 92file.CEL C L 14 46 93file.CEL R L 15 47 94file.CEL C L 15 47 95file.CEL R L 16 48 96file.CEL C L 16 48 On 8/10/2010 4:21 PM, sabrina wrote: > > > > Hi, Karl: > I have a similar question as yours, did you figure it out since your last post? > > Other question I have is about using topTable. In your previous > design, did you use topTable for each contrast, meaning for each particular > contrast, how you find the top list of genes? If that is the case, what is the > best way to sort the genes? Did you use B value or adj.p? Or you just used > decideTests? > > Thanks > > Sabrina > > _______________________________________________ > Bioconductor mailing list > Bioconductor at stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor -- Karl Brand <k.brand at="" erasmusmc.nl=""> Department of Genetics Erasmus MC Dr Molewaterplein 50 3015 GE Rotterdam P +31 (0)10 704 3409 | F +31 (0)10 704 4743 | M +31 (0)642 777 268
ADD COMMENT
0
Entering edit mode
Hi, Karl: Thank you for the detailed explanation. I need to understand the Treatment-contrast parametrization since I always used the design matrix as what you did the first time. I need to digest more about it. Just a quick comment. I noticed that you used decideTests(data, methods="global",adj="BH"...) In Limma's user's guide $10.3, it mentioned that there is no theory prove about the combination of BH and global will correctly control the FDR for negatively correlated contrasts (though simulation is ok), I wonder if you have tried other ways, such as separate topTable etc. Does it make difference? Thanks! For some reason, I posted the same question as you suggested to the mailinglist, it never got posted, I hope this time it works,,, Sabrina On Wed, Aug 11, 2010 at 9:59 AM, Karl Brand <k.brand at="" erasmusmc.nl=""> wrote: > Hi Sabrina, > > Yes, i believe i did figure out a contrast matrix which allowed me to use > limma to indentify DE genes between two factors of interest ('Tissue' and > 'Pperiod' in my case) whilst controlling for a third factor ('Time' in my > case where the samples were collected at different times of the day as part > of another question incorporated in the same experiment but not discussed > here). So to be clear, my goal was to find DE genes between all pairwise > comparisons of 6 data sets, namely, R.S, R.E, R.L, C.S, C.E and C.L and also > for any genes significant for interaction between 'Tissue' and 'Pperiod'. > > In my original post "[BioC] Another limma factorial that needs controlling > and ?pairing" i'd outlined an approach based closly on the first example > discussed in Gorodn Smyth's chapter 23*, pg412. In the end this approach > didn't work for me. Instead i succeeded using the "treatment- contrast" > parametrization (pg414) example. You'll find a representation of the code i > used for my analysis below. I haven't run it, but for the purpose of > illustration, together with my targets file should suffice in explaining > exactly how i used limma. The most diffciult part for me was understanding > the way R sets up this "treatment-contrast" parametrization. To be sure > about how i was setting up the contrasts, i manually calculated the Log2FC > for each contrast and compared this with the limma output. This confirmation > was critical for a dilatante like myself to be sure my costrasts were indeed > what i thought they were. I would advise similar such corss checking of your > contrasts. > > As you'll see from my code, i used 'decide tests' since overall significance > of genes across all contrasts was of greater help to me in understanding the > biology. As to which correction method you use, there is plenty of online > info (BioC archives, limma user guide) to clarify this. In short, "global" > seemed the best balanced method for my question(s). > > Final note on controlling for the paring of tissues (being obtained from the > same animal): i omitted this for simplicty in the code below. In fact Prof. > Smyth clarified out to deal with this using duplicateCorrelation in a later > BioC post if this is relevant for you; and i can also repost if needed. > > hth, > > Karl > > *http://www.statsci.org/smyth/pubs/limma-biocbook-reprint.pdf > > ###not run > library(limma) > targets <- read.delim("RNA_Targets.txt") > Tissue <- factor(targets$Tissue, levels = c("R", "C")) > Pperiod <- factor(targets$Pperiod, levels = c("E", "L", "S")) > Time <- factor(targets$Time) > > design <- model.matrix(~Tissue * Pperiod + Time) > colnames(design) > > fit <- lmFit(rma.pp, design) > cont.matrix <- cbind( > > R.E_R.S = c(0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), > R.E_R.L = c(0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), > R.S_R.L = c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), > > C.E_C.S = c(0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1), > C.E_C.L = c(0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0), > C.S_C.L = c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1), > > R.S_C.S = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1), > R.E_C.E = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), > R.L_C.L = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0), > > R.E_R.S__C.E_C.S = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1),#int > R.E_R.L__C.E_C.L = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0))#int.2 > > fit2 <- contrasts.fit(fit, cont.matrix) > fit3 <- eBayes(fit2) > > #write resutls of contrasts for 0.05 > results <- (decideTests(fit3, method="global", adjust.method="BH", > p.value=0.05)) > write.fit(fit3, results, file="w.int.glob.BH.0.05.csv", sep=",") > > #write resutls of contrasts for 0.01 > results <- (decideTests(fit3, method="global", adjust.method="BH", > p.value=0.01)) > write.fit(fit3, results, digits=10, file="w.int.glob.BH.0.01.csv", sep=",") > > ###end > > "RNA_Targets.txt": > > FileName ? ?Tissue ? ?Pperiod ? ?Time ? ?Animal > 01file.CEL ? ?R ? ?S ? ?1 ? ?1 > 02file.CEL ? ?C ? ?S ? ?1 ? ?1 > 03file.CEL ? ?R ? ?S ? ?2 ? ?2 > 04file.CEL ? ?C ? ?S ? ?2 ? ?2 > 05file.CEL ? ?R ? ?S ? ?3 ? ?3 > 06file.CEL ? ?C ? ?S ? ?3 ? ?3 > 07file.CEL ? ?R ? ?S ? ?4 ? ?4 > 08file.CEL ? ?C ? ?S ? ?4 ? ?4 > 09file.CEL ? ?R ? ?S ? ?5 ? ?5 > 10file.CEL ? ?C ? ?S ? ?5 ? ?5 > 11file.CEL ? ?R ? ?S ? ?6 ? ?6 > 12file.CEL ? ?C ? ?S ? ?6 ? ?6 > 13file.CEL ? ?R ? ?S ? ?7 ? ?7 > 14file.CEL ? ?C ? ?S ? ?7 ? ?7 > 15file.CEL ? ?R ? ?S ? ?8 ? ?8 > 16file.CEL ? ?C ? ?S ? ?8 ? ?8 > 17file.CEL ? ?R ? ?S ? ?9 ? ?9 > 18file.CEL ? ?C ? ?S ? ?9 ? ?9 > 19file.CEL ? ?R ? ?S ? ?10 ? ?10 > 20file.CEL ? ?C ? ?S ? ?10 ? ?10 > 21file.CEL ? ?R ? ?S ? ?11 ? ?11 > 22file.CEL ? ?C ? ?S ? ?11 ? ?11 > 23file.CEL ? ?R ? ?S ? ?12 ? ?12 > 24file.CEL ? ?C ? ?S ? ?12 ? ?12 > 25file.CEL ? ?R ? ?S ? ?13 ? ?13 > 26file.CEL ? ?C ? ?S ? ?13 ? ?13 > 27file.CEL ? ?R ? ?S ? ?14 ? ?14 > 28file.CEL ? ?C ? ?S ? ?14 ? ?14 > 29file.CEL ? ?R ? ?S ? ?15 ? ?15 > 30file.CEL ? ?C ? ?S ? ?15 ? ?15 > 31file.CEL ? ?R ? ?S ? ?16 ? ?16 > 32file.CEL ? ?C ? ?S ? ?16 ? ?16 > 33file.CEL ? ?R ? ?E ? ?1 ? ?17 > 34file.CEL ? ?C ? ?E ? ?1 ? ?17 > 35file.CEL ? ?R ? ?E ? ?2 ? ?18 > 36file.CEL ? ?C ? ?E ? ?2 ? ?18 > 37file.CEL ? ?R ? ?E ? ?3 ? ?19 > 38file.CEL ? ?C ? ?E ? ?3 ? ?19 > 39file.CEL ? ?R ? ?E ? ?4 ? ?20 > 40file.CEL ? ?C ? ?E ? ?4 ? ?20 > 41file.CEL ? ?R ? ?E ? ?5 ? ?21 > 42file.CEL ? ?C ? ?E ? ?5 ? ?21 > 43file.CEL ? ?R ? ?E ? ?6 ? ?22 > 44file.CEL ? ?C ? ?E ? ?6 ? ?22 > 45file.CEL ? ?R ? ?E ? ?7 ? ?23 > 46file.CEL ? ?C ? ?E ? ?7 ? ?23 > 47file.CEL ? ?R ? ?E ? ?8 ? ?24 > 48file.CEL ? ?C ? ?E ? ?8 ? ?24 > 49file.CEL ? ?R ? ?E ? ?9 ? ?25 > 50file.CEL ? ?C ? ?E ? ?9 ? ?25 > 51file.CEL ? ?R ? ?E ? ?10 ? ?26 > 52file.CEL ? ?C ? ?E ? ?10 ? ?26 > 53file.CEL ? ?R ? ?E ? ?11 ? ?27 > 54file.CEL ? ?C ? ?E ? ?11 ? ?27 > 55file.CEL ? ?R ? ?E ? ?12 ? ?28 > 56file.CEL ? ?C ? ?E ? ?12 ? ?28 > 57file.CEL ? ?R ? ?E ? ?13 ? ?29 > 58file.CEL ? ?C ? ?E ? ?13 ? ?29 > 59file.CEL ? ?R ? ?E ? ?14 ? ?30 > 60file.CEL ? ?C ? ?E ? ?14 ? ?30 > 61file.CEL ? ?R ? ?E ? ?15 ? ?31 > 62file.CEL ? ?C ? ?E ? ?15 ? ?31 > 63file.CEL ? ?R ? ?E ? ?16 ? ?32 > 64file.CEL ? ?C ? ?E ? ?16 ? ?32 > 65file.CEL ? ?R ? ?L ? ?1 ? ?33 > 66file.CEL ? ?C ? ?L ? ?1 ? ?33 > 67file.CEL ? ?R ? ?L ? ?2 ? ?34 > 68file.CEL ? ?C ? ?L ? ?2 ? ?34 > 69file.CEL ? ?R ? ?L ? ?3 ? ?35 > 70file.CEL ? ?C ? ?L ? ?3 ? ?35 > 71file.CEL ? ?R ? ?L ? ?4 ? ?36 > 72file.CEL ? ?C ? ?L ? ?4 ? ?36 > 73file.CEL ? ?R ? ?L ? ?5 ? ?37 > 74file.CEL ? ?C ? ?L ? ?5 ? ?37 > 75file.CEL ? ?R ? ?L ? ?6 ? ?38 > 76file.CEL ? ?C ? ?L ? ?6 ? ?38 > 77file.CEL ? ?R ? ?L ? ?7 ? ?39 > 78file.CEL ? ?C ? ?L ? ?7 ? ?39 > 79file.CEL ? ?R ? ?L ? ?8 ? ?40 > 80file.CEL ? ?C ? ?L ? ?8 ? ?40 > 81file.CEL ? ?R ? ?L ? ?9 ? ?41 > 82file.CEL ? ?C ? ?L ? ?9 ? ?41 > 83file.CEL ? ?R ? ?L ? ?10 ? ?42 > 84file.CEL ? ?C ? ?L ? ?10 ? ?42 > 85file.CEL ? ?R ? ?L ? ?11 ? ?43 > 86file.CEL ? ?C ? ?L ? ?11 ? ?43 > 87file.CEL ? ?R ? ?L ? ?12 ? ?44 > 88file.CEL ? ?C ? ?L ? ?12 ? ?44 > 89file.CEL ? ?R ? ?L ? ?13 ? ?45 > 90file.CEL ? ?C ? ?L ? ?13 ? ?45 > 91file.CEL ? ?R ? ?L ? ?14 ? ?46 > 92file.CEL ? ?C ? ?L ? ?14 ? ?46 > 93file.CEL ? ?R ? ?L ? ?15 ? ?47 > 94file.CEL ? ?C ? ?L ? ?15 ? ?47 > 95file.CEL ? ?R ? ?L ? ?16 ? ?48 > 96file.CEL ? ?C ? ?L ? ?16 ? ?48 > > > > On 8/10/2010 4:21 PM, sabrina wrote: >> >> >> >> Hi, Karl: >> I have a similar question as yours, did you figure it out since your last >> post? >> >> Other question I have is about using topTable. In your previous >> design, did you use topTable for each contrast, meaning for each >> particular >> contrast, how you find the top list of genes? If that is the case, what is >> the >> best way to sort the genes? Did you use B value or adj.p? Or you just used >> decideTests? >> >> Thanks >> >> Sabrina >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at stat.math.ethz.ch >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: >> http://news.gmane.org/gmane.science.biology.informatics.conductor > > -- > Karl Brand <k.brand at="" erasmusmc.nl=""> > Department of Genetics > Erasmus MC > Dr Molewaterplein 50 > 3015 GE Rotterdam > P +31 (0)10 704 3409 | F +31 (0)10 704 4743 | M +31 (0)642 777 268 >
ADD REPLY
0
Entering edit mode
k. brand ▴ 420
@k-brand-1874
Last seen 9.9 years ago
Hi Sabrina, On 8/16/2010 3:33 PM, Sabrina AT Case wrote: > Hi, Karl: > Thank you for the detailed explanation. I need to understand the > Treatment-contrast parametrization since I always used the design > matrix as what you did the first time. I need to digest more about > it. > Lacking statistical training, i devoted considerable time to both the theory and practice of this method of performing multifactorial ANOVA. Getting even just a little face to face time with a statistican can also be invaluable. Failing this i suggest searching/posting the general R-help forum for gaining theoretical and practical advice on the subject. > Just a quick comment. I noticed that you used > decideTests(data, methods="global",adj="BH"...) > In Limma's user's guide $10.3, it mentioned that there is no theory > prove about the combination of BH and global will correctly control > the FDR for negatively correlated contrasts (though simulation is ok), > I wonder if you have tried other ways, such as separate topTable etc. > Does it make difference? Thanks! Yes, i also tried method="nestedF" and "hierarchical" and "separate" (analagous to usin topTable() for each of the contrasts) out of interest and can confirm the results btwn methods were moderatley different with my data, but didn't change the main 'message' emerging from my gene lists. Irrespective, I was careful in choosing method="global" as everything described in the user guide (and discussed in the Bio-C forum) indicated that this approach was most appropriate for my questions. hth, karl > > For some reason, I posted the same question as you suggested to the > mailinglist, it never got posted, I hope this time it works,,, > > Sabrina > > On Wed, Aug 11, 2010 at 9:59 AM, Karl Brand<k.brand at="" erasmusmc.nl=""> wrote: >> Hi Sabrina, >> >> Yes, i believe i did figure out a contrast matrix which allowed me to use >> limma to indentify DE genes between two factors of interest ('Tissue' and >> 'Pperiod' in my case) whilst controlling for a third factor ('Time' in my >> case where the samples were collected at different times of the day as part >> of another question incorporated in the same experiment but not discussed >> here). So to be clear, my goal was to find DE genes between all pairwise >> comparisons of 6 data sets, namely, R.S, R.E, R.L, C.S, C.E and C.L and also >> for any genes significant for interaction between 'Tissue' and 'Pperiod'. >> >> In my original post "[BioC] Another limma factorial that needs controlling >> and pairing" i'd outlined an approach based closly on the first example >> discussed in Gorodn Smyth's chapter 23*, pg412. In the end this approach >> didn't work for me. Instead i succeeded using the "treatment- contrast" >> parametrization (pg414) example. You'll find a representation of the code i >> used for my analysis below. I haven't run it, but for the purpose of >> illustration, together with my targets file should suffice in explaining >> exactly how i used limma. The most diffciult part for me was understanding >> the way R sets up this "treatment-contrast" parametrization. To be sure >> about how i was setting up the contrasts, i manually calculated the Log2FC >> for each contrast and compared this with the limma output. This confirmation >> was critical for a dilatante like myself to be sure my costrasts were indeed >> what i thought they were. I would advise similar such corss checking of your >> contrasts. >> >> As you'll see from my code, i used 'decide tests' since overall significance >> of genes across all contrasts was of greater help to me in understanding the >> biology. As to which correction method you use, there is plenty of online >> info (BioC archives, limma user guide) to clarify this. In short, "global" >> seemed the best balanced method for my question(s). >> >> Final note on controlling for the paring of tissues (being obtained from the >> same animal): i omitted this for simplicty in the code below. In fact Prof. >> Smyth clarified out to deal with this using duplicateCorrelation in a later >> BioC post if this is relevant for you; and i can also repost if needed. >> >> hth, >> >> Karl >> >> *http://www.statsci.org/smyth/pubs/limma-biocbook-reprint.pdf >> >> ###not run >> library(limma) >> targets<- read.delim("RNA_Targets.txt") >> Tissue<- factor(targets$Tissue, levels = c("R", "C")) >> Pperiod<- factor(targets$Pperiod, levels = c("E", "L", "S")) >> Time<- factor(targets$Time) >> >> design<- model.matrix(~Tissue * Pperiod + Time) >> colnames(design) >> >> fit<- lmFit(rma.pp, design) >> cont.matrix<- cbind( >> >> R.E_R.S = c(0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), >> R.E_R.L = c(0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), >> R.S_R.L = c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), >> >> C.E_C.S = c(0,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1), >> C.E_C.L = c(0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0), >> C.S_C.L = c(0,0,-1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,1), >> >> R.S_C.S = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1), >> R.E_C.E = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), >> R.L_C.L = c(0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0), >> >> R.E_R.S__C.E_C.S = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1),#int >> R.E_R.L__C.E_C.L = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0))#int.2 >> >> fit2<- contrasts.fit(fit, cont.matrix) >> fit3<- eBayes(fit2) >> >> #write resutls of contrasts for 0.05 >> results<- (decideTests(fit3, method="global", adjust.method="BH", >> p.value=0.05)) >> write.fit(fit3, results, file="w.int.glob.BH.0.05.csv", sep=",") >> >> #write resutls of contrasts for 0.01 >> results<- (decideTests(fit3, method="global", adjust.method="BH", >> p.value=0.01)) >> write.fit(fit3, results, digits=10, file="w.int.glob.BH.0.01.csv", sep=",") >> >> ###end >> >> "RNA_Targets.txt": >> >> FileName Tissue Pperiod Time Animal >> 01file.CEL R S 1 1 >> 02file.CEL C S 1 1 >> 03file.CEL R S 2 2 >> 04file.CEL C S 2 2 >> 05file.CEL R S 3 3 >> 06file.CEL C S 3 3 >> 07file.CEL R S 4 4 >> 08file.CEL C S 4 4 >> 09file.CEL R S 5 5 >> 10file.CEL C S 5 5 >> 11file.CEL R S 6 6 >> 12file.CEL C S 6 6 >> 13file.CEL R S 7 7 >> 14file.CEL C S 7 7 >> 15file.CEL R S 8 8 >> 16file.CEL C S 8 8 >> 17file.CEL R S 9 9 >> 18file.CEL C S 9 9 >> 19file.CEL R S 10 10 >> 20file.CEL C S 10 10 >> 21file.CEL R S 11 11 >> 22file.CEL C S 11 11 >> 23file.CEL R S 12 12 >> 24file.CEL C S 12 12 >> 25file.CEL R S 13 13 >> 26file.CEL C S 13 13 >> 27file.CEL R S 14 14 >> 28file.CEL C S 14 14 >> 29file.CEL R S 15 15 >> 30file.CEL C S 15 15 >> 31file.CEL R S 16 16 >> 32file.CEL C S 16 16 >> 33file.CEL R E 1 17 >> 34file.CEL C E 1 17 >> 35file.CEL R E 2 18 >> 36file.CEL C E 2 18 >> 37file.CEL R E 3 19 >> 38file.CEL C E 3 19 >> 39file.CEL R E 4 20 >> 40file.CEL C E 4 20 >> 41file.CEL R E 5 21 >> 42file.CEL C E 5 21 >> 43file.CEL R E 6 22 >> 44file.CEL C E 6 22 >> 45file.CEL R E 7 23 >> 46file.CEL C E 7 23 >> 47file.CEL R E 8 24 >> 48file.CEL C E 8 24 >> 49file.CEL R E 9 25 >> 50file.CEL C E 9 25 >> 51file.CEL R E 10 26 >> 52file.CEL C E 10 26 >> 53file.CEL R E 11 27 >> 54file.CEL C E 11 27 >> 55file.CEL R E 12 28 >> 56file.CEL C E 12 28 >> 57file.CEL R E 13 29 >> 58file.CEL C E 13 29 >> 59file.CEL R E 14 30 >> 60file.CEL C E 14 30 >> 61file.CEL R E 15 31 >> 62file.CEL C E 15 31 >> 63file.CEL R E 16 32 >> 64file.CEL C E 16 32 >> 65file.CEL R L 1 33 >> 66file.CEL C L 1 33 >> 67file.CEL R L 2 34 >> 68file.CEL C L 2 34 >> 69file.CEL R L 3 35 >> 70file.CEL C L 3 35 >> 71file.CEL R L 4 36 >> 72file.CEL C L 4 36 >> 73file.CEL R L 5 37 >> 74file.CEL C L 5 37 >> 75file.CEL R L 6 38 >> 76file.CEL C L 6 38 >> 77file.CEL R L 7 39 >> 78file.CEL C L 7 39 >> 79file.CEL R L 8 40 >> 80file.CEL C L 8 40 >> 81file.CEL R L 9 41 >> 82file.CEL C L 9 41 >> 83file.CEL R L 10 42 >> 84file.CEL C L 10 42 >> 85file.CEL R L 11 43 >> 86file.CEL C L 11 43 >> 87file.CEL R L 12 44 >> 88file.CEL C L 12 44 >> 89file.CEL R L 13 45 >> 90file.CEL C L 13 45 >> 91file.CEL R L 14 46 >> 92file.CEL C L 14 46 >> 93file.CEL R L 15 47 >> 94file.CEL C L 15 47 >> 95file.CEL R L 16 48 >> 96file.CEL C L 16 48 >> >> >> >> On 8/10/2010 4:21 PM, sabrina wrote: >>> >>> >>> >>> Hi, Karl: >>> I have a similar question as yours, did you figure it out since your last >>> post? >>> >>> Other question I have is about using topTable. In your previous >>> design, did you use topTable for each contrast, meaning for each >>> particular >>> contrast, how you find the top list of genes? If that is the case, what is >>> the >>> best way to sort the genes? Did you use B value or adj.p? Or you just used >>> decideTests? >>> >>> Thanks >>> >>> Sabrina >>> >>> _______________________________________________ >>> Bioconductor mailing list >>> Bioconductor at stat.math.ethz.ch >>> https://stat.ethz.ch/mailman/listinfo/bioconductor >>> Search the archives: >>> http://news.gmane.org/gmane.science.biology.informatics.conductor >> >> -- >> Karl Brand<k.brand at="" erasmusmc.nl=""> >> Department of Genetics >> Erasmus MC >> Dr Molewaterplein 50 >> 3015 GE Rotterdam >> P +31 (0)10 704 3409 | F +31 (0)10 704 4743 | M +31 (0)642 777 268 >> -- Karl Brand Department of Genetics Erasmus MC Dr Molewaterplein 50 3015 GE Rotterdam T +31 (0)10 704 3457 | F +31 (0)10 704 4743 | M +31 (0)642 777 268
ADD COMMENT

Login before adding your answer.

Traffic: 469 users visited in the last hour
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