# Esercitazione: Analisi Dati Espressione genica 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 # OK
# 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
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. Clustering dei pazienti
# funzione per normalizzare i geni rispetto alla mediana e MAD
standardize <- function(z) {
  rowmed <- apply(z, 1, median)
  rowmad <- apply(z, 1, mad)
  rv <- sweep(z, 1, rowmed)
  rv <- sweep(rv, 1, rowmad, "/")
  return(rv)
}

# 5.1 Clustering gerarchico dei dati di espressione dell' oggetto ALLhm (81 probe sets e 15 campioni)
ALLhme <- exprs(ALLhm)
ALLdist1 <- dist(t(standardize(ALLhme)))
ALLhc1 <- hclust(ALLdist1, method="complete")



## plot del relativo dendrogramma:
plot(ALLhc1, xlab="", sub="", main="Clustering ALL1/AF4 e E2A/PBX1 patients (complete)")

# Analisi dell'oggetto ritornato da hclust:
ALLhc1$labels
ALLhc1$height
ALLhc1$merge

###################################################################
# modifica delle etichette dei pazienti per la visualizzazione del dendrogramma
ALLhc1$labels <- paste(ALLs$mol.biol, ALLhc1$labels)
ALLhc1$labels;
windows();
plot(ALLhc1, xlab="", sub="", main="Clustering ALL1/AF4 e E2A/PBX1 patients (complete)")
colnames(exprs(ALLs)) <- paste(ALLs$mol.biol, colnames(exprs(ALLs)))

#I risultati del clustering cambiano se si usa un algoritmo diverso? si provi con il method average ...

# 5.2 Clustering dei pazienti con geni selezionati in base alla loro variabilita' 
# selezione dei probe che hanno valor medio sopra la soglia del rumore nei due gruppi
ALLsub2 <- exprs(ALLs[(s1 | s2), ])
# calcolo del MAD
rowMads <- apply(ALLsub2, 1, mad)
# selezione dei probe con MAD>1.4
ALLsub2 <- ALLsub2[rowMads > 1.4, ]

ALLdist2 <- dist(t(standardize(ALLsub2)))
ALLhc2 <- hclust(ALLdist2)

windows(); plot(ALLhc2, xlab="", sub="", main="Clustering ALL1/AF4 e E2A/PBX1 (selezione MAD)")
# evidenziazione dei cluster nel plot e lista dei cluster ottenuti
list.hc2.2.clusters <- rect.hclust(ALLhc2, k = 2);
list.hc2.4.clusters <- rect.hclust(ALLhc2, k = 4, border=3);


###################################################################
# 6. Clusterizzazione con l'algoritmo k-means
# clusterizzazione con i geni selezionati come in 5.2
exprs.ALLsub2 <- t(standardize(ALLsub2))
kmeans.2 <- kmeans(exprs.ALLsub2, 2);
list.cluster.kmeans.2 <- list(k1 = kmeans.2$cluster[kmeans.2$cluster==1], k2 = kmeans.2$cluster[kmeans.2$cluster==2])

# i cluster individuati sono uguali a quelli trovati dal'algoritmo gerarchcico?

# provare piu' run delle due seguenti linee
kmeans.4 <- kmeans(exprs.ALLsub2, 4);
list.cluster.kmeans.4 <- list(k1 = kmeans.4$cluster[kmeans.4$cluster==1], k2 = kmeans.4$cluster[kmeans.4$cluster==2],
                              k3 = kmeans.4$cluster[kmeans.4$cluster==3], k4 = kmeans.4$cluster[kmeans.4$cluster==4])

# visualizzazione risultati kmeans tramite PCA
PCA<-prcomp(exprs.ALLsub2);
PC12 <- PCA$rotation[,1:2];  # selezionate le prime 2 PC
# proiezione sul sottospazio dlle prime 2 componenti principali
proj.exprs.ALLsub2 <- exprs.ALLsub2 %*% PC12;
windows(); plot(proj.exprs.ALLsub2, col = kmeans.2$cluster)  # visualizzazione 2 cluster
windows(); plot(proj.exprs.ALLsub2, col = kmeans.4$cluster)  # visualizzazione 4 cluster


# clustering sui dati proiettati sulle prime 2 componenti principali
# 2 cluster
kmeans.PCA.2 <- kmeans(proj.exprs.ALLsub2, 2);
list.cluster.kmeans.PCA.2 <- list(k1 = kmeans.PCA.2$cluster[kmeans.PCA.2$cluster==1], k2 = kmeans.PCA.2$cluster[kmeans.PCA.2$cluster==2])

windows();
plot(proj.exprs.ALLsub2, col = kmeans.PCA.2$cluster)
points(kmeans.PCA.2$centers, col = 1:2, pch = 8, cex=2);

