This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Hierarchical Clustering in single-channel agilent microarray experiment

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

r hierarchical clustering agilent
• 6,282 views • link
updated by Kevin Blighe Kevin Blighe • written by Leite Leite

Hi Leite,

I assume your question is about microarray data, but you never specified that. I adapted your title to clarify this.

but always from error

It's very hard for us to figure out what's going wrong if you don't show us the error message.

Cheers,
Wouter

• 0 views • link
written by WouterDeCoster WouterDeCoster

Hey WouterDeCoster,

Sorry for the lack of information, you're correct, its about microarray data analysis.

I've tryed this code for bloxplot:

boxplot(project.NormData, col=cols)
Error in sort.int(x, na.last = na.last, decreasing = decreasing, ...) : 
'x' must be atomic

And this code for Hierarchical Clustering:

plot(project.NormData)
Error in as.matrix(E$E) : object 'E' not found

Thank you so much,

Leite

• 0 views • link
written by Leite Leite

2 answers

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 ( C: Single-color Agilent array analyzing in R )

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)

box

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

violin

Hierarchical clustering (unsupervised on entire dataset - very CPU and memory intensive)

For a simple dendrogram or circular dendrogram, take a look at my threads:

dend

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)

heatmap

ComplexHeatmap

For ComplexHeatmap, see my recent post here: C: how to cluster genes in heatmap

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

• 0 views • link
written by Kevin Blighe Kevin Blighe

ComplexHeatmap

  ColAnn <- data.frame(condition)
  colnames(ColAnn) <- c("Condition")
  ColAnn <- HeatmapAnnotation(df=ColAnn,
    which="col",
    col=list(Condition=c("Wild-type"="yellow", "Knock-out"="royalblue")))

  #Boxplots for rows (genes) and columns (samples)
  boxAnnCol <- HeatmapAnnotation(
    boxplot=anno_boxplot(heat,
      border=FALSE,
      gp=gpar(fill="#CCCCCC"),
      lim=NULL,
      pch=".",
      size=unit(2, "mm"),
      axis=FALSE,
      axis_side=NULL,
      axis_gp=gpar(fontsize=12)),
    annotation_width=unit(c(1, 7.5), "cm"))
  boxAnnRow <- rowAnnotation(
    boxplot=row_anno_boxplot(heat,
      border=FALSE,
      gp=gpar(fill="#CCCCCC"),
      lim=NULL,
      pch=".",
      size=unit(3, "cm"),
      axis=FALSE,
      axis_side="top",
      axis_gp=gpar(fontsize=12)),
    annotation_width=unit(c(3), "cm"))

  myCol <- colorRampPalette(c("violet", "black", "springgreen"))(100)
  myBreaks <- seq(-3, 3, length.out=100)

  hmap <- Heatmap(heat,
            name="Expression Z-score",
            col=colorRamp2(myBreaks, myCol),
            heatmap_legend_param=list(
              color_bar="continuous",
              legend_direction="horizontal",
              legend_width=unit(5,"cm"),
              title_position="topcenter",
              title_gp=gpar(fontsize=15, fontface="bold")),

            #Row annotation configurations
            cluster_rows=TRUE,
            show_row_dend=TRUE,
            row_title="Gene",
            row_title_side="left",
            row_title_gp=gpar(fontsize=16, fontface="bold"),
            row_title_rot=90,
            show_row_names=TRUE,
            row_names_side="left",
            row_names_gp=gpar(fontsize=10),

            #Column annotation configuratiions
            cluster_columns=TRUE,
            show_column_dend=TRUE,
            column_title="Samples",
            column_title_side="top",
            column_title_gp=gpar(fontsize=16, fontface="bold"),
            column_title_rot=0,
            show_column_names=TRUE,
            column_names_gp=gpar(fontsize=14),

            #Dendrogram configurations: columns
            clustering_distance_columns=function(x) as.dist(1-cor(t(x))),
            clustering_method_columns="ward.D2",
            column_dend_height=unit(30,"mm"),

            #Dendrogram configurations: rows
            clustering_distance_rows=function(x) as.dist(1-cor(t(x))),
            clustering_method_rows="ward.D2",
            row_dend_width=unit(30,"mm"),

            #Annotations (row annotation must be added with 'draw' function, below)
            top_annotation_height=unit(0.5,"cm"),
            top_annotation=ColAnn,

            bottom_annotation_height=unit(3, "cm"),
            bottom_annotation=boxAnnCol)

    draw(hmap + boxAnnRow, heatmap_legend_side="left", annotation_legend_side="right")

complexheatmap

• 0 views • link
written by Kevin Blighe Kevin Blighe

Check out clustermap

• 0 views • link
written by Hussain Ather Hussain Ather

Log in to answer this question.

More posts like this