limma single-channel analysis and intraspotCorrelation
1
0
Entering edit mode
Ren Na ▴ 250
@ren-na-870
Last seen 10.3 years ago
Hi, I have a data set of four slides (two pairs of dye swap), slide# cy3 cy5 slide1 wo wr slide2 wr wo slide3 mo mr slide4 mr mo wo: wildtype without treatment wr: wildtype with treatment mo: mutant without treatment mr: mutant with treatment I have been thinking if it is possible to get lists of differentially expressed genes between wr and wo, genes between mr and mo, genes between mr and wr, genes between mo and wo, and genes responding to treatment differently in mutant compared to wild-type. I know factorial design is the better way to go. But based on the data I have now, can I use single channel analysis to get comparison of interest? What I tried and what I will do are like the following, library(limma) targets2 <- readTargets("Targets2.txt") RG<-read.maimages(targets2$FileName, source="spot",wt.fun=wtarea(100)) RG$genes<- readGAL("Mouse24052004_635final2.txt") RG$printer<-getLayout(RG$genes) spottypes<-readSpotTypes("spottypes2.txt") RG$genes$Status<- controlStatus(spottypes, RG$genes) # assign weight 0 to missing spots w<-modifyWeights(RG$weights,status=RG$genes$Status, "miss",0) RG$weights<-w RG<-backgroundCorrect(RG,method="minimum") MA<-normalizeWithinArrays(RG) MA<-normalizeBetweenArrays(MA, method="quantile") targets2.sc <- targetsA2C(targets2) design.sc<-model.matrix(~0+factortargets2.sc$Target)+factor(targets2. sc$channel)) colnamesdesign.sc)<-c("mo","mr","wo","wr","ch") corfit<-intraspotCorrelation(MA,design.sc) # got error from intraspotCorrelation() # what I will do fit<-lmscFit(MA,design.sc,correlation=corfit$consensus) cont.matrix<-makeContrasts(mrvsmo=mr-mo,wrvswo=wr-wo,mrvswr=mr- wr,movswo=mo-wo,diff=(mr-mo)-(wr-wo),levels=design.sc) fit2<-contrast.fit(fit,cont.matrix) fit2<-eBayes(fit2) Then use topTable() to get lists of differentially expressed genes for each of five contrasts. Here are my questions, 1) Is it reasonable? 2) When I run to function intraspotCorrelation, I got the error, > corfit<-intraspotCorrelation(MA,design.sc) Loading required package: statmod Attaching package 'statmod': The following object(s) are masked from package:limma : matvec vecmat Error in chol(ZVZ + lambda * I) : the leading minor of order 1 is not positive definite In addition: Warning message: reml: Max iterations exceeded in: remlscore(y, X, Z) > traceback() 3: chol(ZVZ + lambda * I) 2: remlscore(y, X, Z) 1: intraspotCorrelation(MA, design.sc) > What could cause this error? 3) If I can not get lists of differentially expressed genes between mr and wr, genes between mo and wo, and the genes which respond to treatment differently in mutant compared to wildtype, I still like to get differential expression between mr and mo by doing the following, design<-c(1,-1) fit<-lmFit(MA, design) fit<-eBayes(fit) topTable(fit,adjust="fdr") and also I can do the same thing to get differential expression between wr and wo. But my concern is that if the number of replicates is too small and the result won't be reliable? Thanks for any help and Thank Dr. Gordon Smyth for new limma user's guide Ren [[alternative HTML version deleted]]
GO limma ASSIGN GO limma ASSIGN • 1.3k views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 1 hour ago
WEHI, Melbourne, Australia
>Wed Nov 3 17:30:23 CET 2004 >Hi, >I have a data set of four slides (two pairs of dye swap), >slide# cy3 cy5 >slide1 wo wr >slide2 wr wo >slide3 mo mr >slide4 mr mo > >wo: wildtype without treatment >wr: wildtype with treatment >mo: mutant without treatment >mr: mutant with treatment > >I have been thinking if it is possible to get lists of differentially >expressed genes between wr and wo, genes between mr and mo, genes between >mr and wr, genes between mo and wo, and genes responding to treatment >differently in mutant compared to wild-type. I know factorial design is >the better way to go. But based on the data I have now, can I use single >channel analysis to get comparison of interest? What I tried and what I >will do are like the following, In principle you can, but I think your dataset is too small for what you're trying to achieve. You're trying to estimate 5 linear model coefficients and (implicitly) 2 variance parameters, making a total 7 parameters per gene with only 8 single-channel observations, some of which may be missing. If you remove the dye effect from your model and subset MA so that it contains only those genes which are non-missing and not negative controls, then the fit will probably work. Gordon >library(limma) >targets2 <- readTargets("Targets2.txt") >RG<-read.maimages(targets2$FileName, source="spot",wt.fun=wtarea(100)) >RG$genes<- readGAL("Mouse24052004_635final2.txt") >RG$printer<-getLayout(RG$genes) >spottypes<-readSpotTypes("spottypes2.txt") >RG$genes$Status<- controlStatus(spottypes, RG$genes) ># assign weight 0 to missing spots >w<-modifyWeights(RG$weights,status=RG$genes$Status, "miss",0) >RG$weights<-w >RG<-backgroundCorrect(RG,method="minimum") >MA<-normalizeWithinArrays(RG) >MA<-normalizeBetweenArrays(MA, method="quantile") >targets2.sc <- targetsA2C(targets2) >design.sc<-model.matrix(~0+factortargets2.sc$Target)+factor(targets2 .sc$channel)) >colnamesdesign.sc)<-c("mo","mr","wo","wr","ch") >corfit<-intraspotCorrelation(MA,design.sc) > ># got error from intraspotCorrelation()
ADD COMMENT
0
Entering edit mode
Dear Users, I am following the example on Lab 5: Cluster analysis (June 2003) with my own data. I have filtered my expression set as shown on the example and I get the following > sub <- genefilter(X,ffun) > sum(sub) [1] 1124 I save this subset of genes and then log transform it. But when I type the next command I get the following error: > X <- X[sub,] > X <- log2(X) > RawDataSub <- Raw.Data[,sub] Error in Raw.Data[, sub] : (subscript) logical subscript too long Why do I get this error?? Also, if I have stored the subset expression data as X, why is Raw.Data [,sub] using [,sub] again? I don?t really understand this step, if anyone could explain its purpose. I?m running R 1.9.1 on an XP computer Thanks a lot for your help David
ADD REPLY

Login before adding your answer.

Traffic: 558 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