A Gene Selection Method for Cancer Classification

Abstract

This paper proposes a method to select a set of genes from a large number of genes with the ability of classifying types of diseases. The proposed gene selection method is designed according to correlation analysis and the concept of 95% reference range. The method is very simple and uses the information of all genes. We have used the method in leukemia patients and achieved good classification results.

1. Introduction

In the clinical treatment of cancer, the corresponding treatment methods and measures are based on the correct classification of tumors. The traditional classification methods are based on tumor cell morphology, but sometimes tumors with the same histopathological type have different responses to treatment. So it becomes the current hot research topic to classify the tumor type using genomics [1, 2].

The gene chip experimental technology has provided a strong technology platform for tumor classification in genomics. By the discrimination analysis of gene expression levels in patients with different types of disease, the discrimination function established can be used to assist classification of clinical cases [3, 4].

As the gene chips have a very large number of genes, not all of these genes will provide information on the classification of type. When the expression levels of genes in different types of tissue samples do not change much, these genes have statistically no or only small discrimination capability. These genes are redundant.

Excluding these genes without classification capabilities will help to optimize the gene discrimination function which will be convenient for practical use. Therefore, it is necessary to select the genes with classification capabilities from a large number of genes [5].

Stepwise discrimination analysis is a commonly used statistical method for variable selection. However, for tissue samples containing thousands of genes, the stepwise discrimination analysis module in commonly used statistical software packages such as SAS and SPSS cannot function properly.

This paper presents a gene screening method that can analyze the classification capabilities of genes from thousands of genes and select the genes helpful to gene classification. The present new method has achieved good results in practice.

In the following sections we describe our studies on the statistical method for gene screening with discriminating ability in tumor classification.

In Methods we describe our new method for gene screening.

In Results we provide an application of our method in the study of the classification of leukemia patients.

Some concluding remarks are presented in Discussion.

2. Methods

Let a kind of disease have two subtypes 𝐴 and 𝐵. There are in total 𝑚 cases of the disease. Of these, 𝑛₁ cases are type 𝐴 and 𝑛₂ cases are type 𝐵. One tumor tissue sample was obtained from each of the patients. The expression levels of the 𝑛 genes 𝑔₁,𝑔,…,𝑔_𝑛 of each tissue sample were detected by gene chip.

2.1. The Correlation Coefficients Computation for Each Gene and the Classification Vector

If we list the 𝑛₁ cases of type 𝐴 first and the 𝑛₂ cases of type 𝐵 come in the tail, then the corresponding classification vector will be 𝑐=(𝑐₁,𝑐₂,…,𝑐_𝑚), where 𝑐_𝑗 corresponds to the case 𝑗. If the case 𝑗 is of type 𝐴, then 𝑐_𝑗=1, otherwise 𝑐_𝑗=0, 𝑗=1,2,…𝑚. Therefore, the classification vector has the form 𝑐=(1,1,^𝑛1…1,0,0,…^𝑛2,0).

Let the mean and standard deviation of the score for gene 𝑔_𝑖 in the 𝑛₁ tissue samples of type 𝐴 be 𝜇_𝐴(𝑔_𝑖) and 𝜎_𝐴(𝑔_𝑖), respectively. Similarly, the mean and standard deviation of the score for gene 𝑔_𝑖 in the 𝑛₂ tissue samples of type 𝐵 will be 𝜇_𝐵(𝑔_𝑖) and 𝜎_𝐵(𝑔_𝑖).

The correlation coefficient for gene 𝑔_𝑖 and the classification vector 𝑐 is defined as 𝑃𝑔_𝑖=𝜇,𝑐^𝐴𝑔_𝑖−𝜇_𝐵𝑔_𝑖𝜎_𝐴𝑔_𝑖+𝜎_𝐵𝑔_𝑖. (1)

The greater the absolute value of 𝑃(𝑔_𝑖,𝑐), the stronger the correlation of gene 𝑔_𝑖 and the classification vector 𝑐. In other words, the gene 𝑔_𝑖 has the ability to distinguish between type 𝐴 and type 𝐵.

From formula (1) we can compute 𝑃(𝑔₁,𝑐),𝑃(𝑔₂,𝑐),…,𝑃(𝑔_𝑛,𝑐) for genes 𝑔₁,𝑔,…,𝑔_𝑛.

2.2. Determine the Critical Value of Gene Screening

Let 𝑐^∗₁,𝑐^∗₂,…,𝑐^∗_𝑛 be 𝑛 random permutation vectors of the classification vector 𝑐=(1,…,1,0,…,0).

