Hi,

While finding differential gene expression of Affymetrix data, I am unable
to fix a small problem. The non-differentially regulated data are
showing bigger in the subsequent MA-plot.
Please correct me.

Thanks,
Santana

= = = = = = = == =

 myRMA <- justRMA()

 number_genes <- dim(exprs(myRMA))[1]

###  Filtering out uninteresting data

is_control <- logical(number_genes)
for (row in 1:number_genes)
{
if (charmatch("AFFX", rownames
  (exprs(myRMA))[row],
nomatch=0) == 0) is_control[row] <- TRUE
else is_control[row] <- FALSE
}

myRMA_no_controls <- myRMA[is_control]

# Filtering out the least variable genes, defined by the 90th percentile of
the distribution of CV-values

sd_values <- apply (exprs(myRMA_no_controls), MARGIN=1, FUN="sd")
mean_values <- apply (exprs(myRMA_no_controls), MARGIN=1, FUN="mean")
CV_values <- sd_values/mean_values
quantile_cut_off <- quantile(CV_values, probs=0.9)
myRMA_filtered <- myRMA_no_controls[CV_values>quantile_cut_off]


design <- cbind (Healthy = c(0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1),
                  Diseased=c(1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0))

fit <- lmFit (myRMA_filtered, design)

contrasts_matrix <- makeContrasts(Healthy-Diseased, levels = design)
fit2 <- contrasts.fit (fit, contrasts_matrix)
fit3 <- eBayes (fit2)

number_genes <- dim(exprs(myRMA_filtered))[1]
test_results <- topTable(fit3, number=number_genes, adjust="BH")

# False discovery rate cut off set to 0.05.

FDR_cutoff <- 0.05
p_values <- fit3$p.value
adjusted_p_values <- test_results$adj.P.Val

# Identify significant genes :

significant_genes <- test_results[test_results$adj.P.Val <= FDR_cutoff,]
gene_index <- rownames(significant_genes)

# MA-plot displaying the log fold change between diseased and healthy
samples as a function of the average expression level across all samples.

status <- character (length=number_genes)
status <- rep ("not changing", number_genes)
names (status) <- seq (1,number_genes,1)
status [gene_index] <- "significant"
plotMA (fit3, status=status, col=c("blue","pink"))  # the pink data-ponits
are coming relatively bigger

> sessionInfo()
R version 2.10.1 (2009-12-14)
i386-pc-mingw32

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] hgu133acdf_2.5.0     limma_3.2.1          affycoretools_1.18.0
annaffy_1.18.0       KEGG.db_2.3.5        GO.db_2.3.5          RSQLite_0.8-4
 [8] DBI_0.2-5            AnnotationDbi_1.8.0  affy_1.24.1
 Biobase_2.6.0

loaded via a namespace (and not attached):
 [1] affyio_1.14.0        annotate_1.24.0      biomaRt_2.2.0
 Biostrings_2.14.2    Category_2.12.0      gcrma_2.18.0
genefilter_1.28.0
 [8] GOstats_2.12.0       graph_1.24.1         GSEABase_1.8.0
IRanges_1.4.3        preprocessCore_1.8.0 RBGL_1.22.0          RCurl_1.3-1
[15] splines_2.10.1       survival_2.35-8      tools_2.10.1
XML_2.8-1            xtable_1.5-6

	[[alternative HTML version deleted]]

