Re-analysis of the data presented in Hinoue 2012

My previous attempts to identify CIMP in TCGA colorectal and colon datasets using 27k and 450k arrays failed. The purpose of this exercise was to understand precisely how the analyses are performed in already published papers. I chose the most recent paper on CIMP in CRC from Peter Laird's group at USC: "Genome-scale analysis of aberrant DNA methylation in colorectal cancer" by T. Hinoue et al published in Genome Research in 2011. Why I chose it:

1. They used Illumina 27k platform
2. The descriptions of their analyses are pretty clear, in addition Dr Hinoue has provided me with a Sweave file describing some of their analyses
3. They have the main clinical and technical information available: gender, age, tumor site, tumor stage and batch
4. Dr Hinoue has also provided a file with the mutation status of BRAF, KRAS and P53. Supplementary files that they didn't provide with the paper and weren't uploaded to GEO can be found on their group website here.
5. In the series matrix from GEO (GSE25062) they marked all probes excluded from the analysis as "null" (except those that belong to X and Y chromosomes, which was clear only from the Sweave file) which helped me in early identification of these probes on the platform instead of learning how to do it myself. It appeared later that it is very easy to do with R/Bioconductor package Genomic Ranges. Howto can be found here.

Few important points:

1. They verified the clusters that they obtained using RPMM algorithm (specific for beta-distributed data) by doing logit transformation of the data (into M value) and using R/Bioconductor package ConsensusClusterPlus (link to the paper describing the algorithm. The paper is well written and easy to understand. I found that I like the package based on this publication more than I like pvclust). Therefore since I personally prefer M values I focused on reproducing their results using this package rather than RPMM. 
2. It wasn't clear anywhere whether their results were generated with batch/age/gender adjusted data or only batch adjusted data. Only when I requested age information from Dr Hinoue he mentioned that 25 patients do not have any age information which made me realize that they didn't adjust for it. I wanted to understand the difference in obtained clusters with and without age/gender adjustment. 

Parameters: euclidean distance, hierarchical clustering, complete linkage, 10000 bootstraps (pvclust library)

Hinoue data: 125 tumor patients, batch removed, M value, remove 5060 "null probes" identified from the series matrix. Take 10% the most variable probes. Pvclust to identify cluster stability.

The best results I can get is 69 to 83 AU value for cluster stability even after 10,000 bootstraps. Does it mean that the clusters are weak and they just accepted it as is?

ConsensusCluster plus with the same parameters as pvlcust above

Will using a different package make any difference? Use ConsensusClusterPlus package with the same parameters as pvclust. Hierarchical clustering, evaluate 20 clusters, use 80% of the data for bootstrapping. They claimed that the identified 4 clusters. 

It seems that there is some separation but the size of the clusters is very uneven, nothing like was presented in the paper. 

After receiving the Sweave file I found that 5060 "null" probes don't include probes from XY chromosomes. Then I took the most variable probes I used for clustering and identified that 25% of those are from X or Y chromosomes. 

> icl[["clusterConsensus"]]
      k cluster clusterConsensus
 [1,] 2       1        0.8790997
 [2,] 2       2        0.8825873
 [3,] 3       1        0.8900365
 [4,] 3       2        0.8615545
 [5,] 3       3        0.8763479
 [6,] 4       1        0.7145161
 [7,] 4       2        0.8060535
 [8,] 4       3        0.9781181
 [9,] 4       4        0.9889430
[10,] 5       1        0.8289404
[11,] 5       2        0.7577152
[12,] 5       3        0.8221909
[13,] 5       4        0.7363796
[14,] 5       5        0.9454712
[15,] 6       1        0.8789223
[16,] 6       2        0.7593188
[17,] 6       3        0.7090342
[18,] 6       4        0.7150963
[19,] 6       5        0.9857516
[20,] 6       6        0.9189523

> icl[["clusterConsensus"]]
      k cluster clusterConsensus
 [1,] 2       1        0.9008927
 [2,] 2       2        0.9411635
 [3,] 3       1        0.8680856
 [4,] 3       2        0.8760303
 [5,] 3       3        0.7333630
 [6,] 4       1        0.7745014
 [7,] 4       2        0.7901681
 [8,] 4       3        0.8024310
 [9,] 4       4        0.8247401
[10,] 5       1        0.9153414
[11,] 5       2        0.8141924
[12,] 5       3        0.6573712
[13,] 5       4        0.6708655
[14,] 5       5        0.6319631
[15,] 6       1        0.8695192
[16,] 6       2        0.8568005
[17,] 6       3        0.7128842
[18,] 6       4        0.7407073
[19,] 6       5        0.5585342
[20,] 6       6        0.7426816