We now perform the following three steps of computation for each random permutation vector 𝑐^∗_𝑗, 𝑗=1,2,…,𝑛.(1)Collect the first 𝑛₁ cases in 𝑐^∗_𝑗 to a set denoted as class₁ and the remaining cases in 𝑐^∗_𝑗 to a set denoted as class₂.(2)For 𝑖=1,2,…,𝑛, compute 𝜇_class1(𝑔_𝑖,𝑐^∗_𝑗) and 𝜎_class1(𝑔_𝑖,𝑐^∗_𝑗) corresponding to 𝑔_𝑖 in class₁ and 𝜇_class2(𝑔_𝑖,𝑐^∗_𝑗) and 𝜎_class2(𝑔_𝑖,𝑐^∗_𝑗) corresponding to 𝑔_𝑖 in class₂, respectively.(3)For 𝑖=1,2,…,𝑛, compute 𝑃𝑔_𝑖,𝑐^∗_𝑗=𝜇_class1𝑔_𝑖,𝑐^∗_𝑗−𝜇_class2𝑔_𝑖,𝑐^∗_𝑗𝜎_class1𝑔_𝑖,𝑐^∗_𝑗+𝜎_class2𝑔_𝑖,𝑐^∗_𝑗. (2)

From the computation above, we obtain the correlation coefficients for the 𝑛 gene expression levels and the 𝑛 random permutation vector 𝑐^∗_𝑗, 𝑗=1,2,…,𝑛 as follows: 𝑃𝑔₁,𝑐^∗₁𝑃𝑔₁,𝑐^∗₂⋯𝑃𝑔₁,𝑐^∗_𝑛𝑃𝑔₂,𝑐^∗₁𝑃𝑔₂,𝑐^∗₂⋯𝑃𝑔₂,𝑐^∗_𝑛𝑃𝑔⋮⋮⋮⋮_𝑛,𝑐^∗₁𝑃𝑔_𝑛,𝑐^∗₂⋯𝑃𝑔_𝑛,𝑐^∗_𝑛. (3)

For a given value 𝑟 and vector 𝑐, denote the number of genes having correlation coefficients not less than 𝑟 as 𝑁₁(𝑐,𝑟) and the number of genes having correlation coefficients not greater than −𝑟 as 𝑁₂(𝑐,𝑟).

For all 1≤𝑖, 𝑗≤𝑛, we can define the following formulas: 𝑆₁𝑐^∗_𝑗,𝑔_𝑖=𝑔,𝑟1𝑃_𝑖,𝑐^∗_𝑗𝑔≥𝑟0𝑃_𝑖,𝑐^∗_𝑗𝑆<𝑟,₂𝑐^∗_𝑗,𝑔_𝑖=𝑔,𝑟1𝑃_𝑖,𝑐^∗_𝑗𝑔≤−𝑟0𝑃_𝑖,𝑐^∗_𝑗𝑁>−𝑟,₁𝑐^∗_𝑗,𝑟=^𝑛_𝑖=1𝑆₁𝑐^∗_𝑗,𝑔_𝑖,𝑁,𝑟₂𝑐^∗_𝑗,𝑟=^𝑛_𝑖=1𝑆₂𝑐^∗_𝑗,𝑔_𝑖,𝑟, (4)(5) where 𝑁₁(𝑐^∗_𝑗,𝑟) is the number of genes having correlation coefficients with vector 𝑐^∗_𝑗 not less than 𝑟 and 𝑁₂(𝑐^∗_𝑗,𝑟) is the number of genes having correlation coefficients with vector 𝑐^∗_𝑗 not greater than −𝑟, 𝑗=1,2,…,𝑛.

The right 5% quantile of the 𝑛 items 𝑁₁(𝑐^∗₁,𝑟),𝑁₁(𝑐^∗₂,𝑟),…𝑁₁(𝑐^∗_𝑛,𝑟) is denoted as 𝐿^𝑟₁ and the right 5% quantile of the 𝑛 items 𝑁₂(𝑐^∗₁,𝑟),𝑁₂(𝑐^∗₂,𝑟),…𝑁₂(𝑐^∗_𝑛,𝑟) is denoted as 𝐿^𝑟₂.

By increasing the value of 𝑟 gradually we can get 𝑟₁𝑟₂⋯𝑟_𝑘𝑁₁𝑐,𝑟₁𝑁₁𝑐,𝑟₂⋯𝑁₁𝑐,𝑟_𝑘𝐿^𝑟1₁𝐿^𝑟2₁⋯𝐿^𝑟𝑘₁. (6)

