<p>I am using heatmap.2 from R to cluster microarry data. The figure below was made with </p>

<pre><code>heatmap.2(x, col=redgreen(75), scale="row", ColSideColors=col,  key=TRUE, symkey=FALSE,                             density.info="none",cexRow=1,cexCol=1,margins=c(6,11),  trace="none",srtCol=45)
</code></pre>

<p><img src='http://s24.postimg.org/5hbw04xh1/Sham2_T00d4h.png' alt='enter image description here' /></p>

<p>Looking at the expression of this figure I would think that the lowest three genes would be clustered closer to the upper 6 genes, as both groups down regulated in condition 1 (red) and upregulated in condition 2 (blue). Why are the first 6 genes clustered with the middle four, even though they show opposite patterns? The manpage of heatmap.2 says that it orders columns and rows on mean expression, so do the groups of genes that are clustered together have a similar average expression, and the pattern doesn't matter?</p>

<p>As you can see from the other replies, there are a couple ways of reading your question. The simpler version (answered partly by Mitch) could be rewritten, &quot;Regardless of the fact that the top 6 and bottom 3 genes cluster together, why aren&#39;t the dendrogram branchpoints rotated such that they&#39;re neighboring each other?&quot; The answer to that is found in <code>help(hclust)</code>, quoting:</p>

<p>In hierarchical cluster displays, a decision is needed at each merge to specify which subtree should go on the left and which on the right. ... The algorithm used in &lsquo;hclust&rsquo; is to order the subtree so that the tighter cluster is on the left (the last, i.e., most recent, merge of the left subtree is at a lower value than the last merge of the right subtree). Single observations are the tightest clusters possible, and merges involving two observations place them in order by their observation sequence number.</p>

<p>So that&#39;s why the tree is ordered why it is.</p>

<p>The other understanding of your question (addressed in a shorter version by Irsan) is, &quot;Since the top 6 and bottom 3 genes change in the same direction, why don&#39;t they cluster together?&quot; The answer to that is what Irsan wrote, namely that things are (by default) clustered with the &quot;complete&quot; method according to their euclidean distance. But of course, what does that actually mean? Probably the most important part of that is the <a href="http://en.wikipedia.org/wiki/Euclidean_distance">euclidean distance</a>. The easiest way to visualize that is graphically, so I&#39;ll make some test data that will more-or-less act like yours:</p>

<pre>
<code>library(gplots)
d &lt;- matrix(rnorm(78,c(rep(c(5,5,5,3,3,3),6), rep(c(2,2,2,3,3,3),4), rep(c(1,1,1,0,0,0),3)), 0.2), ncol=6, byrow=T)
</code></pre>

<p>The <code>d</code> matrix will then be something (the numbers are random) this:</p>

<pre>
<code>           [,1]     [,2]      [,3]       [,4]        [,5]        [,6]
 [1,] 5.1712874 4.946509 4.7496928 2.63230256  2.63394566  3.08378751
 [2,] 4.6592866 4.806637 4.9270387 3.00397576  2.96095830  3.05562337
 [3,] 5.1117932 4.982021 4.9432372 2.79170612  2.68319462  3.28125643
 [4,] 5.1939668 5.103573 4.9274765 2.75459710  2.92578504  3.17743957
 [5,] 4.8172162 4.869808 4.6627615 2.69654061  2.78337580  3.19926029
 [6,] 5.0833024 5.280635 5.0151975 2.79977773  3.18152508  3.04505210
 [7,] 2.2524162 2.079043 2.1589240 2.98382139  3.23926676  3.01663328
 [8,] 1.9869665 1.861576 2.1608973 2.94434318  2.50710978  2.75273211
 [9,] 1.8246908 2.119846 2.5934129 3.02386462  2.98645345  2.40102654
[10,] 2.2040333 2.248246 1.9449646 2.84486854  2.85842004  3.10360764
[11,] 1.1539225 1.005955 0.8953501 0.07245758 -0.03334591 -0.14624306
[12,] 0.8860263 0.752870 0.7780058 0.11035305 -0.22448713  0.07403528
[13,] 1.1470910 1.088321 1.0455616 0.15681927  0.06387492 -0.13126294
</code></pre>

<p>and a heatmap like this:</p>

<p><img alt="example heatmap" src="http://i.imgur.com/MVJQRqy.png" /></p>

<p>The Euclidean distance between two rows (<code>row_a</code> and <code>row_b</code>) is defined as <code>sqrt(sum(row_b-row_a)^2)</code>, using the vector functions in R (I only mention that because if you&#39;re used to thinking in, say, C, then this wouldn&#39;t make sense since these function only work with single values). You can get a nice visual of this by just plotting a subset:</p>

<pre>
<code>matplot(t(d), xlab=&quot;Sample&quot;, ylab=&quot;Expression Value&quot;)
</code></pre>

<p>This will look something like the following, with the rows points labeled given their row number (1-C instead of 1-13).</p>

<p><img alt="matplot output" src="http://i.imgur.com/O5DoIMi.png" /></p>

<p>As you can see, the top 6 rows in the matrix are really far away from the bottom 3 on this graph. In fact, the Euclidean distance between them is around <code>21</code>. The distance between the top 6 and middle 4 is only around <code>9</code> and the distance between the middle 4 and bottom 3 around <code>12</code>. These numbers match how closely things look by eye-ball. So, it&#39;s not so much that things are ordered by average expression but how closely they are to each other as measured by Euclidean distance.</p>

<p>The actual clustering procedure isn&#39;t so important in this example. The default is the <code>complete</code> method in <code>hclust</code>. Here, the distance between two clusters <code>A</code> and <code>B</code> is the maximum of the distances between the elements of <code>A</code> and <code>B</code>. One could instead use the <code>single</code> method, which takes the minimum distance between the elements of <code>A</code> and <code>B</code>, among many others. Each of these has its pros and cons, such as the <code>single</code> method being more prone to merging clusters if one has an outlier (in many cases, but see below). N.B., if you ever need to calculate distance matrices in R and use the <code>Ward</code> method, remember to <code>hclust(dist(x)^2, method=&quot;ward&quot;)</code> rather than <code>hclust(dist(x), method=&quot;ward&quot;)</code> to get the right numbers (<strong>most people don&#39;t do this and almost all tutorials that use the <code>ward</code> method are wrong!</strong>).</p>

<p><strong>Thought experiment:</strong> Suppose you did <code>d[8,6] &lt;- 1</code>. Given what I just wrote regarding the <code>complete</code> and <code>single</code> methods, how do you think the results could differ had you used the <code>single</code> method instead? This can be checked by comparing <code>plot(hclust(dist(d)))</code> and <code>plot(hclust(dist(d), method=&quot;single&quot;))</code>.</p>

<p><strong>Second though experiment:</strong> Suppose you instead did <code>d[1,1] &lt;- 6</code>. Which of the <code>single</code> and <code>complete</code> methods produce the results that make intuitive sense to you? Now you know why relying on <code>max</code> and <code>min</code> distances to compute distances between clusters can be problematic.</p>

<p>Not directly related to your question, but if you ever want to make similar heatmaps for RNAseq data, or any other data with a mean-variance relationship, remember that <strong>using Euclidian distance is a really really bad idea</strong>. For the reasons why, have a look at <a href="http://liorpachter.wordpress.com/2014/01/19/why-do-you-look-at-the-speck-in-your-sisters-quilt-plot-and-pay-no-attention-to-the-plank-in-your-own-heat-map/">this blog post</a> by Lior Pachter (auther of cufflinks, among other things).</p>
