A Guide to multiClust

Introduction

Whole transcriptomic profiles are useful for studying the expression levels of thousands of genes across samples. Clustering algorithms are used to identify patterns in these profiles to determine clinically relevant subgroups. Feature selection is a critical integral part of the process. Currently, there are many feature selection and clustering methods to identify the relevant genes and perform clustering of samples. However, choosing the appropriate methods is difficult as recent work demonstrates that no method is the clear winner. Hence, we present an R-package called multiClust that allows researchers to experiment with the choice of combination of methods for gene selection and clustering with ease. In addition, using multiClust, we present the merit of gene selection and clustering methods in the context of clinical relevance of clustering, specifically clinical outcome.

Our integrative R-package contains:

1. A function to read in gene expression data and format appropriately for analysis in R.

2. Four different ways to select the number of genes a. Fixed b. Percent c. Poly d. GMM

3. Four gene ranking options that order genes based on different statistical criteria a. CV_Rank b. CV_Guided c. SD_Rank d. Poly

4. Two ways to determine the cluster number a. Fixed b. Gap Statistic

5. Two clustering algorithms a. Hierarchical clustering b. K-means clustering

6. A function to calculate average gene expression in each sample cluster

7. A function to correlate sample clusters with clinical outcome

Function Workflow

The seven functions listed below should be used in the following order:

1. input_file, a function to read-in a text file containing the matrix of gene probes and samples to analyze.

2. number_probes, a function to determine the number of gene probes to select for in the feature selection process.

3. probe_ranking, a function to rank and select for probes using one of the available ranking options.

4. number_clusters, a function to determine the number of clusters to be used to cluster gene probes and samples.

5. cluster_analysis, a function to perform Kmeans or Hierarchical clustering analysis of the selected gene probe expression data.

6. avg_probe_exp, a function to produce a matrix containing the average expression of each gene probe within each sample cluster.

7. surv_analysis, a function to produce Kaplan-Meier Survival Plots after the discretization of samples into different clusters.

Other Functions Included

• WriteMatrixToFile, a function to write a data matrix to a text file.

• nor.min.max, a function to that uses feature scaling to normalize values in a matrix between 0 and 1.

1. Getting Started

In order to use multiClust, the user will need two text files.

• The first file is a gene probe expression dataset. This file should be a matrix with columns being the samples and the rows being genes or probes.

• The second file contains the patient clinical parameters. This file should be a matrix consisting of two columns. The first column contains the patient survival time and the second column contains the patient survival event occurrence. Survival time should be recorded in months. Survival event occurrence should be indicated by a 0 or 1. A 0 indicates censorship and a 1 indicates an event has occurred.

# Load GEO and biobase libraries
library(GEOquery)
library(Biobase)
library(multiClust)
    
# Obtain GSE series matrix file from GEO website using getGEO function
gse <- getGEO(GEO="GSE2034")
    
# Save the gene expression matrix as an object
data.gse <- exprs(gse)
    
# Save the patient clinical data to an object
pheno <- pData(phenoData(gse))
    
# Write the gene expression and clinical data to text files
WriteMatrixToFile(tmpMatrix=data.gse, tmpFileName="GSE2034.expression.txt",
    blnRowNames=TRUE, blnColNames=TRUE)
    
WriteMatrixToFile(tmpMatrix=pheno, tmpFileName="GSE2034.clinical.txt",
    blnRowNames=TRUE, blnColNames=TRUE)

# Load GEO and biobase libraries
library(GEOquery)
library(Biobase)
library(multiClust)
    
# Obtain GSE series matrix file from GEO website using getGEO function
gse <- getGEO(filename="GSE2034_series_matrix.txt.gz")

# Save the gene expression matrix as an object
data.gse <- exprs(gse)
    
# Save the patient clinical data to an object
pheno <- pData(phenoData(gse))
    
# Write the gene expression and clinical data to text files
WriteMatrixToFile(tmpMatrix=data.gse, tmpFileName="GSE2034.expression.txt",
    blnRowNames=TRUE, blnColNames=TRUE)
    
WriteMatrixToFile(tmpMatrix=pheno, tmpFileName="GSE2034.clinical.txt",
    blnRowNames=TRUE, blnColNames=TRUE)

