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]]