# 4 cluster
kmeans.PCA.4 <- kmeans(proj.exprs.ALLsub2, 4);
list.cluster.kmeans.PCA.4 <- list(k1 = kmeans.PCA.4$cluster[kmeans.PCA.4$cluster==1], k2 = kmeans.PCA.4$cluster[kmeans.PCA.4$cluster==2],
                              k3 = kmeans.PCA.4$cluster[kmeans.PCA.4$cluster==3], k4 = kmeans.PCA.4$cluster[kmeans.PCA.4$cluster==4])

windows();
plot(proj.exprs.ALLsub2, col = kmeans.PCA.4$cluster)
points(kmeans.PCA.4$centers, col = 1:4, pch = 8, cex=2);


							  
							  
##############################################
# 7. Heatmap e clustering
# Visualizzazione congiunta delle heatmpa e del clustering gerarchico rispetto a geni e pazienti 
windows()
heatmap(exprs(ALLhm),  col=hmgreycol, ColSideColors=spcol)
windows()
heatmap(exprs(ALLhm),  col=hmcol, ColSideColors=spcol)


#Heatmap delle distanze: anche da qui si puo' desumere che esistono 2 cluster principali
windows(); heatmap(as.matrix(ALLdist2), sym=TRUE, col=hmcol, distfun=function(x) as.dist(x))
windows(); heatmap(as.matrix(ALLdist2), sym=TRUE, col=hmgreycol, distfun=function(x) as.dist(x))

######################################################

######################################################

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

library(hgu95av2.db);
gene.symbols <- unlist(mget(rownames(exprs(ALLhm)),hgu95av2SYMBOL));
gene.names <- unlist(mget(rownames(exprs(ALLhm)),hgu95av2GENENAME));
gene.ENTREZID <- unlist(mget(rownames(exprs(ALLhm)),hgu95av2ENTREZID));

# visualizare i risultati del clustering con i gene symbol
rownames(exprs(ALLhm)) <- gene.symbols;
windows()
heatmap(exprs(ALLhm),  col=hmgreycol, ColSideColors=spcol)

# salvare la heatmap in formato jpeg:
jpeg("heatmap.jpg")
heatmap(exprs(ALLhm),  col=hmcol, ColSideColors=spcol)
dev.off()

# salvare la heatmap in formato postscript:
postscript("heatmap.ps")
heatmap(exprs(ALLhm),  col=hmcol, ColSideColors=spcol)
dev.off()


###################################################
# 9. Correlazione dell'espressione con la posizione sui cromosomi
###################################################
library("geneplotter")
# Locazione cromosomica dei geni corrispondenti ai probe del chip  Affymetrix hgu95av2
chrLoc <- buildChromLocation("hgu95av2")
## cPlot(chrLoc, main="HG-U95Av2")


#Selezione dei geni con "elevata" espressione 
ALLch <- ALLs[s1|s2, ]
# calcolo dele medie per ciasuc gene separate per i due gruppi
m1  <- rowMeans(exprs(ALLch)[, ALLch$mol.biol=="ALL1/AF4"])
m2  <- rowMeans(exprs(ALLch)[, ALLch$mol.biol=="E2A/PBX1"])

# Calcolo dei decili dei dati combinati
deciles <- quantile(c(m1,m2), probs=seq(0,1,.1))
# trasformazione dei valori in intervalli:
s1dec <- cut(m1, deciles)
s2dec <- cut(m2, deciles)
# nomi dei geni (Affy ID)
gN <- names(s1dec) <- names(s2dec) <- featureNames(ALLch)



#Plot dell' espressione (quantizzata ad intervalli) lungo i cromosomi
# creazione di una palette di colori
colors <- brewer.pal(10, "RdBu")
# suddivsione del layout della finestra in 3 settori
layout(matrix(1:3,nr=1), widths=c(5,5,2))
# plot della locazione dei geni lungo i cromosomi per i campioni ALL1/AF4
cPlot(chrLoc, main="ALL1/AF4")
# plot dei colori lungo i cromosomi in accordo agli intervalli. Si noti che s1dec (un fattore) e' interpretato come indice intero
cColor(gN, colors[s1dec], chrLoc)
# plot della locazione dei geni lungo i cromosomi per i campioni E2A/PBX1
cPlot(chrLoc, main="E2A/PBX1")
# plot dei colori lungo i cromosomi in accordo agli intervalli.
cColor(gN, colors[s2dec], chrLoc)
# plot della legenda dei colori
image(1,1:10,matrix(1:10,nc=10),col=colors, axes=FALSE, xlab="", ylab="")
axis(2, at=(1:10), labels=levels(s1dec), las=1)


#Plot di un singolo cromosoma
par(mfrow=c(1,1))
# Smoothing degli strand positivi e negativi dei cromosomi
msobj <- Makesense(ALLs, "hgu95av2")
# plot del cromosoma 22
plotChr("22", msobj,  col = ifelse(ALLs$mol.biol=="ALL1/AF4", "red", "blue"), log=FALSE, xloc="equispaced")
# plot del cromosoma X
plotChr("X", msobj,  col = ifelse(ALLs$mol.biol=="ALL1/AF4", "red", "blue"), log=FALSE, xloc="equispaced")