# Load necessary libraries
library(GEOquery)
library(Biobase)
library(preprocessCore)
library(multiClust)
  
    
# Obtain GSE series matrix file from GEO website using getGEo function
gse <- getGEO(GEO="GSE2034")
    
# Save the gene expression matrix as an object
data.gse <- exprs(gse)

# Quantile normalization of the dataset
data.norm <- normalize.quantiles(data.gse, copy=FALSE)

# shift data before log scaling to prevent errors from log scaling
    # negative numbers
if (min(data.norm)> 0) {
    } 
    else {
        mindata.norm=abs(min(data.norm)) + .001
        data.norm=data.norm + mindata.norm  
    }

# Log2 scaling of the dataset
data.log <- t(apply(data.norm, 1, log2))

# Write the gene expression and clinical data to text files
WriteMatrixToFile(tmpMatrix=data.log,
    tmpFileName="GSE2034.normalized.expression.txt",
    blnRowNames=TRUE, blnColNames=TRUE)

# Obtain clinical outcome file
clin_file <- system.file("extdata", "GSE2034-RFS-clinical-outcome.txt",
    package="multiClust")

# Load in survival information file
clinical <- read.delim2(file=clin_file, header=TRUE)

# Display first few rows of the file
clinical[1:5, 1:2]

##   RFS_time RFS
## 1    77.88   0
## 2    49.32   1
## 3    130.2   0
## 4     82.8   0
## 5   144.96   0

# Column one represents the survival time in months
# Column two represents the relapse free survival (RFS) event 

2. Loading Your Gene Probe Expression Dataset into R

The first function of multiClust is used to load a gene expression dataset into R and format the matrix so that the probe or gene names are the rownames of the matrix.

# Load necessary libraries
library(multiClust)

## Loading required package: amap

# Obtain gene expression matrix
exp_file <- system.file("extdata", "GSE2034.normalized.expression.txt", package= "multiClust")

# Load the gene expression matrix 
data.exprs <- input_file(input=exp_file)

## [1] "The gene expression matrix has been loaded"

# View the first few rows and columns of the matrix
data.exprs[1:4,1:4]

##              GSM36777    GSM36778    GSM36779    GSM36780
## 1007_s_at 10.65626945 10.82617907 11.12609748 10.47785521
## 1053_at   6.293524052 6.105409137 6.239483636 6.628974931
## 117_at    6.785154455 6.446397235 6.779564325 6.673216272
## 121_at    8.147451246 8.757777784 8.530315085 8.499970427

3. Gene Selection Algorithms

# Example 1: Choosing a fixed gene or probe number
# Load necessary libraries
library(multiClust)

# Obtain gene expression matrix
exp_file <- system.file("extdata", "GSE2034.normalized.expression.txt",
    package="multiClust")

gene_num <- number_probes(input=exp_file,
    data.exp=data.exprs, Fixed=300,
    Percent=NULL, Poly=NULL, Adaptive=NULL)

## [1] "The fixed gene probe number is: 300"

# Example 2: Choosing 50% of the total selected gene probes in a dataset
# Load necessary libraries
library(multiClust)

# Obtain gene expression matrix
exp_file <- system.file("extdata", "GSE2034.normalized.expression.txt",
    package="multiClust")

gene_num <- number_probes(input=exp_file, 
    data.exp=data.exprs, Fixed=NULL,
    Percent=50, Poly=NULL, Adaptive=NULL)

## [1] "The percent gene probe number is: 150"

# Example 3: Choosing the polynomial selected gene probes in a dataset
# Load necessary libraries
library(multiClust)

# Obtain gene expression matrix
exp_file <- system.file("extdata", "GSE2034.normalized.expression.txt",
    package="multiClust")

gene_num <- number_probes(input=exp_file, 
    data.exp=data.exprs, Fixed=NULL,
    Percent=NULL, Poly=TRUE, Adaptive=NULL)

## [1] "The poly gene probe number is: 19"

# Example 4: Choosing the Adaptive Gaussian Mixture Modeling method
# Load necessary libraries
library(multiClust)

