---
title: "sigFeature: Significant feature selection using SVM-RFE & t-statistic"
author: "Pijush Das, Dr. Susanta Roychudhury, Dr. Sucheta Tripathy"
date: 'April 26, 2018'
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{sigFeature}
%\VignetteEngine{knitr::rmarkdown}
%\usepackage[utf8]{inputenc}
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## 1 Introduction
sigFeature is an R package which is able to find out the significant features
using support vector machine recursive feature elimination method (SVM-RFE)
(Guyon, I., et al. 2002) and t-statistic. Feature selection is an important
part dealing with machine learning technology. SVM-RFE is recognized as one
of the most effective filtering methods, which is based on a greedy algorithm
that only finds the best possible combination for classification without
considering the differentially significant features between the classes.
To overcome this limitation of SVM-RFE, the proposed approach is tuned to
find differentially significant features along with notable classification
accuracy. This package is able to enumerate the feature selection of any
two-dimensional (for binary classification) data such as a microarray etc.
This vignette explains the use of the package in a publicly available
microarray dataset.
## 2 Data
Example dataset in this package is provided from patients with squamous
cell carcinoma of the oral cavity (OSCC)(ACCNO: GSE2280 O'Donnell
RK et al. (2005)). Affymetrix Human Genome U133A Array is used for
genome-wide transcription analysis in this dataset. In O'Donnell
RK et al. (2005), the gene expression profiles obtained from primary
squamous cell carcinoma of the oral cavity (OSCC). Samples having
metastasis to lymph nodes (N+) is compared with non-metastatic samples
(N-). A total of 18 OSCCs were characterized for gene expression.
In their analysis, a predictive rule was built using a support vector
machine, and the accuracy of the rule was evaluated using cross-validation
on the original dataset. This was tested against four independent patient
samples. A signature gene set is produced which is able to predict correctly
the four independent patients. Using this method lymph node metastasis can
be predicted by gene expression profiles. The details of the study
can be found at:
<https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE2280>.
For "sigFeature" package evaluation, the microarray dataset (GSE2280) has
been classified into two classes such as lymph node metastatic (N+) and
no lymph node metastatic (N-) according to the TNM staging, provided in
the dataset. After downloading the dataset from Gene Expression Omnibus
(GEO) database firstly, it was normalized using "quantile" normalization
method by using the Bioconductor package "limma". There are altogether 27
samples (19 N+ and 8 N-). In order to reduce the runtime, t statistics was
computed between the N+ and N- samples for each feature (gene). The p
values less than 0.07 was used as a subset for testing. Now the expression
value of the sub-dataset is considered here as "x". The patients without
lymph node metastasis are represented -1, and the patients with lymph
node metastasis are represented as 1. Those -1 and 1 value is incorporated
with "y" as sample labels.
## 3 Example
In this section, we show the R script necessary to conduct the significant
feature selection with "sigFeature" package on the example dataset. First,
we load the "sigFeature" package and the example dataset (as mentioned in
the data section) into the memory of the machine. After that, step by step
execution takes place. Backward feature elimination procedure is a
time-consuming procedure. If we are dealing with the full dataset, we may
have to wait for a long time to process the dataset. Alternatively, we can
take advantage of parallel processing using the R package "parallel" to
reduce the processing time. To reduce the build-up time of this package
sometime, we use to process data in this vignette.
```{r cars}
library(sigFeature)
library(SummarizedExperiment)
data(ExampleRawData, package="sigFeature")
ExampleRawData
```
In this example, the variable "x" is containing the microarray expression
data, and the "y" contains the sample labels. The variable "x" is a
two-dimensional matrix in which the row contains the sample names, and
the column contains the gene/feature names. The variable "y" contains
the sample labels as -1 and 1.
```{r}
x <- t(assays(ExampleRawData)$counts)
y <- colData(ExampleRawData)$sampleLabels
```
Before going into details of the vignette, let us plot the unadjusted
p-value for each features which are estimated by using the expression
values of the two classes. The function named "sigFeaturePvalue()" is
used here to calculate the p-value using t-statistic. In this case,
we assume unequal variances and using the Welch approximation to
estimate the appropriate degrees of freedom.
```{r}
pvals <- sigFeaturePvalue(x,y)
hist(unlist(pvals),breaks=seq(0,0.08,0.0015),col="skyblue",
xlab="p value",ylab="Frequency",main="")
```
The idea of "sigFeature()" (Significant Feature Selection): The recursive
feature elimination that uses support vector machine can rank better
classifiers, but those classifiers may or may not be differently
significant between the classes. Features with notable classification
accuracy and also differentially significant between the classes have
a predominant role in a biological aspect. The algorithm, introduced
here is very efficient in ranking classifiers, as well as in determining
differently significant classes among them.
In data mining and optimization, feature selection is a very active field
of research where the SVM-RFE is distinguished as one of the most
effective methods. SVM-RFE is a greedy method that only hopes to find
the best possible combination for classification. In this proposed
algorithm, we have included a t-statistics along with SVM-RFE. The
primary objective of this algorithm is to enumerate the ranking weights
for all features and sort the features according to weight vectors as the
classification basis. The coefficient and the expression mean difference
between two compared groups of each individual feature is used to calculate
the weight value of that particular gene or the feature. The iteration
process is followed by backward removal of the feature. The iteration
process is continued until there is only one feature remaining in the
dataset. The smallest ranking weight will be removed by the algorithm
and stored in a stack, while remaining the feature variables have a
significant impact for producing weight. Finally, the feature variables
which are stored in the stack will be listed in the descending order
of descriptive different degree.
```{r}
#system.time(sigfeatureRankedList <- sigFeature(x, y))
```
The variable "sigfeatureRankedList" will return the address of the
genes/features after execution. To visualize the output, print command
can be used that will print the address of top 10 ranked genes/features.
```{r}
data(sigfeatureRankedList)
print(sigfeatureRankedList[1:10])
```
For producing the model vector for training the data set by using the
top 1000 genes/features, the R script is given as below. In this case
the package "e1071" is used. For producing the model vector, linear
kernel is chosen to get better performance.
```{r}
library(e1071)
sigFeature.model=svm(x[ ,sigfeatureRankedList[1:1000]], y,
type="C-classification", kernel="linear")
summary(sigFeature.model)
```
The R script below shows the prediction result of the same dataset
that was trained before. The user can also separate the dataset into
two parts e.g. training dataset and testing dataset and can implement
the algorithm.
```{r}
pred <- predict(sigFeature.model, x[ ,sigfeatureRankedList[1:1000]])
table(pred,y)
```
To solve the classification problem with the help of ranking the
features, another algorithm was proposed by Guyon et. al 2002 named
as SVM-RFE. In this algorithm, the dataset has been trained with SVM
linear kernel model to produce a weight vector for each feature. After
that, the feature with smallest rank is recursively removed and stored
in a stack data structure. The iteration process continues till the
last feature remains. This criterion is the w value of the decision
hyperplane given by the SVM. The weight vector is calculated by squiring
the variable w. Finally the features are selected in a list in the
descending order. The function "svmrfeFeatureRanking()" is included in
this package to compare the performance with "sigFeature()".
```{r}
#system.time(featureRankedList <- svmrfeFeatureRanking(x, y))
data(featureRankedList)
print("Top 10 features are printed below:")
print(featureRankedList[1:10])
```
For producing the model vector of the training dataset by using the
top 1000 genes/features, the R script is given below. The linear kernel
is chosen for producing the model vector.
```{r}
RFE.model=svm(x[ ,featureRankedList[1:1000]], y,
type="C-classification", kernel="linear")
summary(RFE.model)
```
The R script below will show the prediction result of the same dataset
which is trained before. The user can also separate the dataset into two
parts e.g. training dataset and testing dataset and can implement the
algorithm.
```{r}
pred <- predict(RFE.model, x[ ,featureRankedList[1:1000]])
table(pred,y)
```
In this example section we will observe the p value of the selected top
100 features. The function named "sigFeaturePvalue()" will calculate
the p value using t-statistic. The expression value will be taken from
the dataset of top 100 features.
```{r}
pvalsigFe <- sigFeaturePvalue(x, y, 100, sigfeatureRankedList)
pvalRFE <- sigFeaturePvalue(x, y, 100, featureRankedList)
par(mfrow=c(1,2))
hist(unlist(pvalsigFe),breaks=50, col="skyblue", main=paste("sigFeature"),
xlab="p value")
hist(unlist(pvalRFE),breaks=50, col="skyblue",
main=paste("SVM-RFE"), xlab="p value")
```
In this section a comparison of p values produced by using top 100
features (using "sigFeature" and "SVM-RFE" algorithm) is shown by
using scatter box plot.
```{r}
mytitle<-'Box Plot'
boxplot(unlist(pvalsigFe), unlist(pvalRFE), main=mytitle,
names=c("sigFeature", "SVM-RFE"),
ylab="p value", ylim=c(min(unlist(pvalsigFe)), max(unlist(pvalRFE))))
stripchart(unlist(pvalsigFe), vertical=TRUE, method="jitter", add=TRUE, pch=16,
col=c('green'))
stripchart(unlist(pvalRFE), vertical=TRUE, at=2, method="jitter", add=TRUE,
pch=16, col=c('blue'))
grid(nx=NULL, ny=NULL, col="black", lty="dotted")
```
In this section a heatmap is shown below with the expression value
of top 20 features generated by using "sigFeature()" function.
```{r}
library("pheatmap")
library("RColorBrewer")
pheatmap(x[ ,sigfeatureRankedList[1:20]], scale="row",
clustering_distance_rows="correlation")
```
In this section a heatmap is shown below with the expression value
of top 20 features generated by using "svmrfeFeatureRanking()"
function.
```{r}
pheatmap(x[ ,featureRankedList[1:20]], scale="row",
clustering_distance_rows="correlation")
```
For estimation of feature selection algorithms, Ambroise and McLachlan
(2002) suggested external cross-validation, in which the feature
selection and validation are performed on different parts of the
sample set, to obtain an unbiased performance estimate. The stratified
cross-validation (Breiman et al.,1984) is slightly different from regular
cross-validation. In k-fold stratified cross-validation, a data set is
randomly partitioned into k equally sized folds such that the class
distribution in each fold is approximately the same as that in the entire
data set. In contrast, regular cross-validation randomly partitions the
data set into k-folds without considering class distributions. Kohavi
(1995) reported that stratified cross-validation has smaller bias and
variance than regular cross-validation. The "sigFeature()" function
thus farther enhances by incorporating external cross-validation method.
In this external k-fold cross-validation procedure, k-1 folds are used
for selecting the feature, and one fold remains untouched which will
later use as a test sample set.
Detail description of the input parameters is given in the help file.
The set.seed function is used here only for the purpose of getting
the exact results when regenerating this vignette from its source
files. The function is given below.
```{r}
#set.seed(1234)
#results = sigFeature.enfold(x, y, "kfold", 10)
data("results")
str(results[1])
```
In this example section a new function is introduced named as
"sigFeatureFrequency()". The main purpose of this function is
to arrange the feature on the basis of its frequency. In the
previous section the dataset is divided into k-folds. Out of
which k-1 folds are randomly taken and used for selecting the
feature. This procedure is repeated k times and the frequency
values are calculated for top 400 genes from the k number of lists.
The "sigFeatureFrequency()" function is used here to rank further
all the k'th feature list according to its frequency. Details
description of this function is given in the help file.
```{r}
FeatureBasedonFrequency <- sigFeatureFrequency(x, results, 400, 400, pf=FALSE)
str(FeatureBasedonFrequency[1])
```
In this section external cross-validation is computed and finally
mean cross validation errors and the standard deviation of the
errors is returnd. The features, which are produced on the basis
of frequency is used for producing the mean cross validation error.
Training and test samples are same which were initially split from
the main sample set in to k-fold. In each iteration k-1 fold are
considered as training sample and remaining one fold is considered
as testing sample. In this external cross validation procedure
feature number is increased one by one and using that feature
training sample is trained. The number of miss classification is
averaged out for each fold change which finally represented as
error rate.
```{r}
#inputdata <- data.frame(y=as.factor(y) ,x=x)
#To run the code given bellow will take huge time. Thus the process
#data is given below.
#featsweepSigFe = lapply(1:400, sigCVError, FeatureBasedonFrequency, inputdata)
data("featsweepSigFe")
str(featsweepSigFe[1])
```
In this section the "PlotErrors()" function draw two plots.
The first plot is representing the average external cross
validation error and the second plot is representing the
standard deviation of the of the miss classifications.
```{r}
PlotErrors(featsweepSigFe, 0, 0.4)
```
Writing the results can be achieved with the "WritesigFeature()"
function. Please refer to the help file of "WritesigFeature()"
for more details on how to save to results to a specific file.
```{r}
#WritesigFeature(results, x)
```
## 4 SessionInfo
As last part of this document, we call the function "sessionInfo()",
which reports the version numbers of R and all the packages used
in this session. It is good practice to always keep such a record
as it will help to trace down what has happened in case that an R
script ceases to work because the functions have been changed in
a newer version of a package.
```{r}
sessionInfo()
```
## 5 References
1. Meyer, D., et al., e1071: Misc Functions of the Department of
Statistics (e1071), TU Wien. 2012.
2. O'Donnell RK, Kupferman M, Wei SJ, Singhal S et al. Gene expression
signature predicts lymphatic metastasis in squamous cell carcinoma of
the oral cavity. Oncogene 2005 Feb 10;24(7):1244-51.
3. Guyon, I., et al. (2002) Gene selection for cancer classification
using support vector machines, Machine learning, 46, 389-422.
4. Ambroise, C., & McLachlan, G. J. (2002). Selection bias in gene
extraction on the basis of microarray gene-expression data.
Proceedings of the national academy of sciences, 99(10), 6562-6566.
5. Breiman, L., Friedman, J., Olshen, R. and Stone, C. (1984) Classiflcation
and Regression Trees. Wadsworth and Brooks, Monterey, CA.
6. Kohavi, R. (1995, August). A study of cross-validation and bootstrap
for accuracy estimation and model selection. In Ijcai (Vol. 14,
No. 2, pp. 1137-1145).