[BioC] Simple affymetrix question (treated vs non-treated)

[Morten Mattingsdal]

Hi Wonjong,
I have some experience with affymetrix analysis in BioC, but im not a 
core member so my views are from a "end-user" point of view:
1. you design seems correct. If you want to double check it, use the 
affylmGUI library. Its a good way to check your design definitions.
2. Positive M values would mean up-regulated in group1

3. Yes. If you want to manually check probes to visually inspect their 
intensity values, use:
Eset=as.matrix(eset)
Eset["Pf.4.224.0_CDS_x_at",]

4. Ok.
5. Yes
6. If you have a relatively large experiment with many replicates, this 
is not necessary since weak spots tend to have a large variance and 
hence get poor p-values. Either way you can filter your data with 
modifications to the following code, where "data" is the return of the 
ReadAffy. This will keep probes with more than 3 present calls

Eset=as.matrix(eset)
calls=as.matrix(mas5calls(data))
calls.sum <- rowSums(calls=="P")
keep <- names(calls.sum[calls.sum>=3])
new=Eset[keep,]

then feed Limma the "new" object.

Your toptable is only 10 probes long. You probably have some 
up-regulated probes as well further down the list ?
To get all probes with positive B value use:
top <- topTable(fit2,coef=1, adjust="BH", sort.by="B", number=20000)
top <- top[top$B >0,]

NB have you done any MA plots or boxplots of your data ? maybe that can 
give you any hints

hope this helps
morten


Wonjong Moon wrote:
> I am trying to analyze Affymetrix data (treated (group1) vs Non-Treated
> (group2))
> I would like to know up-regulated probesets in group1 (treated)?
>
> 1. I would like to know that I set up the design matrix correctly.
> 2. I did group1 (treated)  - group1 (non-treated). So positive M values
> means up-regulated in group1 (treated), is that right? Or I switched
> treated to Non-treated?
> 3. Output numbers look like to have opposite signs.
> 4. If I go down to positive M values, then all P-values are 1.
> 5  Does 4 mean there's no up-regulated probe sets at significant level?
>
> Data and R codes are available http://binf.gmu.edu/wmoon/diff.
>
> 6. How can I remove 'Absent' flagged probesets? Or is it OK to leave
> them?
> Thank you.
>
> Wonjong
>
>
> Target file: 141PD.txt
>
> Name	FileName	Target
> CSA1	1A-1_SA1_141PD.CEL	CSA
> CSA2	2A-1_SA2_141PD.CEL	CSA
> CSA3	3A-1_SA4_141PD.CEL	CSA
> CSA4	4A-1_SA5_141PD.CEL	CSA
> Non5	5A-1_Non1_141PD.CEL	Non
> Non6	6Ar-1_Non2_141PD.CEL	Non
> Non7	7A-1_Non4_141PD.CEL	Non
> Non8	8A-1_Non5_141PD.CEL	Non
>
>
> library(limma) # Loads limma library.
> targets <- readTargets("141PD.txt") 
> library(affy); data <- ReadAffy(filenames=targets$FileName) # Reads CEL
> files 
> eset <- rma(data) # Normalizes data with 'rma' 
> esign <- model.matrix(~ -1+factor(c(1,1,1,1,2,2,2,2))) 
> colnames(design) <- c("group1", "group2") 
> fit <- lmFit(eset, design)
> contrast.matrix <- makeContrasts(group1-group2,levels=design) 
> fit2 <- contrasts.fit(fit, contrast.matrix)
> fit2 <- eBayes(fit2) # Computes moderated t-statistics and log-odds 
> topTable(fit2, coef=1, adjust="fdr", sort.by="M", number=10)
>
>
> output
>                          ID     M    A     t  P.Value adj.P.Val     B
> 21231   Pf.4.224.0_CDS_x_at -4.49 3.99 -37.7 2.59e-10  5.90e-06 10.27
> 21230   Pf.4.223.0_CDS_x_at -4.32 4.10 -30.8 1.32e-09  1.50e-05  9.77
> 22728           X03144.1_at -3.57 4.39 -23.4 1.17e-08  5.74e-05  8.86
> 22101    Pf.7.64.0_CDS_a_at -5.03 4.64 -23.2 1.22e-08  5.74e-05  8.84
> 612        AF306408.1_RC_at -3.51 3.52 -22.6 1.50e-08  5.74e-05  8.73
> 20063 Pf.13_1.84.0_CDS_a_at -5.03 4.54 -22.6 1.51e-08  5.74e-05  8.73
> 20855      Pf.2.36.0_CDS_at -2.58 4.65 -20.2 3.73e-08  1.21e-04  8.24
> 22524     Pf.9.267.0_CDS_at -4.16 4.05 -18.5 7.27e-08  2.04e-04  7.84
> 21350   Pf.5.119.0_CDS_x_at -4.55 5.32 -18.3 8.07e-08  2.04e-04  7.78
> 20078 Pf.13_1.99.0_CDS_x_at -2.36 6.50 -17.9 9.52e-08  2.17e-04  7.67
>
> 	[[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor
>
>
> .
>
>