gene_num <- number_probes(input=exp_file, 
    data.exp=data.exprs, Fixed=NULL,
    Percent=NULL, Poly=NULL, Adaptive=TRUE)

# Load necessary libraries
library(multiClust)

# Obtain gene expression matrix
exp_file <- system.file("extdata", "GSE2034.normalized.expression.txt",
    package="multiClust")

# Load the gene expression matrix 
data.exprs <- input_file(input=exp_file)

## [1] "The gene expression matrix has been loaded"

# Call probe_ranking function
# Select for 500 probes
# Choose genes using the SD_Rank method
ranked.exprs <- probe_ranking(input=exp_file,
    probe_number=300, 
    probe_num_selection="Fixed_Probe_Num",
    data.exp=data.exprs, 
    method="SD_Rank")

## [1] "The selected SD_Rank Gene Expression text file has been written"

# Display the first few columns and rows of the matrix
ranked.exprs[1:4,1:4]

##             GSM36777 GSM36778 GSM36779 GSM36780
## 1405_i_at   2.663816 2.992651 1.833091 3.793433
## 200656_s_at 5.616488 2.372672 7.134110 5.639419
## 200641_s_at 4.165681 2.922955 4.397136 4.934420
## 200670_at   7.606198 6.223078 8.604945 7.279471

4. Cluster Analysis of Selected Genes and Samples

# Load necessary libraries
library(multiClust)

# Call the number_clusters function
# data.exp is the original expression matrix object outputted from
# the input_file function
# User specifies that samples will be separated into 3 clusters 
# with the "Fixed" argument
cluster_num <- number_clusters(data.exp=data.exprs, Fixed=3,
    gap_statistic=NULL)

## [1] "The fixed cluster number is: 3"

# Load necessary libraries
library(multiClust)

# Call the number_clusters function
# data.exp is the original expression matrix object ouputted from
# the input_file function
# User chooses the gap_statistic option by making gap_statistic equal TRUE
# The Fixed argument is also set to NULL
cluster_num <- number_clusters(data.exp=data.exprs, Fixed=NULL,
    gap_statistic=TRUE)

# Load necessary libraries

suppressPackageStartupMessages(library(ctc))
suppressPackageStartupMessages(library(gplots))
suppressPackageStartupMessages(library(dendextend))
library(ctc)
library(gplots)
library(dendextend)
library(graphics)
library(grDevices)
library(amap)
library(multiClust)

# Call the cluster_analysis function
hclust_analysis <- cluster_analysis(sel.exp=ranked.exprs,
    cluster_type="HClust",
    distance="euclidean", linkage_type="ward.D2", 
    gene_distance="correlation",
    num_clusters=3, data_name="GSE2034 Breast", 
    probe_rank="SD_Rank", probe_num_selection="Fixed_Probe_Num",
    cluster_num_selection="Fixed_Clust_Num")

## [1] "Your HClust Sample Dendrogram has been outputted"
## [1] "Your atr, gtr, and cdt files have been outputted for viewing in Java TreeView"
## [1] "A CSV file has been produced containing your sample and cluster assignment information"

# Display the first few columns and rows of the object
head(hclust_analysis)

## GSM36777 GSM36778 GSM36779 GSM36780 GSM36781 GSM36782 
##        1        3        1        1        1        2

# Load necessary libraries
library(multiClust)

# Call the cluster_analysis function
 kmeans_analysis <- cluster_analysis(sel.exp=ranked.exprs,
    cluster_type="Kmeans",
    distance=NULL, linkage_type=NULL, 
    gene_distance=NULL, num_clusters=3,
    data_name="GSE2034 Breast", probe_rank="SD_Rank",
    probe_num_selection="Fixed_Probe_Num",
    cluster_num_selection="Fixed_Clust_Num")

## [1] "A CSV file has been produced containing your sample and cluster assignment information"

 # Display the first few rows and columns of the object
 head(kmeans_analysis)

## GSM36777 GSM36778 GSM36779 GSM36780 GSM36781 GSM36782 
##        2        1        2        3        3        2

5. Obtaining the Average Expression for Each Gene/Probe in Each Cluster

The sixth function in this package avg_probe_exp is used to determine the average expression of each gene/probe for the samples in a particular cluster.

Function Arguments

