Besides being ad-hoc, the codes are almost unreadable. Why you have chosen the cutpoints of 0.5 and 9 ? Regards, Adai On Tue, 2005-09-20 at 10:54 +0200, Morten wrote: > > 2) The result from IN HOUSE data seems to be nice, there are some DE > > genes in the listing. Instead, there is nothing differencially > > expressed genes for the data from OS LAB. So is it POSSIBLE? > > Assuming it is possible 'cause variations which came from operations > > from different lab, what should I do from now on? > > Hi. > > Analyzing dye swap experiments in limma, yeld in my humble experience, very poor > satistics mainly due to gene specific dye bias. In addition when you have > relativily few arrays, you would expect to get poor statistics. > > One rather "ad hoc" method I use (and the statesticians probably dispise the > following. sorry:) > > x <- topTable(fit, adjust="fdr", sort.by="M", number=2000); y <- x[x$P.Value < 1 > & (x$M > 0.5 | x$M < -0.5) & x$A > 9,];y; print("Number of genes in this > list:"); length(y$ID) #which gives you a toplist(y) of genes with M>|0.5| and > A>9 from the limma object "fit" > > this I find usefull to analyse statisticly weak experiemnts (dye swap with few > arrays) to look for interesting genes. > > hope this can be of any help > morten > > _______________________________________________ > Bioconductor mailing list > Bioconductor at stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/bioconductor >