If we draw two curves of (𝑟_𝑖,𝑁₁(𝑐,𝑟_𝑖)) and (𝑟_𝑖,𝐿^𝑟𝑖₁) on the plane, we can see they have one intersection. The abscissa of the intersection is denoted as 𝑟₀ (see Figure 1).

Figure 1 

Two curves of (𝑟𝑖,𝑁1(𝑐,𝑟𝑖)) and (𝑟𝑖,𝐿𝑟𝑖1).

Similarly, for each 𝑟_𝑖 we also have 𝑟₁𝑟₂⋯𝑟_𝑘𝑁₂𝑐,𝑟₁𝑁₂𝑐,𝑟₂⋯𝑁₂𝑐,𝑟_𝑘𝐿^𝑟1₂𝐿^𝑟2₂⋯𝐿^𝑟𝑘₂. (7)

If we draw two curves of (𝑟_𝑖,𝑁₂(𝑐,𝑟_𝑖)) and (𝑟_𝑖,𝐿^𝑟𝑖₂) on the plane, we can see they have one intersection. The abscissa of the intersection is denoted as 𝑡₀ (see Figure 2).

Figure 2 

Two curves of (𝑟𝑖,𝑁2(𝑐,𝑟𝑖)) and (𝑟𝑖,𝐿𝑟𝑖2).

Let 𝑟^∗=max{𝑟₀,|𝑡₀|}. If |𝑃(𝑔_𝑖,𝑐)|≥𝑟^∗, then gene 𝑔_𝑖 is considered to have the ability to distinguish between type 𝐴 and type 𝐵. Therefore, it can be used as the index of the discrimination function for all 1≤𝑖, 𝑗≤𝑛.

3. Results

We have applied our method in the study of the classification of leukemia patients. We obtained a 38×7129 data matrix by taking tissue samples from 38 cases of clinically diagnosed leukemia patients and 7129 gene expression levels being detected for tissue samples of each case. Of the 38 cases, 27 cases had been diagnosed as acute lymphoblastic leukemia (ALL) and 11 cases are acute myeloid leukemia (AML). Before the two-type discrimination analysis, the 7129 genes are screened first using the method presented in this paper.

There are 38 components in the classification vector 𝑐=(1,…,1,0,…,0), whichrepresentsthe original classificationvector. The first 27 components of 𝑐 are 1 and the last 11 components of 𝑐 are 0. The correlation coefficients of the expression level for each gene and the classification vector 𝑐 are computed by formula (1). Table 1 is the frequency distribution table for the 7129 absolute values of the correlation coefficients.

Table 1 

The frequency distribution table for the 7129 correlation coefficients.

We generate 500 random permutation vectors 𝑐^∗_𝑗(𝑗=1,2,…,500) by 500 times of random permutation of the vector 𝑐=(1,…,1,0,…,0). The correlation coefficients 𝑃(𝑔_𝑖,𝑐^∗_𝑗) of the expression level for each gene and the classification vector 𝑐^∗_𝑗 are computed by formula (1) (1≤𝑖≤7129, 1≤𝑗≤500).

For the 6 values of 𝑟 between 0.1∼0.6 and all 𝑗 (1≤𝑗≤500), we compute 𝑁₁(𝑐^∗_𝑗,𝑟) and 𝑁₂(𝑐^∗_𝑗,𝑟) by formula (5) then their right 5% quantiles 𝐿^𝑟₁ and 𝐿^𝑟₂. The computation results are shown in Tables 2 and 3.

Table 2 

The right 5% quantiles of 𝑁₁(𝑐^∗_𝑗,𝑟) (1≤𝑗≤500).

Table 3 

The right 5% quantiles of 𝑁₂(𝑐^∗_𝑗,𝑟) (1≤𝑗≤500).

From the data of Tables 2 and 3, we can draw the corresponding curves as shown in Figures 3 and 4.

Figure 3 

Curves of (𝑟𝑖,𝑁1(𝑐,𝑟𝑖)) and (𝑟𝑖,𝐿𝑟𝑖1) for different value of 𝑟.

Figure 4 

Curves of (𝑟𝑖,𝑁2(𝑐,𝑟𝑖)) and (𝑟𝑖,𝐿𝑟𝑖2) for different value of 𝑟.

