Question: Hello everyone,

I'm trying to put a Hierarchical Clustering and boxplot in this code but always from error, what would be the best package to do this?

#Load limma
library(limma)

#Set-up
targetinfo <- readTargets("Targets.txt",row.names="FileName",sep="")

#Read files
project <- read.maimages(targetinfo,source="agilent", green.only=TRUE)

#Background correction
project.bgc <- backgroundCorrect(project, method="normexp", offset=16)

#Normalize the data with the 'quantile' method for 1-color
project.NormData <-normalizeBetweenArrays(project.bgc,method="quantile")

# load colour libraries
library(RColorBrewer)
# set colour palette
cols <- brewer.pal(8, "Set1")

#Histogram of non-normalized
plotDensities(project.bgc, col=cols, legend=FALSE)

#Histogram of normalized
plotDensities(project.NormData, col=cols, legend=FALSE)

#Create the study design and comparison model
design <- paste(targetinfo$Target, sep="")
design <- factor(design)
comparisonmodel <- model.matrix(~0+design)
colnames(comparisonmodel) <- levels(design)
#Checking the experimental design
design
comparisonmodel

project.fit <- lmFit(project.NormData, comparisonmodel)
project.fit <- lmFit(project.NormData,comparisonmodel)

#Applying the empirical Bayes method
project.fit.eBayes <- eBayes(project.fit)
names(project.fit.eBayes)

#Make individual contrasts and fit the new model
CaseControl <- makeContrasts(CaseControl="D0A-Control", levels=comparisonmodel)
CaseControl.fitmodel <- contrasts.fit(project.fit.eBayes, CaseControl)
CaseControl.fitmodel.eBayes <- eBayes(CaseControl.fitmodel)

#Filtering Results
nrow(topTable(CaseControl.fitmodel.eBayes, coef="CaseControl", number=99999, lfc=2))
probeset.list <- topTable(CaseControl.fitmodel.eBayes, "CaseControl", number=99999, adjust.method="BH", sort.by="P", lfc=2)

#To save results
write.table(probeset.list, "results.txt", sep="\t", quote=FALSE)

Best regards,

Leite

Answer: The expression values for a 2- or single-colour Agilent array are stored in the 'M' or 'E' variable, respectively, i.e., `project.NormData$M` or `project.NormData$E` ( https://www.biostars.org/p/287206/#287267 )

For each of the following functions, I encourage you to devote a full day to understanding what each and every parameter is doing. That is the best way for you to learn.

-------------------------------------------
-------------------------------------

#Box-and-whisker plot

par(mar=c(8,8,5,5), cex=1.0, cex.axis=1.4, cex.lab=1.4)
boxplot(project.NormData$E,
main="Box-and-whisker plot",
xlab="", ylab=bquote(~Log[2]~expression),
names=paste("Sample", c(1:ncol(project.NormData$E))),
col="skyblue",
las=2,
outline=FALSE)

-----------------------------------

#Violin plot (Wouter likes violin plots - maybe he plays a violin)

require(reshape2)
violinMatrix <- reshape2::melt(project.NormData$E)
colnames(violinMatrix) <- c("Gene","Sample","Expression")

library(ggplot2)
ggplot(violinMatrix, aes(x=Sample, y=Expression)) +
geom_violin() +
theme(axis.text.x = element_text(angle=45, hjust=1))

-----------------------------

#Hierarchical clustering (unsupervised on entire dataset - very CPU and memory intensive)
For a simple dendrogram or circular dendrogram, take a look at my threads:

- https://www.biostars.org/p/278032/#278045
- https://www.biostars.org/p/276306/#276358

require(pvclust)
pv <- pvclust(project.NormData$E,
method.dist="euclidean",
method.hclust="ward.D2",
nboot=100)
plot(pv)

--------------------------------

#heatmap.2 (hierarchical clustering dendrogram with heatmap)
For the heatmaps, you usually want to filter your expression matrix for genes that are differentially expressed. You appear to have just fitered out probes that are greater than absolute log (base 2) fold change 2, stored in your *probeset.list* object

#Filter the expression matrix to include only differentially expressed genes
sigmatrix <- project.NormData$E[probeset.list,]

#Scale the filtered expression matrix (convert to Z scale)
heat <- t(scale(t(sigmatrix)))

#Set colour
require(RColorBrewer)
myCol <- colorRampPalette(c("violet", "black", "springgreen"))(100)
myBreaks <- seq(-3, 3, length.out=101)

require("gplots")

#Euclidean distance; Ward's linkage
par(mar=c(1,1,1,1), cex=1.0)
heatmap.2(heat,
col=myCol,
breaks=myBreaks,
main="",
key=T, key.xlab="Expresssion\nZ-score", keysize=1.0,
scale="none",
ColSideColors=condition,
density.info="none",
reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean),
trace="none",
cexRow=1.0, cexCol=1.0,
distfun=function(x) dist(x, method="euclidean"),
hclustfun=function(x) hclust(x, method="ward.D2"),
margins=c(6, 6))
legend("top",
bty="n",
cex=1.0,
title="Condition",
c("Wild-type", "Knock-out"), fill=c("yellow", "royalblue"),
horiz=TRUE)

#1 - Pearson correlation distance; Ward's linkage
par(mar=c(1,1,1,1), cex=1.0)
heatmap.2(heat,
col=myCol,
breaks=myBreaks,
main="",
key=T, key.xlab="Expresssion\nZ-score", keysize=1.0,
scale="none",
ColSideColors=condition,
density.info="none",
reorderfun=function(d,w) reorder(d, w, agglo.FUN=mean),
trace="none",
cexRow=1.0, cexCol=1.0,
distfun=function(x) as.dist(1-cor(t(x))),
hclustfun=function(x) hclust(x, method="ward.D2"),
margins=c(6, 6))
legend("top",
bty="n",
cex=1.0,
title="Condition",
c("Wild-type", "Knock-out"), fill=c("yellow", "royalblue"),
horiz=TRUE)

#ComplexHeatmap

For ComplexHeatmap, see my recent post here: https://www.biostars.org/p/286187/#286507

I have also posted code in a comment following this answer.