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





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)

<a href="https://ibb.co/djsTjb"><img src="https://preview.ibb.co/i0Excw/box.png" alt="box" border="0"></a>

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

#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))

<a href="https://ibb.co/cPotjb"><img src="https://preview.ibb.co/kz3hcw/violin.png" alt="violin" border="0"></a>

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

#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)

<a href="https://ibb.co/dZHaxw"><img src="https://preview.ibb.co/frRDjb/dend.png" alt="dend" border="0"></a>

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

#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)

<a href="https://ibb.co/gGSdjb"><img src="https://preview.ibb.co/fb87cw/heatmap.png" alt="heatmap" border="0"></a>

#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.



