#####                                                      #####
# Practical training: Microarray data analysis                 #
#####                                                      #####

####
#  Part1 : Correction of experimental biases
####

# all graphs will be saved in a PDF file

pdf("normAnalysis.pdf")

# setwd("/home/cyril/session1_microarray/data_forStudent")
# reading of GPR file
GprFile = "dataFile1_normAnalysis.gpr"
library(marray)
rawdata = read.GenePix(GprFile)

# intensity values in red channel
rawdata@maRf
# intensity values in green channel
rawdata@maGf

# visualization of background signals
image(rawdata, xvar = "maRb", main = "Red background (with flags)")
image(rawdata, xvar = "maGb", main = "Green background (with flags)")

# spots that are flaggued will be eliminated from further analyses
# their intensity values are replaced by NA (missing values)
copyRawData <- rawdata
copyRawData@maRb[rawdata@maW < 0] = NA
copyRawData@maGb[rawdata@maW < 0] = NA

# visualization of background signals
image(copyRawData, xvar = "maRb", main = "Red background (without flag)")
image(copyRawData, xvar = "maGb", main = "Green background (without flag)")

# flag location on the slide
MyColor = maPalette(low = "blue", high = "white" , k = 10)
image(rawdata, xvar = "maW", col = MyColor, zlim = c(min(rawdata@maW), max(rawdata@maW)), main = "Location of Flags on the slide")

# spots that are flaggued are excluded from the analysis
rawdataWithoutFlags = rawdata

rawdataWithoutFlags@maRb[rawdataWithoutFlags@maW < 0] = NA
rawdataWithoutFlags@maGb[rawdataWithoutFlags@maW < 0] = NA
rawdataWithoutFlags@maRf[rawdataWithoutFlags@maW < 0] = NA
rawdataWithoutFlags@maGf[rawdataWithoutFlags@maW < 0] = NA

# no background correction is performed here
rawdataWithoutFlags@maGb[] = 0
rawdataWithoutFlags@maRb[] = 0

# comparison of Cy5 and Cy3 signals
plot(rawdataWithoutFlags, main = "MA plot before normalization")
boxplot(rawdataWithoutFlags, main = "Boxplot before normalization")

# three different procedures for signal normalization
rawdataWithoutFlagsNorm <- maNorm(rawdataWithoutFlags, norm = "median", echo = T)
rawdataWithoutFlagsNorm2 <- maNorm(rawdataWithoutFlags, norm = "loess", echo = T)
rawdataWithoutFlagsNorm3 <- maNorm(rawdataWithoutFlags, norm = "printTipLoess", echo = T)

# several plots allow for comparison of the normalization methods
plot(rawdataWithoutFlagsNorm, legend.func = NULL, main = "norm = Median")
plot(rawdataWithoutFlagsNorm2, legend.func = NULL, main = "norm = Loess")
plot(rawdataWithoutFlagsNorm3, legend.func = NULL, main = "norm = printTipLoess")

boxplot(rawdataWithoutFlagsNorm, main = "norm = Median")
boxplot(rawdataWithoutFlagsNorm2, main = "norm = Loess")
boxplot(rawdataWithoutFlagsNorm3, main = "norm = printTipLoess")

plot(density(maM(rawdataWithoutFlagsNorm2),na.rm = T), lwd = 2, col = 2, main = "Distribution of log(Ratio)")
lines(density(maM(rawdataWithoutFlags), na.rm = T), lwd = 2)
abline(v = 0)
legend(2,0.8,c("Before normalization","After normalization"), fill = c(1,2))

# close PDF image file
dev.off()


#####
#   Part2: Search for differentially expressed genes
#####

# Replicated data
dataFile <- "dataFile_diffAnalysis.txt"
data <- as.matrix(read.table(dataFile, row.names = 1, header = T))

# library LIMMA
library(limma)

# Linear model estimation
fit <- lmFit(data)

# Bayesian statistics
limma.res <- eBayes(fit)

# selection of differentially expressed genes (p-vlaue threshold = 0.01 here)
allgenes.limma <- topTable(limma.res, number = nrow(data))
siggenes.limma <- allgenes.limma[allgenes.limma[,4] < 0.01,]
dim(siggenes.limma[siggenes.limma[,2] > 0,])[1]
dim(siggenes.limma[siggenes.limma[,2] < 0,])[1]

write.table(siggenes.limma[siggenes.limma[,2] > 0,],
            row.names = T, quote = F, sep = "\t",
            file = "limma_up_signif_genes.txt")

write.table(siggenes.limma[siggenes.limma[,2] < 0,],
            row.names = T, quote = F, sep = "\t",
            file = "limma_low_signif_genes.txt")

# Volcano plot
volcanoplot(limma.res)