1. The first argument “sel.exp” is an object containing the numeric selected gene/prone expression matrix. This object is an output of the probe_ranking function.

2. The second argument “samp_cluster” is an object vector containing the samples and the cluster number they belong to. This is an output of the cluster_analysis function.

3. The third argument “data_name” is a string indicating the cancer type and name of the dataset being analyzed.

4. The fourth argument “cluster_type” is a string indicating the type of clustering method used in the cluster_analysis function. “Kmeans” or “HClust” are the two options.

5. The fifth argument “distance” is a string describing the distance metric used for HClust in the cluster_analysis function. Options include one of “euclidean”, “maximum”, manhattan“,”canberra“,”binary“, or”minkowski“.

6. The sixth argument “linkage_type” is a string describing the linkage metric used in the HClust method in the cluster_analysis function. Options include “ward.D2”, “average”, “complete”, “median”, “centroid”, “single”, and “mcquitty”.

7. The seventh argument “probe_rank” is a string indicating the feature selection method used in the probe_ranking function. Options include “CV_Rank”, “CV_Guided”, “SD_Rank”, and “Poly”.

8. The eighth argument “probe_num_selection” is a string indicating the way in which probes were selected in the number_probes function. Examples include “Fixed_Probe_Num”, “Percent_Probe_Num”, “Poly_Probe_Num”, and “Adaptive_Probe_Num”.

9. The ninth argument “cluster_num_selection” is a string indicating how the number of clusters were determined in the number_clusters function. Examples include “Fixed_Clust_Num” and “Gap_Statistic”.

Example

Below is an example of how to use the avg_probe_exp function after performing Kmeans clustering analysis with the cluster_analysis function:

# Load necessary libraries
library(multiClust)

# Call the avg_probe_exp function
avg_matrix <- avg_probe_exp(sel.exp=ranked.exprs,
    samp_cluster=kmeans_analysis,
    data_name="GSE2034 Breast", cluster_type="Kmeans", distance=NULL,
    linkage_type=NULL, probe_rank="SD_Rank",
    probe_num_selection="Fixed_Probe_Num",
    cluster_num_selection="Fixed_Clust_Num")

## [1] "Your matrix containing the average gene probe expression in each cluster is finished"

# Display the first few rows and columns of the matrix
head(avg_matrix)

##             Cluster 1 Cluster 2 Cluster 3
## 1405_i_at    3.444507  3.295349  3.325252
## 200656_s_at  3.824532  4.635851  5.794335
## 200641_s_at  3.116710  4.061541  5.177569
## 200670_at    6.648387  7.571051  7.107753
## 200632_s_at  5.995347  5.470008  6.080703
## 200605_s_at  3.659888  5.028460  5.137496

Function Outputs

The outputs from this function are a text file containing the matrix of average expression of each probe/gene in each of the different clusters. In addition the object “avg_matrix” would also contain the average expression matrix.

6. Clinical Analysis of Selected Gene Probes and Samples

After dividing the selected gene probes and samples into their respective clusters, the surv_analysis function can be used to produce Kaplan-Meier Survival Plots. This function uses a Cox Proportional Hazard Model to portray the survival probabilities of the patients within each cluster over time.

Function Arguments

1. The first argument “samp_cluster” is an object vector containing the samples and the cluster number they belong to. This object is an output of the cluster_analysis function.

2. The second argument “clinical” is a string indicating the name of the text file containing patient clinical information. This file should be a matrix consisting of two columns. The first column contains the patient survival time information in months. The second column indicates the occurrence of a censorship (0) or an event (1).

3. The third argument “survival_type” is a string specifying the type of survival event being analyzed. Examples include “Disease-free survival (DFS)”, “Overall Survival (OS)”, “Relapse-free survival (RFS)”, etc.

4. The fourth argument “data_name” is a string indicating the name to be used to label the Kaplan-Meier Survival Plot.

5. The fifth argument “cluster_type” is a string indicating the type of clustering method used in the cluster_analysis function. “Kmeans” or “HClust” are the two options.

6. The sixth argument “distance” is a string describing the distance metric uses for HClust in the cluster_analysis function. Options include one of “euclidean”, “maximum”, manhattan“,”canberra“,”binary“, or”minkowski“.

