--- title: "A Guide to multiClust" author: "Nathan Lawlor" date: '`r Sys.Date()`' output: html_document: toc: yes toc_depth: 6 pdf_document: toc: yes vignette: | %\VignetteIndexEntry{"A Guide to multiClust"} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ## 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. ### 1.1 Obtaining a Gene Expression Dataset and Clinical Information In this example, a gene expression dataset, GSE2034, will be obtained from the Gene Expression Omnibus (GEO) website . The following code can be used to extract the gene expression data and clinical information from the series matrix text zip file. When using this code, be sure that you have a strong internet connection. The GSE datasets are very large! ```{r, eval=TRUE, message=FALSE} # Load GEO and biobase libraries library(GEOquery) library(Biobase) library(multiClust) library(preprocessCore) library(ctc) library(gplots) library(dendextend) library(graphics) library(grDevices) library(amap) library(survival) ``` ```{r, eval=FALSE, message=FALSE} # 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[[1]]) # Save the patient clinical data to an object pheno <- pData(phenoData(gse[[1]])) # 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) ``` In the event that an error is produced when trying to obtain GSE datasets from the NCBI GEO FTP site, please use this alternative method: 1. On the NCBI GEO website (), simply search for the desired dataset (e.g."GSE2034") and scroll down to the "Analyze with GEO2R" section and click on "Series Matrix Files". 2. From there, download the series matrix txt.gz file to your computer. 3. Move the zipped file to a folder of your choice. 4. In R, set your working directory to the directory with the GSE file. 5. In R, read in the downloaded GSE series matrix file. After downloading the GSE series matrix file and moving it to a directory the following code can be used. ```{r, eval=FALSE} # 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[[1]]) # Save the patient clinical data to an object pheno <- pData(phenoData(gse[[1]])) # 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) ``` #### 1.2 Normalization of Gene Expression Datasets Some datasets on GEO are already normalized using Robust Multichip Average (RMA) procedure available in the 'affy' R pacakge or quantile normalization and log2 scaling using the 'preprocessCore' R package. For datasets that are not normalized, the following code can be used to extract the gene expression matrix from the series matrix text file, quantile normalize the data, log2 scale the data, and lastly write the data to a text file. ```{r, eval=FALSE} # 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) ``` #### 1.3 Formatting the Patient Clinical Information As mentioned in the beginning of this section, the patient survival information file should be a matrix consisting of two columns. However, when the patient clinical information is written to a text file, it is generally a matrix of more than two columns with all of the clinical parameters. It is recommended that this text file be loaded in R or Microsoft Excel to obtain the two columns of survival information needed. Below is an example of what the patient survival information file should look like: ```{r, eval=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] # 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. ### 2.1 Loading Text Files Containing Gene Expression Matrix In this example, the GSE2034 dataset is loaded into the R console using the `input_file` function. **Function Arguments** * This function has one variable called "input" which is the string name of the gene expression text file to load. **Example** ```{r, eval=TRUE} # 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) # View the first few rows and columns of the matrix data.exprs[1:4,1:4] ``` In this example, the gene expression matrix is stored into the object "data.exprs". --- ## 3. Gene Selection Algorithms ### 3.1 Determining the Number of Desired Probes or Genes The second function of this package, `number_probes` is a function used to specify the number of desired probes or genes the user would like to select for in their gene selection analysis. **Function Arguments** * This function has multiple arguments. The first argument "input" represents the string name of the expression matrix text file to be loaded into R. * The second argument "data.exp" is the object containing the expression matrix * The last four arguments "Fixed", "Percent", "Poly", and "Adaptive" are the four different methods that determine how the number of genes are selected. The default option in the function is set to "Fixed". 1. **Fixed** is an option where the user provides a positive integer to specify desired number of genes or probes to select for. The default value of "Fixed" is set to 1000. If the user writes: ```{r, eval=TRUE} # Example 1: Choosing a fixed gene or probe number # 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) ``` As a result, 300 genes from the dataset will be selected during the gene or probe selection process. 2. **Percent** is another option where the user provides a positive integer between 0 and 100 indicating the percentage of total genes or probes to select from the dataset. For example, if the dataset has 20,000 total probes and the user specifies: ```{r, eval=TRUE} # Example 2: Choosing 50% of the total selected gene probes in a dataset # 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) ``` As a result, 50% of the total dataset probes, or 300 probes, will be selected. 3. **Poly** is the third option which fits three second degree polynomial functions to the gene expression dataset based on the dataset's mean and standard deviation. When the user specifies: ```{r, eval=TRUE} # Example 3: Choosing the polynomial selected gene probes in a dataset # 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) ``` 4. **Adaptive** is the last option which implements Gaussian mixture modeling (GMM) to determine the appropriate number of genes or probes to select for. When the user specifies: ```{r, eval=FALSE} # Example 4: Choosing the Adaptive Gaussian Mixture Modeling method gene_num <- number_probes(input=exp_file, data.exp=data.exprs, Fixed=NULL, Percent=NULL, Poly=NULL, Adaptive=TRUE) ``` The "Mclust" function from the package mclust is used to apply Gaussian mixture modeling to the inputted gene expression matrix. In addition, files containing information about the dataset's mean, variance, mixing proportion, and gaussian assignment are also outputted. It should be noted that when choosing one of the three methods of gene or probe number selection, the other methods should be set equal to "NULL". Furthermore, the **Adaptive** option has a very long computational time. **Function Output** This function returns an object containing the number of genes/probes to select for during the feature selection process. #### 3.2 Choosing a Gene Selection Algorithm The third function of this package, `probe_ranking` is used to select the most informative gene probes within a dataset. This function allows the user choose from one of four different probe ranking algorithms. **Function Arguments** * This function has multiple arguments. The first argument "input" represents the string name of the expression matrix text file to be loaded into R. * The second argument "probe_number" is an output of the `number probes` function. For this argument, the user specifies the number of genes to select for. * The third argument "probe_num selection" is a string specifying which type of probe number selection method was used for the `number_probes` function. * The fourth argument "data.exp" is the object containing the expression matrix. * The last argument "method" is a string that specifies which of the four probe ranking methods to use. **List of the Five Probe Ranking Methods** 1. __CV_Rank__ is a gene probe ranking method that selects for probes using the coefficient of variation of the entire dataset (CV). 2. __CV_Guided__ is a gene probe ranking method that uses the coefficient of variation (CV) for each probe to select for probes. 3. __SD_Rank__ is a gene probe ranking method that selects for probes with the highest standard deviation within the dataset. 4. __Poly__ is a ranking method that fits three second degree polynomial functions of mean and standard deviation to the dataset to select the most variable probes in the dataset. This is the only method that does not use a probe number input from the `number_probes` function. When using this ranking method "probe_num" can be set to NULL. **Example** Below is an example of how to use the `probe_ranking` function: ```{r, eval=TRUE} # 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) # 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") # Display the first few columns and rows of the matrix ranked.exprs[1:4,1:4] ``` **Function Outputs** The output of this function is a text file containing a selected gene or probe expression matrix. In addition, the selected expression matrix is stored into the chosen variable "ranked.exprs". --- ## 4. Cluster Analysis of Selected Genes and Samples ### 4.1 Determining the Number of Clusters to Divide Samples Into The fourth function is this package, `number_clusters`, is a function used to designate the number of clusters that samples will be assigned to. **Function Arguments** 1. The "data.exp" argument represents an object containing the original gene expression matrix to be used for clustering of genes and samples. This object is an output of the `input_file` function. 2. The "Fixed" argument is a positive integer used to represent the number of clusters the samples and probes will be divided into. 3. The "gap_statistic" argument is a logical indicating whether to calculate the optimal number of clusters to divide the samples into using a gap statistic function `clus_Gap` from the package `cluster`. **Note** The user should only choose either the "Fixed" or "gap_statistic" option, not both. When using the gap_statistic option, change the argument to TRUE and "Fixed" to NULL. **Fixed Cluster Number Example** Below is an example of how to use the `number_clusters` function with a fixed number of clusters: ```{r, eval=TRUE} # 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) ``` The output from this example will be the object "cluster_num" with a value of 3. **Gap statistic Example** The `number_clusters` function can be called in a similar way to use the gap_statistic option: ```{r, eval=FALSE} # 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) ``` The gap_statistic option has a very long computational time and can take up to several hours. The output of this function will be a positive integer indicating the optimal number of clusters to divide samples into. #### 4.2 Kmeans or Hierarchical Clustering of Genes/Probes and Samples The fifth function in this package, `cluster_analysis` is a function used to perform Kmeans or Hierarchical clustering analysis of genes/probes and samples after undergoing a feature selection process in the `probe_ranking` function. **Function Arguments:** 1. The first argument "sel.exp", is an object containing the numeric selected gene expression matrix. This object is an output of the `probe_ranking` function. 2. The second argument "cluster_type" is a string indicating the type of clustering method to use. *Kmeans* or *HClust* are the two options. 3. The third argument "distance" is a string describing the distance metric to use for Hierarchical clustering via the `dist` function. `dist` uses a default distance metric of Euclidean distance. Options include one of "euclidean", "maximum", manhattan", "canberra", "binary", or "minkowski". Kmeans does not use a distance metric. 4. The fourth argument, "linkage_type" is a string describing the linkage metric to be used for hierarchical clustering in the `hclust` function. The default is "ward.D2", however other options include "average", "complete", "median", "centroid", "single", and "mcquitty". Kmeans does not use a linkage metric. 5. The fifth argument, "gene_distance" is a string describing the distance measure to be used for the `Dist` function when performing hierarchical clustering of genes. Options include one of "euclidean", "maximum", "manhattan", "canberra", "binary", "pearson", "abspearson", "correlation", "abscorrelation", "spearman" or "kendall". The default of gene_distance is set to "correlation". The argument can be set to NULL when Kmeans clustering is used. 6. The sixth argument "num_clusters" is a positive integer to specify the number of clusters samples will be divided into. This number is determined by the `number_clusters` function. 7. The seventh argument "data_name" is a string indicating the cancer type and/or name of the dataset being analyzed. This name will be used to label the sample dendrograms and heatmap files. 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". 9. 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_Probe_Num". 10. 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". **Function Outputs** * A CSV file containing the sample names and their respective cluster. An object containing a vector of the sample names and their cluster number is returned. * When hierarchical clustering is chosen as the cluster method, a pdf file of the sample dendrogram as well as atr, gtr, and cdt files for viewing in Java TreeView are outputted. In the Java TreeView heatmap, the samples are clustered by the method indicated by the "linkage_type" argument while genes are clustered by the method indicated by the "gene_distance" argument. The default of sample clustering for the Java TreeView heatmap is "ward.D2" and the default for gene clustering is centered pearson "correlation". **Hierarchical Clustering Example** Below is an example of how to use the `cluster_analysis` function using the Hierarchical clustering option: ```{r, eval=TRUE, fig.show='hold', warning=FALSE} # 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") # Display the first few columns and rows of the object head(hclust_analysis) ``` **Function Outputs** The outputs from this example would be an object "hclust_analysis" containing a vector of the sample names and cluster number. In addition, a sample dendrogram pdf file would be written. Lastly, the atr, gtr, and cdt files are outputted to view a heatmap of the genes and samples in Java TreeView. **Note** The `cluster_analysis` function does not display the sample dendrogram and heatmap in the console. To show what these images would look like, the following examples are provided below. **Java TreeView Heatmap Example** Below is an example of a heatmap produced in Java TreeView from the atr, gtr, abnd cdt files produced by this function. In this heatmap, rows represent genes and columns represent samples. ```{r, eval=TRUE, echo=FALSE, fig.show='hold', fig.width=8, fig.height=8, fig.align='center', warning=FALSE} # retrieve file path for image fpath <- system.file("extdata", "GSE2034_Breast.SD_Rank.ward.D2.Fixed_Probe_Num.Fixed_Clust_Num.jpg", package="multiClust") knitr::include_graphics(path = fpath) ``` **Sample Dendrogram Example** Below is an example of a sample dendrogram outputted by this function: ```{r, eval=TRUE, echo=FALSE, fig.show='hold', fig.width=8, fig.height=8, fig.align='center', warning=FALSE} # Produce a sample dendrogram using the selected gene/probe expression object # In the function, cluster_analysis the dendrogram is not displayed in # the console. # Make the selected expression object numeric colname <- colnames(ranked.exprs) ranked.exprs <- t(apply(ranked.exprs, 1,as.numeric)) colnames(ranked.exprs) <- colname # Normalization of the selected expression matrix before clustering norm.sel.exp <- t(apply(ranked.exprs, 1, nor.min.max)) hc <- hclust(dist(x=t(norm.sel.exp), method="euclidean"), method="ward.D2") dc <- as.dendrogram(hc) # Portray the samples in each cluster as a different color dc <- set(dc, "branches_k_color", value=1:3, k=3) # Make a vector filled with blank spaces to replace sample names with # blank space vec <- NULL len <- length(labels(dc)) for (i in 1:len) { i <- "" vec <- c(vec, i) } # Remove sample labels from dendrogram dc <- dendextend::set(dc, "labels", vec) # Plot the sample dendrogram plot(dc, main=paste("GSE2034 Breast", "Sample Dendrogram", sep=" ")) # Display orange lines on the plot to clearly separate each cluster rect.dendrogram(dc, k=3, border='orange') ``` **Kmeans Clustering Example** Below is an example of how to use the `cluster_analysis` function using the Kmeans clustering option: ```{r, eval=TRUE} # 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") # Display the first few rows and columns of the object head(kmeans_analysis) ``` **Function Outputs** The output from this example would solely be an object "kmeans_analysis" containing the vector of sample names and cluster number. There would be no sample dendrogram file. --- ## 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: ```{r, eval=TRUE} # 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") # Display the first few rows and columns of the matrix head(avg_matrix) ``` **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** ```{r, eval=TRUE, fig.show='hold'} # 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") # Display the survival P value surv ``` **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** ```{r, eval=TRUE, echo=FALSE, fig.show='hold', fig.width=8, fig.height=8, fig.align='center', warning=FALSE} # Obtain clinical outcome file clin_file <- system.file("extdata", "GSE2034-RFS-clinical-outcome.txt", package="multiClust") # Produce a Kaplan-Meier plot using the sample cluster information # from the cluster_analysis function # Read in survival data text file sample.anns <- read.delim2(file=clin_file, header=TRUE, stringsAsFactors=FALSE) time <- sample.anns[,1] event <- sample.anns[,2] time <- as.numeric(time) event <- as.numeric(event) event_type <- "RFS" time_type <- 'Months' # Determine max number of clusters in samp_cluster samp_cluster <- kmeans_analysis Number_Clusters <- max(samp_cluster) # Produce Kaplan Meier survival plots lev <- '1' for(i in 2:Number_Clusters) { lev <- c(lev, paste('',i, sep='')) } groupfac <- factor(samp_cluster, levels=lev) cox <- coxph(Surv(time, event==1) ~ groupfac) csumm <- summary(cox) cox.p.value <- c(csumm$logtest[3]) p <- survfit(Surv(time,event==1)~strata(samp_cluster)) plot (p,xlab= time_type ,ylab='Survival Probability',conf.int=FALSE, xlim=c(0,100), col=1:Number_Clusters, main=paste("GSE2034 Breast", "Kmeans", "SD_Rank", "RFS", sep =' ')) legend('bottomleft',legend=levels(strata(samp_cluster)),lty=1, col=1:Number_Clusters, text.col=1:Number_Clusters, title='clusters') legend('topright', paste('Pvalue =', signif(csumm$logtest[3],digits=2), sep ='')) ``` --- ## 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, . 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