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",ctrlfiles="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}}

