Hi all, I was just wanted someone to have a look at my rcode to see if it is correct. This analysis (Illumina) arrays is comparing CAN vs No-CAN via a common reference. I have 8 pts with CAN, 13 pts without CAN, 3 ref along with 2 with fibrosis only (FO) and 2 with no histology (NB). The contrast i wanted to make is CAN vs No-CAN via CAN vs Ref and NoCAN vs Ref. Ignoring the FO and NB. I have put this code together from the Limma and beadarry guides and general internet published documents. I get a diff exp list at the end, I just wanted to check to see if the analysis is correct, so I can trust the final diff exp list. Sorry, if the code i messy, I have taught myself to use Limma and R. Thank you, Matt Rcode: *library(limma)* *library(beadarray)* <>*setwd("d:/Illumina_results/expset")* <>* BSData = readBeadSummaryData("expset-QCgood-4.txt", qcFile=NULL, sampleSheet=NULL,sep="\t",skip=8, ProbeID="PROBE_ID", columns = list(exprs = "AVG_Signal", se.exprs="BEAD_STDERR", NoBeads = "Avg_NBEADS", Detection="Detection Pval"), dec=".")* <>* BSData.qspline = normaliseIllumina(BSData, method="qspline", transform=NULL)* *BSData.norm = BSData.qspline* <>*design = model.matrix(~-1+factor(c(1,3,4,1,1,1,2,2,3,3,1,2,1,2,1,1,1,2,1,2,4,5, 1,5 ))) * *colnames(design) = c("NC", "C", "Ref", "FO", "NB")* <>*design[**1:24**,1:5]* <>* fit = lmFit(log2(exprs(BSData.norm)), design)* *con.matrix = makeContrasts (CvsNC = C-NC, CvsRef = C-Ref, NCvsRef = NC-Ref, levels = design) * <>*fit2 = contrasts.fit(fit, con.matrix)* *ebFit = eBayes(fit2)* *options(digits=3)* *tt = topTable(ebFit, n=dim(fit)[1], adjust="fdr")* *write.table(tt, "d:/Illumina_results/expset/3mCAN_b.txt ", row.names=F, quote=F, sep="\t")* -- Matthew Vitalone B.Sc (Hons) PhD(progress) NHMRC Centre for Clinical Research Excellence in Renal Medicine Centre For Transplant and Renal Research Transplant Laboratory (Room 2175), Clinical Sciences Westmead Millennium Institute Darcy Road, Westmead, NSW, 2145. Sydney, Australia. Phone: (+61-2) 9845 8906 Mobile: 0416 041783 Hosp. Page: 27147 Fax: (+61-2) 9633 9351 Email: matthew_vitalone at wmi.usyd.edu.au _______________________________________________________________ This electronic message and any attachments may be confidential. If you are not the intended recipient of this message would you please delete the message and any attachments and advise the sender. Western Sydney Area Health Services (WSAHS) uses virus scanning software but excludes any liability for viruses contained in any email or attachment. This email may contain privileged and confidential infor...{{dropped:10}}

On Sun, Feb 8, 2009 at 10:26 PM, Matthew Vitalone < matthew_vitalone@wmi.usyd.edu.au> wrote: > Hi all, > > I was just wanted someone to have a look at my rcode to see if it is > correct. This analysis (Illumina) arrays is comparing CAN vs No-CAN via a > common reference. I have 8 pts with CAN, 13 pts without CAN, 3 ref along > with 2 with fibrosis only (FO) and 2 with no histology (NB). The contrast i > wanted to make is CAN vs No-CAN via CAN vs Ref and NoCAN vs Ref. Ignoring > the FO and NB. I have put this code together from the Limma and beadarry > guides and general internet published documents. I get a diff exp list at > the end, I just wanted to check to see if the analysis is correct, so I can > trust the final diff exp list. Sorry, if the code i messy, I have taught > myself to use Limma and R. > Hi, Matt. I think the code looks pretty good. And, just for future reference, you do not (and should not, as it increases the variance used for testing for significance and, therefore, decreases power) go through a "common reference" when making comparisons, if you can help it. In your code, you have not done so, as your first contrast is cancer-non- cancer. The other two contrasts might be of interest for other reasons, but for the purposes of comparing cancer and non-cancer, your first contrast does that. Sean > > Thank you, > Matt > > Rcode: > > *library(limma)* > > *library(beadarray)* > > <>*setwd("d:/Illumina_results/expset")* <>* > > BSData = readBeadSummaryData("expset-QCgood-4.txt", qcFile=NULL, > sampleSheet=NULL,sep="\t",skip=8, ProbeID="PROBE_ID", columns = list(exprs = > "AVG_Signal", se.exprs="BEAD_STDERR", NoBeads = "Avg_NBEADS", > Detection="Detection Pval"), dec=".")* <>* > > BSData.qspline = normaliseIllumina(BSData, method="qspline", > transform=NULL)* > > *BSData.norm = BSData.qspline* > > <>*design = > model.matrix(~-1+factor(c(1,3,4,1,1,1,2,2,3,3,1,2,1,2,1,1,1,2,1,2,4, 5,1,5 > ))) * > > *colnames(design) = c("NC", "C", "Ref", "FO", "NB")* > > <>*design[**1:24**,1:5]* <>* > > fit = lmFit(log2(exprs(BSData.norm)), design)* > > *con.matrix = makeContrasts (CvsNC = C-NC, CvsRef = C-Ref, NCvsRef = > NC-Ref, levels = design) > > * > > <>*fit2 = contrasts.fit(fit, con.matrix)* > > *ebFit = eBayes(fit2)* > > *options(digits=3)* > > *tt = topTable(ebFit, n=dim(fit)[1], adjust="fdr")* > > *write.table(tt, "d:/Illumina_results/expset/3mCAN_b.txt ", row.names=F, > quote=F, sep="\t")* > > > -- > Matthew Vitalone B.Sc (Hons) PhD(progress) > NHMRC Centre for Clinical Research Excellence in Renal Medicine > Centre For Transplant and Renal Research > Transplant Laboratory (Room 2175), Clinical Sciences > Westmead Millennium Institute > Darcy Road, Westmead, NSW, 2145. > Sydney, Australia. > > Phone: (+61-2) 9845 8906 > Mobile: 0416 041783 > Hosp. Page: 27147 > Fax: (+61-2) 9633 9351 > Email: matthew_vitalone@wmi.usyd.edu.au_____________________________ __________________________________ > > This electronic message and any attachments may be confidential. If you > are not the intended recipient of this message would you please delete the > message and any attachments and advise the sender. Western Sydney > Area Health Services (WSAHS) uses virus scanning software but excludes any > liability for viruses contained in any email or attachment. > > This email may contain privileged and confidential infor...{{dropped:10}} > > _______________________________________________ > Bioconductor mailing list > Bioconductor@stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > [[alternative HTML version deleted]]