Hey Kevin Blighe,
Thank you very much for your reply,
What worries me about the different results found is exactly the lack of reproducibility, because each one will use a specific code in the R with specific packages which are often not mentioned in the methods.
Between the example I quote the two articles use the same data and get different results (in numbers of genes), I reproduced the two tests and got different results (in number of genes). But one thing that might be equal between their findings are the functions and pathways of them can be generally the same.
As I am studying R, was this doubt, is "normal" this difference of results?
Original article method: log2fold change (FC) | > 1 and FDR <0.05 were considered the SDRs.
In this paper they found 208 genes, using my code I found 288.
My R code: #Install the core bioconductor packages, if not already installed
source("http://bioconductor.org/biocLite.R")
biocLite()
#Install additional bioconductor libraries, if not already installed
>biocLite("GEOquery")
>biocLite("affy")
>biocLite("gcrma")
>biocLite("affyPLM")
>biocLite("RColorBrewer")
>biocLite("limma")
>biocLite("hgu133plus2.db")
>biocLite("annotate")
>biocLite("simpleaffy")
#Load the necessary libraries
>library(GEOquery)
>library(affy)
>library(gcrma)
>library(affyPLM)
>library(RColorBrewer)
>library(hgu133plus2.db)
>library(annotate)
>library(simpleaffy)
#Downloading the file
>getGEOSuppFiles("**MyData**")
#Before starting to work, you have to unzip the files.
>untar("**MyData**/**MyData**_RAW.tar", exdir="data")
>cels <- list.files("data/", pattern = "[gz]")
>sapply(paste("data", cels, sep="/"), gunzip)
>cels
#Loading Files
#Phenodata must have names have to be the same as file names
>celfiles <- read.affy(covdesc="phenodata.txt", path="data")
#Normalising the data
>celfiles.rma <- rma(celfiles)
#Checking the quality
#Load colour libraries
>library(RColorBrewer)
#Set colour palette
>cols <- brewer.pal(8, "Set1")
#Plot a boxplot of unnormalised intensity values
>boxplot(celfiles, col=cols)
#Now a normalised boxplot
>boxplot(celfiles.rma, col=cols)
#Plot a density vs log intensity histogram for the unnormalised data
>hist(celfiles, col=cols)
#Plot a density vs log intensity histogram for the normalised data
>hist(celfiles.rma, col=cols)
#The first step of the analysis is to filter non-informative data
>celfiles.filtered <- nsFilter(celfiles.rma, require.entrez=FALSE, remove.dupEntrez=FALSE)
#Find differentially expressed probes
>samples <- celfiles.rma$Target
#Convert into factors
>samples <- as.factor(samples)
# Check factors have been assigned
>samples
#Set up the experimental design
>design <- model.matrix(~0 + samples)
>colnames(design) <- c("disease", "control")
#Inspect the experiment design
>design
#At this point we have normalized filtered data and a description of the data and the samples and the experimental design.
#Fit the linear model to the filtered expression set
>fit <- lmFit(exprs(celfiles.filtered$eset), design)
#Set up a contrast matrix to compare disease v control
>contrast.matrix <- makeContrasts(control-disease, levels=design)
#Now, the contrast matrix is combined with the linear fit model by set of probes.
>fit.con <- contrasts.fit(fit, contrast.matrix)
>fit.eb <- eBayes(fit.con)
>probeset.list <- topTable(fit.eb, coef=1, number=50000, adjust.method="BH",sort.by="P", lfc=1)
#Annotating the results with associated gene symbols
>gene.symbols <- getSYMBOL(rownames(probeset.list), "hgu133plus2")
>results <- cbind(probeset.list, gene.symbols)
>head(results)
#To save results
> write.table(results, "results.txt", sep="\t", quote=FALSE)
PS: I manually filter the FDR <0.05 in the results table.
Best regards
Leite