Dr. Shi, Thank you very much! I am doing it now as you suggested. The reason I asked for my protocols is that both the microarray data and the exon array data were from the same patients, but I got ~6000 significant probes (p.adjust()<0.05) from the exon array data analysis whereas only 12 significant probes/genes from the microarray data analysis. I am thinking that maybe I need to filter some probes for the affy exon array data, or maybe I need to do something to translate the ~6000 exon probes to genes as I am doing gene level analysis on exon array data. Could you give me some guidance on how to do it or whether I need to do that? Thanks, Xiayu From: Wei Shi [mailto:shi@wehi.EDU.AU] Sent: Sunday, July 07, 2013 5:35 PM To: Rao,Xiayu Subject: Re: illumina microarray and affymetrix exon array data Your analysis looks fine. But you should copy to the bioc mailing list when you post next time. On Jul 6, 2013, at 7:36 AM, Rao,Xiayu wrote: Dr. Shi, It is the first time our group analyze array data. I know that you are an expert in this field (Thank you soooo much for answering my previous question in bioconductor). Is it possible that you can have a look at my protocol and let me know if it is correct or not? I know my request may not be very specific, and I feel sorry for that, but it would be great if I can have input from you. Thank you very much! Main purpose: compare gene expression between tumor and normal samples (paired samples) (1) Illumina microarray - genomestudio "HumanHT-12 v4 Gene Expression BeadChip" - single channel : library(limma) x<-read.ilmn(files="irina_second_first_sample0222-analysis3.txt",ctrlf iles="control_probe_file_2.txt",other.columns="Detection") targets <- readTargets("Targets-3.txt") sample subject type 1 S652-N S652 N (Normal) 2 S652-T S652 T (Tumor) 3 S610-N S610 N 4 S610-T S610 T 5 S570-N S570 N 6 S570-T S570 T 7 S623-N S623 N 8 S623-T S623 T 9 S548-N S548 N 10 S548-T S548 T 11 S540-N S540 N 12 S540-T S540 T 13 S530-N S530 N 14 S530-T S530 T 15 S495-N S495 N 16 S495-T S495 T 17 S465-N S465 N 18 S465-T S465 T 19 S401-N S401 N 20 S401-T S401 T 21 S532-N S532 N 22 S532-T S532 T 23 S532-M S532 M (one metastasis sample) 24 S458-N S458 N 25 S458-T S458 T 26 S455-N S455 N 27 S455-T S455 T 28 S430-N S430 N 29 S430-T S430 T 30 S405-N S405 N 31 S405-T S405 T y <- neqc(x) expressed <- apply(y$other$Detection < 0.05,1,any) y <- y[expressed,] subject <- factor(targets$subject) type <- factor(targets$type, levels=c("N","T","M")) design <- model.matrix(~subject+type) fit <- lmFit(y, design) fit <- eBayes(fit) topTable(fit, coef="typeT",n=20) topTable(fit, coef="typeM",n=20) #heatmap 1: compare between Tumor and Normal ind1 <- p.adjust(fit$p.value[, 16], method = "BH") <0.05 z<-y[ind1,] z2<-z$E[1:12,] #remove NAs and retrieve $E expression numeric matrix heatmap.2(z2, col=greenred(75), scale="row", key=T, density.info="none", trace="none") #heatmap 2: compare between tumor subtypes based on the Tumor/Normal values of significant genes Calculation: 2^(value of Tumor - value of Normal) heatmap.2(z2_TN, col=greenred(75), scale="row", key=T, density.info="none", trace="none") (2) Affymetrix exon array (examine gene expression only) - Gene ST 1.x series of arrays - single channel: The samples are also paired with each patient having Tumor and Normal samples library(affy) mydata <- ReadAffy() pData(mydata) sample 6401_4C5-T Plus_(HuGene-1_0-st-v1).CEL 1 6402_4C2-N Plus_(HuGene-1_0-st-v1).CEL 2 6403_405D7-T Plus_(HuGene-1_0-st-v1).CEL 3 6404_405E2-N Plus_(HuGene-1_0-st-v1).CEL 4 6405_430A4-T Plus_(HuGene-1_0-st-v1).CEL 5 6406_430A3-N Plus_(HuGene-1_0-st-v1).CEL 6 6407_465A4-T Plus_(HuGene-1_0-st-v1).CEL 7 6408_464H8-N Plus_(HuGene-1_0-st-v1).CEL 8 6409_518F1-T Plus_(HuGene-1_0-st-v1).CEL 9 6410_518E8-N Plus_(HuGene-1_0-st-v1).CEL 10 6411_495E7-T Minus_(HuGene-1_0-st-v1).CEL 11 6412_495E5-N Minus_(HuGene-1_0-st-v1).CEL 12 6413_540E3-T Minus_(HuGene-1_0-st-v1).CEL 13 6414_540E2-N Minus_(HuGene-1_0-st-v1).CEL 14 6415_545F2-T Minus_(HuGene-1_0-st-v1).CEL 15 6416_545F1-N Minus_(HuGene-1_0-st-v1).CEL 16 6417_548A6-T Minus_(HuGene-1_0-st-v1).CEL 17 6418_548A4-N Minus_(HuGene-1_0-st-v1).CEL 18 6419_615B5-T Minus_(HuGene-1_0-st-v1).CEL 19 6420_615B8-N Minus_(HuGene-1_0-st-v1).CEL 20 6421_P12A3-T Minus_(HuGene-1_0-st-v1).CEL 21 6422_P12A2-N Minus_(HuGene-1_0-st-v1).CEL 22 image(data) hist(mydata[,1:22]) boxplot(mydata) eset <- rma(mydata) boxplot(exprs(eset)) #after normalization library(limma) subject <- c("S4","S4","S405","S405","S430","S430","S465","S465","S518 ","S518","S495","S495","S540","S540","S545","S545","S548","S548","S615 ","S615","SP12","SP12") type <- rep(c("T","N"),11) subject <- factor(subject) type <- factor(type) design <- model.matrix(~subject+type) fit <- lmFit(eset, design) fit <- eBayes(fit) topTable(fit, coef="typeT",n=20,adjust="BH") Heatmap1 > ind1 <- p.adjust(fit$p.value[, 12], method = "BH") <0.05 > z <- eset[ind1,] > z2 <- exprs(z) > heatmap.2(z2, col=greenred(75), scale="row", key=T, density.info="none", trace="none") Thanks a lot! Best regards, Xiayu ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:9}}