Entering edit mode
santana sarma
▴
80
@santana-sarma-3163
Last seen 10.2 years ago
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]]