7. The seventh argument “linkage_type” is a string describing the linkage metric used in the cluster_analysis function. Options include “ward.D2”, “average”, “complete”, “median”, “centroid”, “single”, and “mcquitty”.

8. The eighth argument “probe_rank” is a string indicating the feature selection method used in the probe_ranking function. Options include “CV_Rank”, “CV_Guided”, “SD_Rank”, and “Poly”.

The ninth argument “probe_num_selection” is a string indicating the way in which probes were selected in the number_probes function. Options include “Fixed_Probe_Num”, “Percent_Probe_Num”, “Poly_Probe_Num”, and “Adaptive_Probee_Num”.

The tenth argument “cluster_num_selection” is a string indicating how the number of clusters were determined in the number_clusters function. Options include “Fixed_Clust_Num” and “Gap_Statistic”.

Below is an example of how to use the surv_analysis function:

Example

# Load necessary libraries
library(gplots)
library(graphics)
library(grDevices)
library(multiClust)
library(survival)

# Obtain clinical outcome file
clin_file <- system.file("extdata", "GSE2034-RFS-clinical-outcome.txt",
    package="multiClust")

# Call the avg_probe_exp function
surv <- surv_analysis(samp_cluster=kmeans_analysis, clinical=clin_file,
    survival_type="RFS", data_name="GSE2034 Breast", 
    cluster_type="Kmeans", distance=NULL,
    linkage_type=NULL, probe_rank="SD_Rank",
    probe_num_selection="Fixed_Probe_Num", 
    cluster_num_selection="Fixed_Clust_Num")

## [1] "Your Kaplan Meier Survival Plot has been finished"

# Display the survival P value
surv

##    pvalue 
## 0.4061254

Function Outputs

This function outputs a pdf Kaplan-Meier Survival Plot of the patient survival time in months vs. the survival probability. In this example, the object “surv” will also contain the Cox Survival P value.

Note

The surv_analysis function also does not display the Kaplan-Meier plot in the console, but rather outputs the plot as a pdf file. Below is an example of what the plots would look like.

Kaplan-Meier Survival Plot

7. References

1. Davis, S. and Meltzer, P. S. GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics, 2007, 14, 1846-1847

2. Orchestrating high-throughput genomic analysis with Bioconductor. W. Huber, V.J. Carey, R. Gentleman, …, M. Morgan Nature Methods, 2015:12, 115.

3. Benjamin Milo Bolstad. preprocessCore: A collection of pre-processing functions. R package version 1.30.0.

4. Chris Fraley, Adrian E. Raftery, T. Brendan Murphy, and Luca Scrucca (2012) mclust Version 4 for R: Normal Mixture Modeling for Model-Based Clustering, Classification, and Density Estimation Technical Report No. 597, Department of Statistics, University of Washington

5. Antoine Lucas and Laurent Gautier (2005). ctc: Cluster and Tree Conversion.. R package version 1.42.0. http://antoinelucas.free.fr/ctc

6. Gregory R. Warnes, Ben Bolker, Lodewijk Bonebakker, Robert Gentleman, Wolfgang Huber Andy Liaw, Thomas Lumley, Martin Maechler, Arni Magnusson, Steffen Moeller, Marc Schwartz and Bill Venables (2015). gplots: Various R Programming Tools for Plotting Data. R package version 2.17.0. http://CRAN.R-project.org/package=gplots

7. Therneau T (2015). _A Package for Survival Analysis in S. version 2.38, <URL: http://CRAN.R-project.org/package=survival>.

8. Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M., Hornik, K.(2015). cluster: Cluster Analysis Basics and Extensions. R package version 2.0.3.

9. Tal Galili (2015). dendextend: Extending R’s Dendrogram Functionality. R package version 1.0.1. http://CRAN.R-project.org/package=dendextend

10. R Core Team (2015). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

11. Saldanha, A. J. “Java Treeview–extensible Visualization of Microarray Data.”Bioinformatics 20.17 (2004): 3246-248.

12. Antoine Lucas (2014). amap: Another Multidimensional Analysis Package. R package version 0.8-14. http://CRAN.R-project.org/package=amap