From Figures 3 and 4 we can read 𝑟₀=0.44 and 𝑡₀=−0.49.

Therefore, 𝑟^∗=max{0.44,0.49}=0.49≈0.5.

There are in total 893 genes satisfying |𝑃(𝑔_𝑖,𝑐)|≥0.5.

By a two-type discrimination analysis for the tissue samples of the 38 leukemia patients using the 893 gene expression levels as variables, we can build a discrimination function. The 38 patients were identified and classified by using the discrimination function (discrimination function retrospective assessment). The miscarriage of justice was 0.

We have established a prospective evaluation of the discrimination function.

The data are taken from the website of the Broad Institute of MIT [6]. There are total 34 cases of leukemia patients (of which 20 cases of ALL and 14 cases of AML). The 893 gene data were substituted into the discrimination function and classified by type. The miscarriage of justice was 0.02.

Based on the above assessment, we believe the discrimination function established by selecting 893 genes with distinguishing capability from the 7129 genes using our method can be a good discrimination function for classifying leukemia patients and it will provide a good reference for the effective treatment.

4. Discussion

In the statistical methods of classification, the stepwise discrimination analysis is mainly used for variable selection. As the number of data in gene microarrays can be very large, the stepwise discrimination analysis module in commonly used statistical software packages such as SAS and SPSS would not function properly. We have tried to filter genes with classification ability for the whole sample of 7129 genes of the 27 cases of lymphoblastic leukemia and the 11 cases acute myeloid leukemia. The computer program crashed when a discriminant analysis or principal component analysis method was applied since the number of genes was too large. Therefore, we cannot perform discriminant analysis or principal component analysis for the data set on a personal computer.

Therefore, on such a large number of gene chip data for screening, using stepwise discrimination analysis to filter genes with classification ability in a personal computer is infeasible.

A common solution to this problem is to divide the large number of gene data into several groups of genes. The genes with classification ability in each group are selected by the stepwise discrimination analysis of the gene expression levels within each group. Finally, these discrimination functions of each group are combined to build a discrimination function for the whole of genes.

However, this method is also inadequate because the links between genes are separated artificially by gene group division. As tumors are diseases with multigene combined effects, separating the links between genes will reduce the classification ability of the final selected genes. It will in turn affect the subsequent analysis of the classification accuracy on new samples and the results are also not easy to explain [7].

In addition, how many groups of genes are to be divided into is also subjective and this will directly affect the final result for screening of the genes.

The principle idea behind our random permutation vectors method is very similar to a statistical approach, called Randomization Test [8], which is widely used in many applications. The application of the method implies that we have to enumerate all possible combinations of the elements in vector 𝑐 and this is often a very difficult task. In the cases of this paper, there are total 38!/(27!×11!)=1203322288 different combinations if we divide the 38 cases of leukemia patients into two groups of 27 and 11 cases, respectively. This huge number of combinations is really a restriction for us to apply the Randomization Test method to our cases. Therefore, we use the Monte Carlo sampling method further to the vector 𝑐 to generate 500 random combinations. These 500 combinations are 500 random samples of all possible combinations of elements in vector 𝑐. Although the results of 500 samples do not produce an exact answer, it can be close to the exact answer [9]. In order to make the results closer to the exact answer, we may increase the number of random samples. For example, we may increase the number of random samples to 1000 or more in our cases.

Our presented method applies to selecting the genes with the ability of classifying types of diseases from a large number of genes. Compared to the stepwise discrimination analysis by groups, the new method has an obvious advantage of the full usage of information of all genes.

The new method has a low computational complexity and is very practical in practice. The main costs of computation are in the correlation coefficient computation when we need a random permutation vector from the vector 𝑐 and this is not a difficult task for common personal computers.

In the real analysis circumstance, if there are too many genes with their absolute value greater than 𝑟^∗, then the value of 𝑟^∗ can be adjusted to 𝑟^∗+𝑎. The value of 𝑎 can be adjusted according to the actual situation. The feasibility of the adjusted value of 𝑎 can be checked by a retrospective assessment of the discrimination function established on the selected genes.

Conflict of Interests

The authors have declared that no conflict of interests exists.

Acknowledgments

The authors would like to thank the anonymous referees for their many constructive comments and suggestions for enhancing the quality of the paper. The work is made available under the Creative Commons CC0 public domain dedication. This work was partially funded by the Foundation of Science and Technology of Fengze under Grants nos. 2009FZ24 and 2010FZ02 and the Haixi Project of Fujian under Grant no. A099.