TP de la séance 4, Clustering

Introduction

Choisir son environnement de travail

Vous pouvez choisir de travailler soir sur le cluster core de l’IFB soit sur les ordinateurs de Paris-Diderot.

1. Sur le cluster de l’IFB

   • ouvrez une connexion au serveur RStudio https://rstudio.cluster.france-bioinformatique.fr/ et identifiez-vous

2. Sur les machines de la salle d’ordinateurs de Paris-Diderot

   • vous devez avoir la commande suivante dans votre .bashrc :

source /opt/sdv/anaconda/etc/profile.d/conda.sh`

• puis vous devez lancer l’environnement conda adéquat :

conda activate

• et enfin lancer le serveur Rstudio au moyen de la commande bash: rstudio

Dossier partagé contenant les données

Les données sont dans un répertoire partagé, dont le chemin dépend du serveur auquel vous êtes connectézs. Nous allons définir une variable data.folder qui indiquera le chemin de ce dossier partagé, en fonction du serveur.

Sur le serveur Rstudio de l’IFB-core-cluster, les données sont dans le répertoire /shared/projects/du_bii_2019/data/module3/seance4/BIC/.

## Liste des fichiers de données
data.folder <- "/shared/projects/du_bii_2019/data/module3/seance4/BIC/"

Contenu du dossier de données

Utilisez les commandes R suivantes:

• list.files() pour vérifier le contenu du dosser data.folder,
• file.size() pour calculer la taille de ces fichiers.

Astuces:

• list.files() retourne par défaut le nom de fichier, mais avec l’option full.names=TRUE vous obtiendrez le chemin complet.
• Calculez la taille des fichiers en bytes et en Megabytes (\(1Mb = 1024 \cot 1024 \cdot b\)), sachant que pour chaque conversion il faut diviser par 1024.
• Vous pouvez consulter notre solution à l’aide du code suivant (cliquer sur Code pour l’afficher).

# Return file sizes (in bytes)
message("Listing files in data folder: ", data.folder)

data.files <- list.files(path = data.folder) # List the data files
print(data.files)

 [1] "BIC_counts_all-genes.tsv.gz"                    "BIC_counts_filtered-genes.tsv.gz"               "BIC_edgeR_DEG_table.tsv"                        "BIC_edgeR_DEG_table.tsv.gz"                     "BIC_log2-norm-counts_all-genes.tsv.gz"          "BIC_log2-norm-counts_edgeR_DEG_fdr_0.01.tsv"   
 [7] "BIC_log2-norm-counts_edgeR_DEG_fdr_0.01.tsv.gz" "BIC_log2-norm-counts_edgeR_DEG_top_1000.tsv.gz" "BIC_log2-norm-counts_filtered-genes.tsv"        "BIC_log2-norm-counts_filtered-genes.tsv.gz"     "BIC_pheno.tsv.gz"                               "BIC_sample-classes.tsv.gz"                     

## Full path
data.path <- list.files(path = data.folder, full.names = TRUE) # List the data files
print(data.files)

 [1] "BIC_counts_all-genes.tsv.gz"                    "BIC_counts_filtered-genes.tsv.gz"               "BIC_edgeR_DEG_table.tsv"                        "BIC_edgeR_DEG_table.tsv.gz"                     "BIC_log2-norm-counts_all-genes.tsv.gz"          "BIC_log2-norm-counts_edgeR_DEG_fdr_0.01.tsv"   
 [7] "BIC_log2-norm-counts_edgeR_DEG_fdr_0.01.tsv.gz" "BIC_log2-norm-counts_edgeR_DEG_top_1000.tsv.gz" "BIC_log2-norm-counts_filtered-genes.tsv"        "BIC_log2-norm-counts_filtered-genes.tsv.gz"     "BIC_pheno.tsv.gz"                               "BIC_sample-classes.tsv.gz"                     

data.file.sizes <- file.size(data.path) # Get the size of the data files, in bytes

## Add file names 
## (for some strange reason, file.size returns a vector
## with no names, which is not very convenient)
names(data.file.sizes) <- data.files

## Compute file sizes in megabytes
data.file.Mb <- signif(digits = 2, data.file.sizes / (1024^2))
kable(data.frame(data.file.Mb))

Lire le tableau de valeurs d’expression

Nous allons maintenant lire le fichier d’expression. Pour cela, nous concaténons le chemin du dossier de données et le nom du fichier d’expressiion (BIC_log2-norm-counts_edgeR_DEG_top_1000.tsv.gz).

Ce fichier contient le comptages de lectures RNA-seq par gène, avec une sélection des gènes déclarés positifs pour letest de comparaison de moyennes multiples (voir-ci-dessus). Par ailleurs, nous avons arbitrairement appliqué un seuil supplémentaire en n’exportant que les 1000 gènes les plus significatifs, pour éviter de passer trop de temps sur le clustering hiérachique (complexité quadratique).

# Define the path of the expression file, by concatenating the data folder and the file name
# Note that file.path is convenient because it automatically used the appropriate parameter to separate the elements of a path (/ on Unix, \\ on Windows)
BIC.expr.file <- file.path(data.folder, "BIC_log2-norm-counts_edgeR_DEG_top_1000.tsv.gz")

# Load expression
message("Loading expression file\t", BIC.expr.file)
BIC.expr <- read.table(file = BIC.expr.file, header = TRUE)

Mesure de la taille des données

Prenez le temps d’identifier

• la taille du jeu de données
• le nombre d’individus
• le nombre de variables

Remarque : Classiquement, en analyse de données, les individus sont les lignes du tableau de données, les colonnes sont les variables.

Pour des raisons historiques, en analyse transcriptomique les données sont toujours fournies avec

• 1 ligne = 1 gène
• 1 échantillon biologique = 1 colonne

Cette convention a été établie en 1997, lors des toutes premières publications sur le transcriptome de la levure. Dans ces études, l’objet d’intérêt (l’“éindividu”) était le gène, et les variables étaient ses mesures d’expression dans les différentes conditions testées.

Pour l’anlayse de tissus cancéreux, on considère au contraire que l’“objet” d’intérêt est l’échantillon prélevé sur le patient, et les variables sont les mesures d’expression des différents gènes chez un patient.

Ce qui implique de faire attention, et éventuellement de travailler sur la matrice transposée (fonction t en R) pour utiliser correctement les fonctions classiques.

## Les données d'expression
dim(BIC.expr)

[1] 1000  819

## Les noms des lignes correspondent à des gènes
head(rownames(BIC.expr))

[1] "ENSG00000001617.11" "ENSG00000002834.17" "ENSG00000003989.16" "ENSG00000004838.13" "ENSG00000005302.17" "ENSG00000006016.10"

## Les noms des colonnes correspondent à des échantillons
head(colnames(BIC.expr))

[1] "X1AB92ADA.637E.4A42.A39A.70CEEEA41AE3" "DA98A67C.F11F.41D3.8223.1161EBFF8B58"  "X06CCFD0F.7FB8.471E.B823.C7876582D6FC" "A33B2F42.6EC6.4FB2.8BE5.542407A0382E"  "D021A258.8713.4383.9DCA.45E2F54A0411"  "C705FA90.D9AA.4949.BACA.1C022A14CB03" 

## transposons la table d'expression
message("Transcribing expression matrix")
BIC.expr.transposed <- t(BIC.expr)

Charger les étiquettes de classes des échantillons

Le fichier BIC_sample-classes.tsv.gzcontient les étiquettes de classes des échantillons.

BIC.sample.classes <- read.table(
  file.path(data.folder, "BIC_sample-classes.tsv.gz"),
  header = TRUE)

kable(BIC.sample.classes[1:10,])

Chaque échantillon a été assigné à une classe selon la combinaison de 3 marqueurs immunologiques:

• Estrogen Receptor 1 (ER1)
• Progesterone Receptor 1 (PR1)
• Human epidermal growth factor receptor 2 (Her2)

Utilisez

• La fonction R summary() pour compter le nombre de patientes positives / négatives pour chacun de ces trois marqueurs.

• La fonction R table() pour calculer le nombre d’échantillons de chaque type de cancer.

• La fonction R table() pour calculer une table de contiingence des marqueurs ER1 et PR1

summary(BIC.sample.classes)

       cancer.type        ER1            PR1            Her2    
 Basal.like  :131   Negative:184   Negative:267   Negative:631  
 HER2pos     : 41   Positive:635   Positive:552   Positive:188  
 Luminal.A   :422                                               
 Luminal.B   :118                                               
 Unclassified:107                                               

table(BIC.sample.classes$cancer.type)

  Basal.like      HER2pos    Luminal.A    Luminal.B Unclassified 
         131           41          422          118          107 

table(BIC.sample.classes$ER1, BIC.sample.classes$PR1)

          
           Negative Positive
  Negative      172       12
  Positive       95      540

table(BIC.sample.classes$ER1, 
      BIC.sample.classes$PR1, 
      BIC.sample.classes$Her2)

, ,  = Negative

          
           Negative Positive
  Negative      131        8
  Positive       70      422

, ,  = Positive

          
           Negative Positive
  Negative       41        4
  Positive       25      118

Projection ACP des échantillons

Nous allons réaliser une ACP sans mise à l’échelle.

## Transformation ACP
BIC.prcomp <- prcomp(BIC.expr.transposed, center = FALSE, scale. = FALSE)
names(BIC.prcomp) # check the name of the fields of the PCA object

[1] "sdev"     "rotation" "center"   "scale"    "x"       

par(mfrow = c(1,2))
## Plot de l'écart-type sur les premières composantes
barplot(BIC.prcomp$sdev[1:10], col = "#FFEEDD",
        main = "Standard deviation per PC")

## Visualisation des données
plot(BIC.prcomp, main = "Variance per PC", col="#BBDDFF")

par(mfrow=c(1,1))

Définissez une couleur pour chaque classe, et assignez à chaque échantillon la couleur correspondant à sa classe. Dessinez ensuite un nuage de points avec les coordonnées de chaque échantillon dans les 1ère et 2ème composantes (PC2 vs PC1)

# Assign a color to each cancer type
classes <- unique(BIC.sample.classes$cancer.type)
class.colors <- rainbow(n = length(classes))
names(class.colors) <- classes
data.frame(class.colors)

             class.colors
Luminal.A       #FF0000FF
Unclassified    #CCFF00FF
HER2pos         #00FF66FF
Luminal.B       #0066FFFF
Basal.like      #CC00FFFF

# Associate a color to each sample accordinig to its cancer type
sample.colors <- class.colors[BIC.sample.classes$cancer.type]
table(BIC.sample.classes$cancer.type, sample.colors)

              sample.colors
               #0066FFFF #00FF66FF #CC00FFFF #CCFF00FF #FF0000FF
  Basal.like           0         0         0         0       131
  HER2pos              0         0         0        41         0
  Luminal.A            0       422         0         0         0
  Luminal.B          118         0         0         0         0
  Unclassified         0         0       107         0         0

Question: comment interprétez-vous les barplots des écarts-types et variances pour les premières comosantes ? A discuter pendant le cours.

Clustering hiérarchique

## Compute Spearman correlation coefficient
BIC.cor <- cor(BIC.expr, method = "spearman")

# Check the dimensions of the correlation matrix
# (should be N x N, where N is the number of samples)
dim(BIC.cor)

[1] 819 819

## Derive a dissimilarity value from the corrrelation
BIC.dist <- as.dist(1 - BIC.cor)

## Run hierarchical clustering on the expression data
## (use default parameters)
BIC.hclust.complete <- hclust(BIC.dist)

## Plot the resulting tree
plot(BIC.hclust.complete, labels = F, hang = -1, las = 1, col = "grey",
     main = "Spearman dissimilarity, complete linkage")

##  Run hierarchical clustering with Ward agglomeration
BIC.hclust.ward <- hclust(BIC.dist, method = "ward.D2")

## Plot the resulting tree
plot(BIC.hclust.ward, labels = F, hang = -1, las = 1, col = "#88BB66", 
     main = "Spearman dissimilarity, Ward linkage")

ClassDiscovery::plotColoredClusters(
  BIC.hclust.ward, 
  labs=BIC.sample.classes$cancer.type, col = "grey",
  cols=sample.colors, cex=0.2)

legend("topright",
       legend=sort(classes), 
       col = class.colors,
       text.col = class.colors, cex=0.9)

Combinaison d’un arbre et d’une carte de température

La fonction R heatmap() permet de représenter à la fois les arbres produit par le clustering hiérarchique, et les profils d’expression.

Par défaut, elle effectue simultanément un clustering hiérarchique sur les lignes et sur les colonnes, ce qui permet de distinguer non seulement les groupes d’échantillons biologiques, mais également ceux de gènes.

Attention: la fonction heatmap() effectue par défaut un clustering hiérarchique sur les lignes et colonnes de votre matrice d’expression.

heatmap(as.matrix(BIC.expr), 
        distfun = function(x) as.dist(1 - cor(t(x), method = "spearman")), 
        hclustfun = function(x) hclust(x, method = "ward.D2"),
        labRow = NA, labCol = NA)

Heat map of the expression matrix clustered by genes (rows) and samples (columns).

# Compute mean value per gene
gene.means <- apply(BIC.expr, 1, mean)
gene.sd <-  apply(BIC.expr, 1, sd)

# Compute gene-wise centered expression values
BIC.expr.genes.centered <- BIC.expr - gene.means

# Compute gene-wise centred + scaled expression values
BIC.expr.genes.standardized <- BIC.expr.genes.centered / gene.sd

summary(unlist(BIC.expr.genes.standardized))

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-7.00569 -0.63658 -0.03651  0.00000  0.60201 12.86463 

hist(unlist(BIC.expr.genes.standardized), breaks = 100,
     main = "Gene-wise standardized expression values",
     xlab = "z-score", ylab = "Number of measures", 
     col = "#DDEEFF", xlim = c(-7,7))

# Draw a vertical grid
abline(v = -7:7, col = "#BBBBBB", lty = "dotted")

#Mark center
abline(v = 0, col = "#008800", lwd = 2)

# Mark dispersion (mean +- 1 sdev)
arrows(x0 = -1, y0 = 50000, x1 = 1, y1 = 50000, length = 0.05, angle = 30, code = 3, lwd = 2, col = "#008800")

Histogram of expression values after gene-wise standardization (centering and scaling).

## define a Blue - White - Red palette
frenchflag.palette <- colorRampPalette(c('dark blue','white','dark red'))

## Define a green - black - red palette
GBR.palette <- colorRampPalette(c('green','black','red'))

heatmap(as.matrix(BIC.expr.genes.standardized),
        zlim = c(-4,4),
        distfun = function(x) as.dist(1 - cor(t(x), method = "spearman")), 
        hclustfun = function(x) hclust(x, method = "ward.D2"),
        labRow = NA, labCol = NA,  # DO not print the labels (unreadable anyway)
        col = GBR.palette(100),
        scale = "none")

Heat map of the expression matrix clustered by genes (rows) and samples (columns).

# Define a color for the heatmap
expr.colors <- heat.colors(n = 100)

# Compute dendrogram with the hclust() function
# Choose a custom dissimilirity measure
# and agglomeration rule

# compute dissimilarity between samples
sample.dist <- as.dist(1 - cor(BIC.expr, method = "spearman")) 

# Run hierarchical clustering on samples
sample.clust = hclust(sample.dist, method = "complete")

# Convert the clustering result in a tree object that can be used with plot() and / or heatmap()
sample.tree <- as.dendrogram(sample.clust)

# Define colors for the cancer type
cancer.type <- unique(BIC.sample.classes$cancer.type)
type.cancer.colors <- rainbow(n = length(cancer.type))
names(type.cancer.colors) <- cancer.type

# Define colors for ER1 marker
ER1.classes <- unique(BIC.sample.classes$ER1)
ER1.colors <- rainbow(n = length(ER1.classes))
names(ER1.colors) <- ER1.classes
print(ER1.colors)

   Positive    Negative 
"#FF0000FF" "#00FFFFFF" 

# Define colors for PR1 marker
PR1.classes <- unique(BIC.sample.classes$PR1)
PR1.colors <- heat.colors(n = length(PR1.classes))
names(PR1.colors) <- PR1.classes
print(PR1.colors)

   Positive    Negative 
"#FF0000FF" "#FFFF00FF" 

# Define colors for Her2 marker
Her2.classes <- unique(BIC.sample.classes$Her2)
Her2.colors <- topo.colors(n = length(Her2.classes))
names(Her2.colors) <- Her2.classes

## Detiine annotations for the heatmap
## combining the 3 markers + cancer type
annot.tumeur.column = HeatmapAnnotation(
  df = BIC.sample.classes,
  col = list(cancer.type = type.cancer.colors,
             ER1 = ER1.colors, 
             PR1 = PR1.colors,
             Her2 = Her2.colors)
)

# A first heatmap
my.heatmap <- ComplexHeatmap::Heatmap(
  as.matrix(BIC.expr.genes.standardized), 
  name = "TCGA Breast Invasive Cancer", 
  # col = frenchflag.palette(100),
  col = brewer.pal(11,"RdBu"),
  column_title = "Samples",
  row_title = "Genes",
  cluster_columns = sample.tree,
  show_column_names = FALSE, 
  show_row_names = FALSE,
  bottom_annotation = annot.tumeur.column
)
draw(my.heatmap)

kmeans

1. faire un premier kmeans, par exemple, en prenant le nombre de groupe trouvé sur le hclust

2. faire une boucle pour trouver le nombre optimal de cluster, en calculant l’inertie intra totale en fonction du nombre de groupe kmeans()$totss [faire une boucle pour i allant de 1 à 10 for (i in 1:10) {}]

3. refaire le kmeans avec ce nombre optimal

4. visualiser ces groupes par exemple sur une projection des données dans le plan par PCA, à l’aide de la fonction prcomp().

Astuce: dans le résultat de prcomp(), les coordonnées des points se trouvent dans le champs x .

## Run k-means clustering with 20 centers
BIC.kmeans <- kmeans(BIC.expr.transposed, centers=20)

## Report the table of the clusters
table(BIC.kmeans$cluster)

  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
115  13  62  57   7  33  16  32   5  23  48  36  70  82  32   9  70  74  10  25 

T1 = Sys.time() # take time at the beginning of the task
I.intra = numeric(length=20)
I.intra[1] = kmeans(BIC.expr.transposed, centers=2)$totss
for (i in 2:20) {
  message("Running k-means with ", i, " centers")
  kmi <- kmeans(BIC.expr.transposed, centers=i)
  I.intra[i] <- kmi$tot.withinss
}

# Plot a curve showint the intra-cluster variance as a function of the number of clusters 
plot((1:20)-0.5, I.intra, type="h", lwd=2, col = "blue", 
     xlab = "k", ylab = "Intra-cluster variance",
     main = "k-mean: impact of k on the intra-cluster variance")

Intra-cluster variance plot. for a series of k-mean clustering with increasing values of k.

# Measure elapsed time
T2 = Sys.time()
Tdiff = difftime(T2,T1) ## Measure elapsed time

## Run k-means clustering with 2 clusters
BIC.kmeans2 <- kmeans(BIC.expr.transposed, centers = 2)
table(BIC.kmeans2$cluster) ## Cluster sizes

  1   2 
661 158 

## Run k-means clustering with 3 clusters
BIC.kmeans3 <- kmeans(BIC.expr.transposed, centers = 3)
table(BIC.kmeans3$cluster) ## Cluster sizes

  1   2   3 
244 435 140 

## Run k-means clustering with 10 clusters
BIC.kmeans10 <- kmeans(BIC.expr.transposed, centers = 10)
table(BIC.kmeans10$cluster) ## Cluster sizes

  1   2   3   4   5   6   7   8   9  10 
143  79  16 127  63  36  57 104  58 136 

## Cut the tree at some arbitrary levels to get clusters

## 2 clusters to see if the first subdivision corresponds to one of the 3 markers (ER1, PR1, Her2)
BIC.cutree2 <- cutree(BIC.hclust.ward, k = 2) 

## 5 clusters to see if the match the cancer types defined by biologists
BIC.cutree5 <- cutree(BIC.hclust.ward, k = 5)


## Define sample colors reflecting the cluster membership
hclust.k2.colors <- BIC.cutree2
kmeans.k2.colors <- BIC.kmeans2$cluster
kmeans.k10.colors <- BIC.kmeans10$cluster

## Define characters reflecting markers
pch.cancer.type <- as.numeric(BIC.sample.classes$cancer.type)
pch.er1 <- as.numeric(BIC.sample.classes$ER1)
pch.pr1 <- as.numeric(BIC.sample.classes$PR1)
pch.her2 <- as.numeric(BIC.sample.classes$Her2)

## Compare clusters and markers on the PC plot, to evaluate whether 
## the components capture relevant information
par(mfrow=c(1,2))
plot(BIC.prcomp$x[,1:2], 
     col = hclust.k2.colors, 
     pch = pch.er1,
     las = 1, cex = 0.7, 
     main = "Hclust 2 clusters versus markers")
legend("topright", cex = 0.8, legend = c("ER1+", "ER1-"), pch = )

plot(BIC.prcomp$x[,1:2], 
     col = kmeans.k10.colors, 
     pch = pch.er1,
     las = 1, cex = 0.7, 
     main = "K-means 10 clusters versus markers")

Intra-cluster variance plot. for a series of k-mean clustering with increasing values of k.

par(mfrow=c(1,2))
plot(BIC.hclust.ward, labels=FALSE, hang=-1, main = "Ward")

rect.hclust(BIC.hclust.ward, k=3)

BIC.cutree3 <- cutree(BIC.hclust.ward, k=3)

plot(BIC.hclust.ward, labels = F, hang = -1, main = "complete")

rect.hclust(BIC.hclust.ward, k = 10)

BIC.cutree10 <- cutree(BIC.hclust.ward, k = 10)

table(BIC.cutree3, BIC.kmeans$cluster)

           
BIC.cutree3   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20
          1 114   2  59   0   0   0  16  32   0   0  27  31  70  42   0   0   4  65  10   0
          2   0   0   0   0   7  33   0   0   5  23   0   0   0   0  31   9   0   0   0  24
          3   1  11   3  57   0   0   0   0   0   0  21   5   0  40   1   0  66   9   0   1

table(BIC.cutree10, BIC.kmeans$cluster)

            
BIC.cutree10   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20
          1  107   0  55   0   0   0   2   5   0   0  21   1   3  40   0   0   4  23   0   0
          2    0   0   0   0   7  33   0   0   0  13   0   0   0   0   3   9   0   0   0   1
          3    0   0   0  51   0   0   0   0   0   0   0   0   0   0   1   0   3   0   0   1
          4    1  11   3   6   0   0   0   0   0   0  21   5   0  40   0   0  63   9   0   0
          5    0   0   0   0   0   0  14   1   0   0   0   0   0   0   0   0   0   1   0   0
          6    0   0   0   0   0   0   0   0   0  10   0   0   0   0  28   0   0   0   0  23
          7    7   0   2   0   0   0   0  26   0   0   6  30   0   2   0   0   0  41   0   0
          8    0   2   2   0   0   0   0   0   0   0   0   0  66   0   0   0   0   0   0   0
          9    0   0   0   0   0   0   0   0   0   0   0   0   1   0   0   0   0   0  10   0
          10   0   0   0   0   0   0   0   0   5   0   0   0   0   0   0   0   0   0   0   0

par(mfrow = c(1,1))

clues::adjustedRand(BIC.cutree3, BIC.kmeans3$cluster)

     Rand        HA        MA        FM   Jaccard 
0.6742643 0.3286149 0.3297773 0.6058162 0.4341984 

clues::adjustedRand(BIC.cutree10, BIC.kmeans10$cluster)

     Rand        HA        MA        FM   Jaccard 
0.8481182 0.4189739 0.4228339 0.5141232 0.3377504 

clues::adjustedRand(BIC.cutree3, BIC.kmeans10$cluster)

     Rand        HA        MA        FM   Jaccard 
0.6586272 0.2323173 0.2344287 0.4545910 0.2348100 

clues::adjustedRand(BIC.cutree10, BIC.kmeans3$cluster)

     Rand        HA        MA        FM   Jaccard 
0.6745688 0.2536610 0.2557085 0.4747715 0.2824229 

## Correspondence between markers
kable(table(BIC.sample.classes$ER1, BIC.sample.classes$PR1))

fisher.test(table(BIC.sample.classes$ER1, BIC.sample.classes$PR1))

    Fisher's Exact Test for Count Data

data:  table(BIC.sample.classes$ER1, BIC.sample.classes$PR1)
p-value < 0.00000000000000022
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
  42.94805 165.19652
sample estimates:
odds ratio 
  80.90786 

# cut the tree at different levels (2, 5)
plot(BIC.hclust.ward, labels = FALSE, hang = -1)
rect.hclust(BIC.hclust.ward, k = 5, border = "blue")
rect.hclust(BIC.hclust.ward, k = 2, border = "red")
legend("topright", 
       legend = c("2 clusters", "5 clusters"), 
       col = c("red", "blue"), lwd = 2, cex = 0.7)

## hclust 2 clusters versus each marker
table(BIC.cutree2, BIC.sample.classes$Her2)

           
BIC.cutree2 Negative Positive
          1      508      179
          2      123        9

table(BIC.cutree2, BIC.sample.classes$ER1)

           
BIC.cutree2 Negative Positive
          1       67      620
          2      117       15

table(BIC.cutree2, BIC.sample.classes$PR1)

           
BIC.cutree2 Negative Positive
          1      142      545
          2      125        7

## K-means 2 clusters versus each marker
table(BIC.kmeans2$cluster, BIC.sample.classes$ER1)

   
    Negative Positive
  1       45      616
  2      139       19

table(BIC.kmeans2$cluster, BIC.sample.classes$PR1)

   
    Negative Positive
  1      119      542
  2      148       10

table(BIC.kmeans2$cluster, BIC.sample.classes$Her2)

   
    Negative Positive
  1      497      164
  2      134       24

## 5 clusters
table(BIC.cutree5, BIC.sample.classes$Her2)

           
BIC.cutree5 Negative Positive
          1      337       54
          2      123        9
          3       16       40
          4       96       63
          5       59       22

table(BIC.cutree5, BIC.sample.classes$ER1)

           
BIC.cutree5 Negative Positive
          1        4      387
          2      117       15
          3       42       14
          4        4      155
          5       17       64

table(BIC.cutree5, BIC.sample.classes$cancer.type)

           
BIC.cutree5 Basal.like HER2pos Luminal.A Luminal.B Unclassified
          1          3       0       298        51           39
          2        104       7         1         0           20
          3         11      28         2         5           10
          4          1       2        80        47           29
          5         12       4        41        15            9

cluster.vs.cancertype <- table(
  BIC.cutree5, 
  BIC.sample.classes$cancer.type)

kable(cluster.vs.cancertype, caption = "Hierarchical clusters (c = 5) versus cancer type. ")

## clustering versus cancer.type
## Compute the confusion table and adjusted RAND index (ARI)

## 2 clusters versus each marker
table(BIC.cutree2, BIC.sample.classes$ER1)

           
BIC.cutree2 Negative Positive
          1       67      620
          2      117       15

clues::adjustedRand(
  BIC.cutree2, 
  as.numeric(BIC.sample.classes$ER1))

     Rand        HA        MA        FM   Jaccard 
0.8195844 0.5810808 0.5813049 0.8707027 0.7688294 

table(BIC.cutree2, BIC.sample.classes$PR1)

           
BIC.cutree2 Negative Positive
          1      142      545
          2      125        7

clues::adjustedRand(
  BIC.cutree2, 
  as.numeric(BIC.sample.classes$PR1))

     Rand        HA        MA        FM   Jaccard 
0.7019742 0.3692338 0.3696219 0.7755577 0.6244903 

table(BIC.cutree2, BIC.sample.classes$Her2)

           
BIC.cutree2 Negative Positive
          1      508      179
          2      123        9

clues::adjustedRand(
  BIC.cutree2, 
  as.numeric(BIC.sample.classes$Her2))

       Rand          HA          MA          FM     Jaccard 
 0.53388801 -0.07618062 -0.07559919  0.66226224  0.49369937 

## 5 clusters versus cancer type
table(BIC.cutree5, BIC.sample.classes$cancer.type)

           
BIC.cutree5 Basal.like HER2pos Luminal.A Luminal.B Unclassified
          1          3       0       298        51           39
          2        104       7         1         0           20
          3         11      28         2         5           10
          4          1       2        80        47           29
          5         12       4        41        15            9

clues::adjustedRand(
  BIC.cutree5, 
  as.numeric(BIC.sample.classes$cancer.type))

     Rand        HA        MA        FM   Jaccard 
0.7105600 0.3330860 0.3348268 0.5451923 0.3743410 

## Negative control: compute the same stat with randomly permuted values
table(BIC.cutree5, BIC.sample.classes$cancer.type)

           
BIC.cutree5 Basal.like HER2pos Luminal.A Luminal.B Unclassified
          1          3       0       298        51           39
          2        104       7         1         0           20
          3         11      28         2         5           10
          4          1       2        80        47           29
          5         12       4        41        15            9

clues::adjustedRand(
  sample(BIC.cutree5), 
  as.numeric(BIC.sample.classes$cancer.type))

      Rand         HA         MA         FM    Jaccard 
0.57235104 0.01463136 0.01720336 0.32763791 0.19572626 

table(BIC.kmeans3$cluster, BIC.sample.classes$cancer.type)

   
    Basal.like HER2pos Luminal.A Luminal.B Unclassified
  1         22      31       111        41           39
  2          1       0       310        77           47
  3        108      10         1         0           21

clues::adjustedRand(
  BIC.kmeans3$cluster,
  as.numeric(BIC.sample.classes$cancer.type))

     Rand        HA        MA        FM   Jaccard 
0.6658905 0.2828718 0.2843799 0.5447041 0.3720078 

table(BIC.kmeans10$cluster, BIC.sample.classes[,1])

    
     Basal.like HER2pos Luminal.A Luminal.B Unclassified
  1           2       0       112        17           12
  2          12       4        39        15            9
  3           0       0        13         3            0
  4           0       0        99        21            7
  5           0       0        45         8           10
  6          23       1         1         2            9
  7           0       0        36        12            9
  8          85       6         0         0           13
  9           7      30         2         8           11
  10          2       0        75        32           27

clues::adjustedRand(BIC.kmeans10$cluster,
                    as.numeric(BIC.sample.classes[,1]))

     Rand        HA        MA        FM   Jaccard 
0.6788707 0.1385942 0.1418883 0.3302313 0.1726480 

# Visualisation


## hclust (2 groupes) et HER2+/-
plot(BIC.prcomp$x[,1:2], 
     col = hclust.k2.colors,
     pch = as.numeric(BIC.sample.classes$ER1), # character indicates ER1 status
     las = 1, cex = 0.7, 
     main = "Clusters versus markers")
grid()
texte.legend <- c("ER1-", "ER1+", "gr1", "gr2")
legend("topright", texte.legend, 
       col = c(1, 1, 1, 2), 
       pch = c(1, 2, NA, NA), 
       text.col = c(1, 1, 1, 2))

## Draw a confusion table of cluster versus ER1
kable(table(BIC.cutree2, BIC.sample.classes$ER1), caption = "hclust clusters versus ER1 marker")

Lister son environnement

A la fin de tout travail d’analyse, il est important de conserver une trace précise et complète de l’environnement précis utilisé pour produire les résultats.

## Print the complete list of libraries + versions used in this session
sessionInfo()

R version 3.5.2 (2018-12-20)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.3

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] ClassDiscovery_3.3.9  oompaBase_3.2.6       cluster_2.0.7-1       ComplexHeatmap_1.20.0 RColorBrewer_1.1-2    clues_0.5.9           FactoMineR_1.41       knitr_1.21           

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0           magrittr_1.5         MASS_7.3-51.1        leaps_3.0            mclust_5.4.2         scatterplot3d_0.3-41 colorspace_1.4-0     lattice_0.20-38      rjson_0.2.20         highr_0.7            stringr_1.4.0        tools_3.5.2          circlize_0.4.5       xfun_0.5            
[15] htmltools_0.3.6      yaml_2.2.0           digest_0.6.18        oompaData_3.1.1      GlobalOptions_0.1.0  shape_1.4.4          evaluate_0.13        rmarkdown_1.11       stringi_1.3.1        compiler_3.5.2       flashClust_1.01-2    GetoptLong_0.1.7