Dear Gordon and other users, Based on Limma usersguide and the hint you sent me, "the #of columns in designmatrix is equal to the #of RNA sources-1 for two-color direct comparisons". This is what I did using command modelMatrix. ## something i don't understand ## 1) why my corfit$consensus is NaN? Am I supposed to get a number which indicates a relation between two block arrays(first 2 and last 2). 2) If ignoring NaN, the result of lmFit has a message "not estimable: KO UV" which is the coefficient I am most interested in ("KO UV/WT control"). What is happening here? Again, my TargetsGenePix is >TargetsGenePix SlideNumber FileName Cy3 Cy5 FG WT rep1 Sep2005 3215L FG WT rep1 Sep2005.gpr WT control WT UV FG WT rep2 Sep2005 3202OL FG WT rep2 Sep2005.gpr WT UV WT control FG KO rep2 Sep2005 3215OL FG KO rep2 Sep2005.gpr KO control KO UV FG KO rep1 Sep2005 3202L FG KO rep1 Sep2005.gpr KO UV KO control design <- modelMatrix(TargetsGenePix, ref="WT control") > design KO control KO UV WT UV FG WT rep1 Sep2005 0 0 1 FG WT rep2 Sep2005 0 0 -1 FG KO rep2 Sep2005 -1 1 0 FG KO rep1 Sep2005 1 -1 0 corfit <- duplicateCorrelation(MAgpr, design, ndups=1, block=c(1,1,2,2)) > corfit$consensus [1] NaN > fit <- lmFit(MAgpr, design) Coefficients not estimable: KO UV Thanks for your time Kevin > > Date: Mon, 19 Sep 2005 15:39:04 -0500 >> From: kevin Lin <khlin at="" odin.mdacc.tmc.edu=""> >> Subject: [BioC] two pair dye-swap (replicates) conducted in different >> labs >> To: bioconductor at stat.math.ethz.ch >> Message-ID: <a06002001bf54d1fd7091@[143.111.133.119]> >> Content-Type: text/plain; charset="us-ascii" ; format="flowed" >> >> Dear BioC, >> >> I have technical replicates are in dye-swap pairs with four hybs >> conducted in house and four hybs conducted in commercial lab. Same >> RNA samples were used. >> >> Please see belows, I did separated analyses using Limma pipeline >> commands. My Qs is >> >> 1) For IN HOUSE analysis, the corfit$consensus is supposed to be a >> negative number according to limma usersguide (the technical >> replicates are dye-swap and should vary in opposite directions). Why >> am I getting a positive number in the dye-swap design? Is there >> something wrong with my design matrix > >Yes. Your experiment has four distinct RNA sources, and at least >two fold different changes to >estimate for each gene. How do you expect to be able to do this >with a design matrix with only >one column? > >Gordon > >> or I misunderstood points in >> the usersguide? >> >> 2) The result from IN HOUSE data seems to be nice, there are some DE >> genes in the listing. Instead, there is nothing differencially >> expressed genes for the data from OS LAB. So is it POSSIBLE? >> Assuming it is possible 'cause variations which came from operations >> from different lab, what should I do from now on? >> >> 3) What do I expect to get if I pull those 8 arrays ( yet, I have not >> done it) and do the analysis together? Even though, two results have >> so much difference in separate analyses (if all this based on correct >> analysis). >> >> I am confused. >> >> I appreciate if anyone could point me out. >> >> Kevin >> >> ***** commands and outputs below ****** >> >> ##################################################################### >> ## IN HOUSE: Two pairs of dye-swap experiments conducted in FG >> ## 1st, 2nd and 3rd, 4th are dye-swap pairs. >> ##################################################################### >> library(limma) >> TargetsGenePix <- readTargets("FGtargets.txt") >> >>> TargetsGenePix >> >> SlideNumber FileName >> Cy3 Cy5 >> FG WT rep1 Sep2005 3215L FG WT rep1 Sep2005.gpr WT > > control WT UV >> FG WT rep2 Sep2005 3202OL FG WT rep2 Sep2005.gpr WT UV >> WT control >> FG KO rep2 Sep2005 3215OL FG KO rep2 Sep2005.gpr KO >> control KO UV >> FG KO rep1 Sep2005 3202L FG KO rep1 Sep2005.gpr KO >> UV KO control >> >> SpotTypes <- readSpotTypes() >> RGgpr <- read.maimages(TargetsGenePix$FileName, >> >>source="genepix",wt.fun=wtflags(w=0),annotation=c("Block","Row","Col umn","Name", >> "controltype")) >> RGgpr$printer <- getLayout(RGgpr$genes) >> RGgpr$genes$Status <- controlStatus(SpotTypes, RGgpr$genes) >> MAgpr <- normalizeWithinArrays(RGgpr,method="loess") >> design <- c(1,-1,1,-1) >> corfit <- duplicateCorrelation(MAgpr, design, ndups=1, block=c(1,1,2,2)) >> >> ## > corfit$consensus >> ## [1] 0.3108369 >> >> fit <- lmFit(MAgpr, design, block=c(1,1,2,2), correlation=corfit$consensus) >> fit <- eBayes(fit) >> top200 <- topTable(fit,n=200,adjust="fdr") >> >> Block Row Column Name controltype Status M A t >> P.Value B >> 7828 41 15 4 793067 false cDNA 1.87 8.24 13.23 >> 0.0341 4.80 >> 2919 16 1 7 1228244 false cDNA -1.89 9.80 -12.42 >> 0.0341 4.45 >> 5178 27 20 2 443884 ignore ignore 1.76 11.55 11.77 >> 0.0341 4.15 >> 2035 11 13 3 1265839 false cDNA -2.52 8.32 -10.92 >> 0.0413 3.73 >> 8538 45 8 2 335555 false cDNA 1.45 10.24 10.10 >> 0.0440 3.27 >> 10978 57 22 2 463982 false cDNA 1.45 13.32 9.78 >> 0.0440 3.08 >> 834 5 7 2 335736 false cDNA 1.46 12.89 9.53 >> 0.0440 2.93 >> 6258 33 11 2 334906 false cDNA -1.34 10.56 -9.45 >> 0.0440 2.87 >> 11112 58 14 8 314112 false cDNA 1.33 11.48 9.24 >> 0.0440 2.74 >> 6651 35 12 3 1264958 false cDNA -1.99 8.27 -9.22 >> 0.0440 2.73 >> 490 3 12 2 334575 false cDNA -1.33 10.15 -9.22 >> 0.0440 2.73 >> 567 3 21 7 480745 ignore ignore -1.48 11.24 -9.07 >> 0.0446 2.63 >> 9662 51 2 6 1111090 false cDNA 1.31 6.09 8.77 >> 0.0508 2.43 >> 6394 34 4 2 1197400 false cDNA -1.28 8.58 -8.66 >> 0.0509 2.35 >> 2407 13 9 7 1197119 false cDNA -1.36 10.43 -8.33 >> 0.0553 2.12 >> 9303 49 7 7 1196439 ignore ignore -1.25 10.91 -8.25 >> 0.0553 2.06 >> 8365 44 10 5 440614 false cDNA 1.34 7.89 8.23 >> 0.0553 2.05 >> 3547 19 8 3 1382084 false cDNA 1.18 8.25 8.21 >> 0.0553 2.03 >> 1929 10 24 1 790765 false cDNA -1.26 7.36 -7.98 >> 0.0581 1.85 >> 9803 51 20 3 1096050 false cDNA 1.25 10.02 7.90 >> 0.0581 1.80 >> >> >>#################################################################### ######### >> ## OS LAB: Two pairs of dye-swap experiments conducted >> >>#################################################################### ######### >> TargetsGenePix <- readTargets("OStargets.txt") >> >>>TargetsGenePix >> >> SlideNumber FileName Cy3 >>Cy5 >> OS WT rep1 Sep2005 2630L OS WT rep1 Sep2005.gpr WT control WT UV >> OS WT rep2 Sep2005 2630OL OS WT rep2 Sep2005.gpr WT UV WT control >> OS KO rep1 Sep2005 2631OL OS KO rep1 Sep2005.gpr KO UV KO control >> OS KO rep2 Sep2005 2631L OS KO rep2 Sep2005.gpr KO control KO UV >> >> snip... >> >> design <- c(1,-1,1,-1) >> corfit <- duplicateCorrelation(MAgpr, design, ndups=1, block=c(1,1,2,2)) >> ## > corfit$consensus.correlation >> ## [1] -0.2315571 >> fit <- lmFit(MAgpr, design, block=c(1,1,2,2), correlation=corfit$consensus) >> fit <- eBayes(fit) >> top200 <- topTable(fit,n=200,adjust="fdr") >> >> Block Row Column Name controltype Status M A t >> P.Value B >> 895 5 14 7 1023591 false cDNA 1.378 5.45 4.67 >> 1 -3.67 >> 10580 55 20 4 750226 false cDNA 1.091 6.03 4.49 >> 1 -3.69 >> 7293 38 20 5 944409 false cDNA 1.039 6.39 4.44 >> 1 -3.70 >> 10117 53 10 5 336470 false cDNA -1.046 8.07 -4.02 > > 1 -3.77 >> 3790 20 14 6 894195 false cDNA -0.786 7.14 -3.90 >> 1 -3.79 >> 10077 53 5 5 718360 false cDNA 2.410 5.63 7.22 >> 1 -3.79 >> 6469 34 13 5 573770 false cDNA 0.982 6.86 3.74 >> 1 -3.81 >> 7476 39 19 4 765454 false cDNA -1.191 6.15 -3.69 >> 1 -3.82 >> 1981 11 6 5 719193 false cDNA -0.850 8.22 -3.69 >> 1 -3.82 >> 6796 36 6 4 905129 false cDNA -0.748 6.25 -3.63 >> 1 -3.84 >> 7015 37 9 7 1197111 false cDNA 0.748 7.90 3.61 >> 1 -3.84 >> 6361 33 24 1 538420 false cDNA 1.122 6.14 3.55 >> 1 -3.85 >> 252 2 7 4 948928 false cDNA -0.759 7.71 -3.54 >> 1 -3.85 >> 9467 50 3 3 1449999 false cDNA 1.321 6.08 5.54 >> 1 -3.85 >> 2053 11 15 5 622182 false cDNA 1.414 6.22 5.16 >> 1 -3.88 >> 10450 55 4 2 634988 false cDNA 0.815 5.97 3.41 >> 1 -3.88 >> 940 5 20 4 762123 false cDNA -0.699 6.28 -3.40 >> 1 -3.88 >> 9575 50 16 7 1020833 false cDNA 0.752 7.69 3.39 >> 1 -3.88 >> 4335 23 10 7 1149905 false cDNA 1.215 6.47 4.98 >> 1 -3.89 >> 3292 17 24 4 777258 false cDNA 0.770 6.25 3.36 >> 1 -3.89 > >_______________________________________________ >Bioconductor mailing list >Bioconductor at stat.math.ethz.ch >https://stat.ethz.ch/mailman/listinfo/bioconductor