# Esercitazione: Analisi Dati Espressione Differenziale in pazienti ALL


######################################
# 1. Caricamento dati ALL (Chiaretti et al., 2004)
library("ALL")
data("ALL")

# Informazioni sulle "etichette" dei pazienti
varLabels(ALL)
varMetadata(ALL)

########################################################################
# 2. Selezione di due sottogruppi:  ALL1/AF4 ed E2A/PBX1 
# the poor prognosis of adult patients is partly due to
# the presence of specific genetic translocations, such as the BCR/
# ABL or ALL1/AF4 gene rearrangements.
# The prevalence of the E2A-PBX1 fusion gene is one of the highest that has been described thus far in childhood acute
# lymphoblastic leukemia.
selSamples <- ALL$mol.biol %in% c("ALL1/AF4", "E2A/PBX1")
selSamples
# operazione di subsetting:
ALLs <- ALL[, selSamples]
ALLs
ALLs$mol.biol # ma le classi rimaste sono piu' di 2 ...
# riduzione a 2 soli fattori:
ALLs$mol.biol <- factor(ALLs$mol.biol)
ALLs$mol.biol #
# Abbiamo 15 campioni nell' oggetto ALLs

###########################################
# 3. Selezione dei geni differenzialmente espressi nei 2 gruppi.
# Vengono selezionati geni con espressione media maggiore di 100 nei 2 gruppi
# e con p-value (t-test) minore di 0.0002
library("genefilter")
meanThr <- log2(100)
g       <-   ALLs$mol.biol
# Selezione dei campioni ALL1/AF4 con media > meanThr  
s1 <- rowMeans(exprs(ALLs)[, g=="ALL1/AF4"]) > meanThr
# Selezione dei campioni E2A/PBX1 con media > meanThr  
s2 <- rowMeans(exprs(ALLs)[, g=="E2A/PBX1"]) > meanThr
# Selezione dei geni differenziamente espressi a livello significativita'=0.0002
res.t <- rowttests(ALLs, g);
s3 <- res.t$p.value < 0.0002
# Selezione dei geni che soddisfano la prima o la seconda condizione e la terza:
selProbes <- (s1 | s2) & s3;
AffyID.selected <- names(selProbes[selProbes==TRUE]);
ALLhm <- ALLs[selProbes,]
# selezionati 81 geni !

#################################################
# 4. Visualizzare i risultati: le heatmap
# costruzione di una palette costruita per interpolazione
library("RColorBrewer")
# preparazione di una palette "divergente" di colori
mypal<- brewer.pal(10, "RdBu")
display.brewer.pal(n=10, name="RdBu")

hmcol <- colorRampPalette(brewer.pal(10, "RdBu"))(256)
# si ottengono codici colori in HEX:
hmcol
# assegnamento dei colori "goldenrod", "skyblue" ai pazienti ALL1/AF4 e E2A/PBX1
spcol <- ifelse(ALLhm$mol.biol=="ALL1/AF4", "goldenrod", "skyblue")
heatmap(exprs(ALLhm), Rowv=NA, Colv=NA, col=hmcol, ColSideColors=spcol)



# Color palette toni di grigio ...
greypal<- brewer.pal(9, "Greys")
display.brewer.pal(n=9, name="Greys")
hmgreycol <- colorRampPalette(brewer.pal(9, "Greys"))(256)
windows()
heatmap(exprs(ALLhm), Rowv=NA, Colv=NA, col=hmgreycol, ColSideColors=spcol, reorderfun=NULL)



#######################################################################################
5. Come ottenere i gene.symbols e gene ENTREZ ID dagli identificatory Affymetrix dei probe set

library(hgu95av2.db);
AffyID.selected <- names(selProbes[selProbes==TRUE]);
gene.symbols <- unlist(mget(AffyID.selected,hgu95av2SYMBOL));
gene.names <- unlist(mget(AffyID.selected, hgu95av2GENENAME));
gene.ENTREZID <- unlist(mget(AffyID.selected, hgu95av2ENTREZID));


# Esercizio:
# Costruire una tabella con AffyID, gene symbol, gene name, EntrezID dei probe set selezionati
# Individuare i KEGG pathway associati ai geni selezionati (usare hgu95av2.db e KEGG.db).