# The Classification ofSini Decoction Pattern in Traditional Chinese Medicine by Gene Expression Profiling

Authors: Hung-Tsu Cheng; Chaang-Ray Chen; Chia-Yang Li; Chao-Ying Huang; Wun-Yi Shu; Ian C. Hsu
Journal: Evidence-Based Complementary and Alternative Medicine
 (2016)

---

## Abstract

We investigated the syndromes of theSini decoction pattern (SDP), a common ZHENG in traditional Chinese medicine (TCM). The syndromes of SDP were correlated with various severeYang deficiency related symptoms. To obtain a common profile for SDP, we distributed questionnaires to 300 senior clinical TCM practitioners. According to the survey, we concluded 2 sets of symptoms for SDP: (1) pulse feels deep or faint and (2) reversal cold of the extremities. Twenty-four individuals from Taipei City Hospital, Linsen Chinese Medicine Branch, Taiwan, were recruited. We extracted the total mRNA of peripheral blood mononuclear cells from the 24 individuals for microarray experiments. Twelve individuals (including 6 SDP patients and 6 non-SDP individuals) were used as the training set to identify biomarkers for distinguishing the SDP and non-SDP groups. The remaining 12 individuals were used as the test set. The test results indicated that the gene expression profiles of the identified biomarkers could effectively distinguish the 2 groups by adopting a hierarchical clustering algorithm. Our results suggest the feasibility of using the identified biomarkers in facilitating the diagnosis of TCM ZHENGs. Furthermore, the gene expression profiles of biomarker genes could provide a molecular explanation corresponding to the ZHENG of TCM.


---

## Body

## 1. Introduction

Traditional Chinese medicine (TCM) has been adopted to treat the Chinese for thousands of years on the basis of some classic Chinese medicine textbooks, includingHuang Ti Nei Jing (“The Inner Cannon of Huangti”) [1],Shénnóng Běn Cǎo Jīng [2], andShang Han Za Bing Lun (“Treatise on Cold Damage Disorders”) [3].Shang Han Za Bing Lun is a Chinese medical treatise written by Zhang Zhongjing before AD 220, at the end of the Han dynasty. The treatise is the oldest comprehensive clinical textbook in the world containing principles, methods, formulas, and medicine and is the clinical literature of TCM theory with practice. To obtain the appropriate decoction for a particular ZHENG requires a diagnostic system, namely, the Decoction ZHENG proposed inShang Han Za Bing Lun. Through the methods of inspection, listening, smelling, inquiry, and palpation, physicians collect information on patients to identify their TCM ZHENGs according to their personal experience and then treat these patients with herbs, which is called decoction. The methods applied in TCM have been criticized as insufficiently scientific. To overcome this weakness, we conducted systems biology research to differentiate TCM ZHENGs.Systems biology-based diagnostic principles can be used to support the relationship between TCM and current biomedicine [4, 5]. From the perspective of gene expression profile, TCM ZHENGs differentiation is closely related to gene polymorphisms and differences in gene expression profile. Therefore, applying advanced sequencing techniques and cDNA microarray studies can help to clarify the biological basis of TCM ZHENGs [4]. For example, by using microarray, real-time polymerase chain reaction, and enzyme-linked immunosorbent assay technologies,Shen Yang Xu ZHENG has been determined to be involved with the gene level [6] and may be primarily attributed to the insufficient activity of the mitogen-activated protein kinase (MAPK) pathway [7]. Thecold ZHENG was discovered to be possibly caused by a physiological imbalance and disorder of metabolite processes by using microarray technology [8]. Human mRNAs in peripheral blood mononuclear cells (PBMCs) were extracted for microarray experiments to analyze differently expressed gene expressions for TCMcold andhot ZHENG differentiation [9–11].In this study, we investigated the syndromes ofSini decoction pattern (SDP).Sini decoction was first mentioned inShang Han Za Bing Lun Sini decoction consisting of three herbs: freshAconitum carmichaelii, dryZingiber officinale, andGlycyrrhiza uralensis, which was characterized as a remedy to have essential effect of recuperating the patients from collapse. It was used to treat the syndrome of displaying coldness on the extremities, pulse feels deep or faint, continuous diarrhea with undigested cereal, profuse perspiration, abdominal fullness and distention, vomiting, and lethargy [12]. Based on the clinical applications ofSini decoction, SDP was then developed by Sun [13]. The syndromes of SDP were correlated with numerous symptoms, such as a slow pulse, coldness on the extremities, continuous diarrhea with undigested cereal, and severe diarrhea, but no common profile exists. According to our survey of 300 senior clinical TCM practitioners, we obtained a common profile for SDP: (1) pulse feels deep or faint and (2) reversal cold of the extremities. We collected SDP patients with these 2 symptoms. In this study, 12 patients with SDP and 12 non-SDP individuals were recruited from Taipei City Hospital, Linsen Chinese Medicine Branch, Taiwan.We extracted the total mRNA of PBMCs from the 24 individuals for microarray experiments. Twelve individuals (including 6 SDP patients and 6 non-SDP individuals) were used as the training set to identify the biomarkers that were screened using the hypothesis test for distinguishing the SDP and non-SDP groups. The remaining 12 individuals (including 6 SDP patients and 6 non-SDP individuals) were used as the test set to examine whether the identified biomarkers can distinguish the SDP and non-SDP groups. The test results indicated that the gene expression profiling by the biomarkers could effectively distinguish the 2 groups by using a hierarchical clustering algorithm.Our results suggested the feasibility of using biomarkers extracted from the gene expressions of PBMCs in facilitating the diagnosis of TCM ZHENGs. Furthermore, the gene expression profiles of biomarkers could provide the possible molecular mechanisms of related TCM ZHENGs.


## 2. Materials and Methods

### 2.1. Chinese Medicine Terminology

All Chinese medicine terms used in this paper were in accordance with theInternational Acupuncture Nomenclature (IAN) proposed by WHO in 1991; if not available in IAN, the standard nomenclature proposed in theWHO International Standard Terminologies on Traditional Medicine in the Western Pacific Region [14] was adapted. These translated Chinese medicine terms are presented in italics in this paper.


### 2.2. Human Participants

The syndromes of SDP were associated with a slow, deep, and faint pulse that disappeared, coldness on the extremities, severe diarrhea, continuous diarrhea with undigested cereal, generalized pain, abdominal fullness, cold-fluid retention on the diaphragm, profuse sweating, severe spasms in the feet, tight clamping of the extremities, and reversal cold in both the hands and feet; however, no common profile exists. To obtain a common profile for SDP, we distributed questionnaires to 300 senior clinical TCM practitioners. Our questionnaire listed 41 symptoms related to SDP. In the survey, each senior clinical TCM practitioner assigned a score to each symptom, indicating the correlation between the symptom and SDP from high to low (5, 4, 3, 2, 1, and 0). We identified symptoms that collectively constituted a common SDP profile according to the frequency with which participating TCM practitioners assigned a score of 5 to the symptoms; the frequency threshold was set at higher than 35% of the returned surveys. We recruited SDP patients with the aforementioned symptoms. In this study, 12 SDP patients and 12 non-SDP individuals from Taipei City Hospital, Linsen Chinese Medicine Branch, Taiwan, were recruited. All of the individuals provided informed consent. The definition of healthy subject was diagnosed by western medicine rather different from TCM. It is difficult to diagnose healthy subject. Therefore, we consider substituting non-SDP patients for healthy subject. We try our best to collect non-SDP patients who are alike healthy subjects.


### 2.3. Blood Sample Collection

We phlebotomized 30 mL of human peripheral blood from each donor and then separated the blood into three 10 mL plastic whole blood tubes with spray-coated K2EDTA (Becton Dickenson, Franklin Lakes, NJ, USA). The blood was transferred to a 50 mL centrifuge tube. The PBMCs were then prepared from the 30 mL of peripheral blood through Ficoll-Hypaque density gradient centrifugation. The PBMCs were added to an RNAlater RNA Stabilization Reagent (Qiagen, Valencia, CA, USA) and then placed into a bottle with liquid nitrogen. All of these procedures were performed immediately after the blood was collected at the laboratory of Taipei City Hospital, Linsen Chinese Medicine Branch, Taiwan. The isolated PBMCs were frozen, stored in a liquid nitrogen freezer, and transported to the microarray experimental laboratory at National Tsing Hua University, Hsinchu, Taiwan, for further microarray experiments.


### 2.4. RNA Extraction and Microarray Hybridization

We extracted total RNA by using the RNeasy Mini Kit according to the manufacturer’s protocol (Qiagen, Valencia, CA, USA). Total RNA quantity and quality were assessed using an ultraviolet absorption spectrophotometer at 260 and 280 nm. Subsequently,α- and β-globin mRNA were reduced from a portion of the total RNA samples by using the GLOBINclear™ human kit (Ambion, Austin, TX, USA), according to the manufacturer’s instructions, at the recommended start quantity of 10 μg of total RNA. The quality of total RNA was evaluated using the Agilent 2100 bioanalyzer with the RNA 6000 Nano LabChip kit (Agilent Technologies, Palo Alto, CA, USA). Before reverse transcription, each RNA sample was spiked with a mixture of Arabidopsis mRNAs, which we called doping control. Fluorescence-labeled (Cy3, Cy5-3DNA Capture Reagent) cDNA was conducted using the 3DNA Array 900 labeling kit according to the manufacturer’s protocols (Genisphere, Hatfield, PA, USA). Reverse transcription was performed using SuperScript II (Invitrogen, Carlsbad, CA, USA). Hybridization was performed at 65°C in a water bath for 16 to 18 h, and arrays were washed according to the manufacturer’s protocol (Corning Life Sciences, New York, NY, USA). The arrays were scanned using GenePix 4000B scanner (Axon Instruments, Foster City, CA, USA).


### 2.5. Microarray Fabrication

A total of 9600 human cDNA clones were purchased from Incyte Genomics (Wilmington, DE, USA). These cDNAs were randomly selected from the IMAGE library. Resequencing was performed on these clones, and 7334 cDNAs were verified to have corrected DNA sequences. These 7334 sequence-verified human cDNAs (Incyte Genomics) and 10 Arabidopsis cDNAs (SpotReport™ cDNA Array Validation System; Stratagene, La Jolla, CA, USA), serving as spike-in controls, and one housekeeping gene (β-actin), serving as a positive control, were arrayed on Corning UltraGAPS slides (Corning, Corning, NY, USA). These 7334 human cDNAs and 96 Arabidopsis cDNAs and housekeeping genes were spotted in quadruplicate on each array making a total of 32248 spots to enhance the statistical confidence in the gene expression data. Each array had 32,448 spots. The arrays were postprocessed according to the Corning UltraGAPS Coated Slides instruction manual.


### 2.6. Experimental Design

Loop design, a type of statistical microarray experimental design, was adopted to construct an optimal design. The basic principles of optimal design are a balance among the factors, approximately equal sampling of varieties, and minimal distance between pairs of varieties [15]. Samples from 24 individuals were hybridized on 48 arrays based on the 2 loop designs presented in Figure 1. The individuals of samples L1A1, L1A2, L1A3, L1C1, L1C2, L1C3, L2A1, L2A2, L2A3, L2C1, L2C2, and L2C3 were SDP patients, and those of samples L1B1, L1B2, L1B3, L1D1, L1D2, L1D3, L2B1, L2B2, L2B3, L2D1, L2D2, and L2D3 were non-SDP individuals. Each arrow in the figure indicates a microarray hybridization experiment. The arrow heads and tails represent samples labeled with Cy5 (3DNA Capture Reagent) and Cy3 (3DNA Capture Reagent), respectively. The black solid circles stand for the SDP samples and the open circles stand for non-SDP samples. We used 12 individuals in Loop 1 experiments as the training set and the remaining 12 individuals in Loop 2 experiments as the test set.Figure 1
Experimental designs of Loop 1 (for the training set) and Loop 2 (for the test set). Each arrow indicates a microarray hybridization experiment. The arrow heads and tails represent samples labeled with Cy5 and Cy3, respectively. The black solid circles stand for the SDP samples and the open circles stand for non-SDP samples. Loop 1 is the experimental design for the training set of 12 individuals. Loop 2 is the experimental design for the test set of the other 12 individuals. The annotations of samples in Loop 1 and Loop 2 are shown below each loop.


### 2.7. Microarray Data Analysis and Statistical Analysis

Microarray data preprocessing, normalization, and statistical analysis were performed using the bioinformatics software suite Tsing Hua Engine for Microarray Experiment [16].Spot-screening rules were applied to screen invalid spots on these arrays. The spot-screening rules for the data of Loop 1 (the training set) were as follows: (1) exclude the spots defined as “flag bad” or “absent” in all GenePix Results (GPR), (2) exclude the spots with diameters less than 100μm, (3) exclude the spots with a coefficient of variation (CV) of pixel intensity of over 100% in both channels, and (4) exclude the spots with signal to noise ratios (SNRs) in both channels were less than 2 in the training set microarray experiments. The SNR is defined as (
S
-
B
)
/
B (S: mean of pixel intensities of the signal; B: median pixel intensity of the background). The logarithm of the ratios for all valid spots on each array was normalized using global lowess normalization [17]. After the data were preprocessed, the normalized log ratios of the fluorescence intensity of each cDNA spot were analyzed using a log linear model, which was described in our previous study [18], to obtain the least-squares estimates λ
1
^, λ
2
^, λ
3
^, λ
4
^, λ
5
^, λ
6
^, λ
7
^, λ
8
^, λ
9
^, λ
10
^, λ
11
^, and λ
12
^ for the relative expressions among all of the samples. The relative expressions of sample i are defined as (expression of sample i)/(geometric average of the expression of all samples), for each gene, where i runs from 1 to 12. For the training set microarray data, 4951 genes satisfied the selection criteria of the spot-screening rules and log linear model. The selection criteria for the data of Loop 2 (the test set) are described in the subsequent subsection. The microarray data is available at GEO (GSE67090).


### 2.8. Searching Biomarkers to Distinguish SDP

We used Loop 1 experiments as the training set and Loop 2 experiments as the test set. For the 4951 genes of Loop 1 (the training set), 192 genes were identified using anF test with null hypothesis H
0: (λ
1
^
+
λ
3
^
+
λ
5
^
+
λ
7
^
+
λ
9
^
+
λ
11
^
)
-
(
λ
2
^
+
λ
4
^
+
λ
6
^
+
λ
8
^
+
λ
10
^
+
λ
12
^
)
=
0, at a Bonferroni-adjusted significance level of 0.05/4951, where λ
1
^, λ
3
^, λ
5
^, λ
7
^, λ
9
^, and λ
11
^ were the least-squares estimates for the relative expressions of the SDP samples and λ
2
^, λ
4
^, λ
6
^, λ
8
^, λ
10
^, and λ
12
^ were the least-squares estimates for the relative expressions of the non-SDP samples.The 192 biomarker genes were then used to distinguish the SDP and non-SDP groups of the test set. We applied the selection criteria to filter bad spots to perform an effective normalization and a log linear model. In order to include additional available genes to be tested, we relaxed the selection criteria of the spot-screening rules for the data of Loop 2 (the test set) as follows: (1) exclude the spots defined as “flag bad” or “absent” in all GPR files, (2) exclude the spots with diameters less than 50μm, (3) exclude the spots with a CV of pixel standard deviation of over 150% in both channels, and (4) exclude the spots with SNRs in either channel were less than 1 in the test microarray experiments. The logarithm of the ratios for all valid spots on each array was normalized under the same normalization conditions as those for the training set. We then used hierarchical clustering to verify whether the gene expression profiles of the 192 biomarker genes could classify the 12 test samples into 2 groups and how well they could distinguish the SDP patients and non-SDP individuals.


### 2.9. Functional Annotation and Enrichment Analysis

Functional annotation and enrichment analysis were accessed through the database for annotation, visualization, and integrated discovery (DAVID) [19], which performed a modified Fisher’s exact test to identify overrepresented functions. DAVID is a widely employed web tool that provides a rapid means to understand the biological meaning of a long list of genes of interest. For any gene list, DAVID can map and convert identifiers (IDs) and accessions among different gene IDs. A long gene list can be reduced into biological theme-related groups of genes, for example, gene ontology (GO) [20] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway [21].The Entrez gene IDs of 192 biomarker genes were used as the gene ID for bioinformatics analysis. Biological processes of GO terms were used for functional annotation and enrichment analysis of these biomarker genes.


## 2.1. Chinese Medicine Terminology

All Chinese medicine terms used in this paper were in accordance with theInternational Acupuncture Nomenclature (IAN) proposed by WHO in 1991; if not available in IAN, the standard nomenclature proposed in theWHO International Standard Terminologies on Traditional Medicine in the Western Pacific Region [14] was adapted. These translated Chinese medicine terms are presented in italics in this paper.


## 2.2. Human Participants

The syndromes of SDP were associated with a slow, deep, and faint pulse that disappeared, coldness on the extremities, severe diarrhea, continuous diarrhea with undigested cereal, generalized pain, abdominal fullness, cold-fluid retention on the diaphragm, profuse sweating, severe spasms in the feet, tight clamping of the extremities, and reversal cold in both the hands and feet; however, no common profile exists. To obtain a common profile for SDP, we distributed questionnaires to 300 senior clinical TCM practitioners. Our questionnaire listed 41 symptoms related to SDP. In the survey, each senior clinical TCM practitioner assigned a score to each symptom, indicating the correlation between the symptom and SDP from high to low (5, 4, 3, 2, 1, and 0). We identified symptoms that collectively constituted a common SDP profile according to the frequency with which participating TCM practitioners assigned a score of 5 to the symptoms; the frequency threshold was set at higher than 35% of the returned surveys. We recruited SDP patients with the aforementioned symptoms. In this study, 12 SDP patients and 12 non-SDP individuals from Taipei City Hospital, Linsen Chinese Medicine Branch, Taiwan, were recruited. All of the individuals provided informed consent. The definition of healthy subject was diagnosed by western medicine rather different from TCM. It is difficult to diagnose healthy subject. Therefore, we consider substituting non-SDP patients for healthy subject. We try our best to collect non-SDP patients who are alike healthy subjects.


## 2.3. Blood Sample Collection

We phlebotomized 30 mL of human peripheral blood from each donor and then separated the blood into three 10 mL plastic whole blood tubes with spray-coated K2EDTA (Becton Dickenson, Franklin Lakes, NJ, USA). The blood was transferred to a 50 mL centrifuge tube. The PBMCs were then prepared from the 30 mL of peripheral blood through Ficoll-Hypaque density gradient centrifugation. The PBMCs were added to an RNAlater RNA Stabilization Reagent (Qiagen, Valencia, CA, USA) and then placed into a bottle with liquid nitrogen. All of these procedures were performed immediately after the blood was collected at the laboratory of Taipei City Hospital, Linsen Chinese Medicine Branch, Taiwan. The isolated PBMCs were frozen, stored in a liquid nitrogen freezer, and transported to the microarray experimental laboratory at National Tsing Hua University, Hsinchu, Taiwan, for further microarray experiments.


## 2.4. RNA Extraction and Microarray Hybridization

We extracted total RNA by using the RNeasy Mini Kit according to the manufacturer’s protocol (Qiagen, Valencia, CA, USA). Total RNA quantity and quality were assessed using an ultraviolet absorption spectrophotometer at 260 and 280 nm. Subsequently,α- and β-globin mRNA were reduced from a portion of the total RNA samples by using the GLOBINclear™ human kit (Ambion, Austin, TX, USA), according to the manufacturer’s instructions, at the recommended start quantity of 10 μg of total RNA. The quality of total RNA was evaluated using the Agilent 2100 bioanalyzer with the RNA 6000 Nano LabChip kit (Agilent Technologies, Palo Alto, CA, USA). Before reverse transcription, each RNA sample was spiked with a mixture of Arabidopsis mRNAs, which we called doping control. Fluorescence-labeled (Cy3, Cy5-3DNA Capture Reagent) cDNA was conducted using the 3DNA Array 900 labeling kit according to the manufacturer’s protocols (Genisphere, Hatfield, PA, USA). Reverse transcription was performed using SuperScript II (Invitrogen, Carlsbad, CA, USA). Hybridization was performed at 65°C in a water bath for 16 to 18 h, and arrays were washed according to the manufacturer’s protocol (Corning Life Sciences, New York, NY, USA). The arrays were scanned using GenePix 4000B scanner (Axon Instruments, Foster City, CA, USA).


## 2.5. Microarray Fabrication

A total of 9600 human cDNA clones were purchased from Incyte Genomics (Wilmington, DE, USA). These cDNAs were randomly selected from the IMAGE library. Resequencing was performed on these clones, and 7334 cDNAs were verified to have corrected DNA sequences. These 7334 sequence-verified human cDNAs (Incyte Genomics) and 10 Arabidopsis cDNAs (SpotReport™ cDNA Array Validation System; Stratagene, La Jolla, CA, USA), serving as spike-in controls, and one housekeeping gene (β-actin), serving as a positive control, were arrayed on Corning UltraGAPS slides (Corning, Corning, NY, USA). These 7334 human cDNAs and 96 Arabidopsis cDNAs and housekeeping genes were spotted in quadruplicate on each array making a total of 32248 spots to enhance the statistical confidence in the gene expression data. Each array had 32,448 spots. The arrays were postprocessed according to the Corning UltraGAPS Coated Slides instruction manual.


## 2.6. Experimental Design

Loop design, a type of statistical microarray experimental design, was adopted to construct an optimal design. The basic principles of optimal design are a balance among the factors, approximately equal sampling of varieties, and minimal distance between pairs of varieties [15]. Samples from 24 individuals were hybridized on 48 arrays based on the 2 loop designs presented in Figure 1. The individuals of samples L1A1, L1A2, L1A3, L1C1, L1C2, L1C3, L2A1, L2A2, L2A3, L2C1, L2C2, and L2C3 were SDP patients, and those of samples L1B1, L1B2, L1B3, L1D1, L1D2, L1D3, L2B1, L2B2, L2B3, L2D1, L2D2, and L2D3 were non-SDP individuals. Each arrow in the figure indicates a microarray hybridization experiment. The arrow heads and tails represent samples labeled with Cy5 (3DNA Capture Reagent) and Cy3 (3DNA Capture Reagent), respectively. The black solid circles stand for the SDP samples and the open circles stand for non-SDP samples. We used 12 individuals in Loop 1 experiments as the training set and the remaining 12 individuals in Loop 2 experiments as the test set.Figure 1
Experimental designs of Loop 1 (for the training set) and Loop 2 (for the test set). Each arrow indicates a microarray hybridization experiment. The arrow heads and tails represent samples labeled with Cy5 and Cy3, respectively. The black solid circles stand for the SDP samples and the open circles stand for non-SDP samples. Loop 1 is the experimental design for the training set of 12 individuals. Loop 2 is the experimental design for the test set of the other 12 individuals. The annotations of samples in Loop 1 and Loop 2 are shown below each loop.


## 2.7. Microarray Data Analysis and Statistical Analysis

Microarray data preprocessing, normalization, and statistical analysis were performed using the bioinformatics software suite Tsing Hua Engine for Microarray Experiment [16].Spot-screening rules were applied to screen invalid spots on these arrays. The spot-screening rules for the data of Loop 1 (the training set) were as follows: (1) exclude the spots defined as “flag bad” or “absent” in all GenePix Results (GPR), (2) exclude the spots with diameters less than 100μm, (3) exclude the spots with a coefficient of variation (CV) of pixel intensity of over 100% in both channels, and (4) exclude the spots with signal to noise ratios (SNRs) in both channels were less than 2 in the training set microarray experiments. The SNR is defined as (
S
-
B
)
/
B (S: mean of pixel intensities of the signal; B: median pixel intensity of the background). The logarithm of the ratios for all valid spots on each array was normalized using global lowess normalization [17]. After the data were preprocessed, the normalized log ratios of the fluorescence intensity of each cDNA spot were analyzed using a log linear model, which was described in our previous study [18], to obtain the least-squares estimates λ
1
^, λ
2
^, λ
3
^, λ
4
^, λ
5
^, λ
6
^, λ
7
^, λ
8
^, λ
9
^, λ
10
^, λ
11
^, and λ
12
^ for the relative expressions among all of the samples. The relative expressions of sample i are defined as (expression of sample i)/(geometric average of the expression of all samples), for each gene, where i runs from 1 to 12. For the training set microarray data, 4951 genes satisfied the selection criteria of the spot-screening rules and log linear model. The selection criteria for the data of Loop 2 (the test set) are described in the subsequent subsection. The microarray data is available at GEO (GSE67090).


## 2.8. Searching Biomarkers to Distinguish SDP

We used Loop 1 experiments as the training set and Loop 2 experiments as the test set. For the 4951 genes of Loop 1 (the training set), 192 genes were identified using anF test with null hypothesis H
0: (λ
1
^
+
λ
3
^
+
λ
5
^
+
λ
7
^
+
λ
9
^
+
λ
11
^
)
-
(
λ
2
^
+
λ
4
^
+
λ
6
^
+
λ
8
^
+
λ
10
^
+
λ
12
^
)
=
0, at a Bonferroni-adjusted significance level of 0.05/4951, where λ
1
^, λ
3
^, λ
5
^, λ
7
^, λ
9
^, and λ
11
^ were the least-squares estimates for the relative expressions of the SDP samples and λ
2
^, λ
4
^, λ
6
^, λ
8
^, λ
10
^, and λ
12
^ were the least-squares estimates for the relative expressions of the non-SDP samples.The 192 biomarker genes were then used to distinguish the SDP and non-SDP groups of the test set. We applied the selection criteria to filter bad spots to perform an effective normalization and a log linear model. In order to include additional available genes to be tested, we relaxed the selection criteria of the spot-screening rules for the data of Loop 2 (the test set) as follows: (1) exclude the spots defined as “flag bad” or “absent” in all GPR files, (2) exclude the spots with diameters less than 50μm, (3) exclude the spots with a CV of pixel standard deviation of over 150% in both channels, and (4) exclude the spots with SNRs in either channel were less than 1 in the test microarray experiments. The logarithm of the ratios for all valid spots on each array was normalized under the same normalization conditions as those for the training set. We then used hierarchical clustering to verify whether the gene expression profiles of the 192 biomarker genes could classify the 12 test samples into 2 groups and how well they could distinguish the SDP patients and non-SDP individuals.


## 2.9. Functional Annotation and Enrichment Analysis

Functional annotation and enrichment analysis were accessed through the database for annotation, visualization, and integrated discovery (DAVID) [19], which performed a modified Fisher’s exact test to identify overrepresented functions. DAVID is a widely employed web tool that provides a rapid means to understand the biological meaning of a long list of genes of interest. For any gene list, DAVID can map and convert identifiers (IDs) and accessions among different gene IDs. A long gene list can be reduced into biological theme-related groups of genes, for example, gene ontology (GO) [20] and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway [21].The Entrez gene IDs of 192 biomarker genes were used as the gene ID for bioinformatics analysis. Biological processes of GO terms were used for functional annotation and enrichment analysis of these biomarker genes.


## 3. Results

### 3.1. Questionnaires to Obtain a Common Profile for SDP

We distributed questionnaires to 300 senior clinical TCM practitioners, and 43 questionnaires were returned. Our questionnaire listed 41 symptoms related to SDP. In the survey, each senior clinical TCM practitioner assigned a score to each symptom, indicating the correlation between the symptom and SDP, from high to low (5, 4, 3, 2, 1, and 0). Based on the frequency with which participating TCM practitioners assigned a score of 5 to the symptoms, we identified symptoms that collectively constituted a common SDP profile (Figure2). The frequency threshold was set at higher than 35% of the 43 returned surveys (i.e., symptoms assigned a score of 5 by at least 16 surveys). This process revealed 7 symptoms to form our symptom profile of SDP patients. By grouping related symptoms, we obtained a common profile for SDP as 2 sets of symptoms: (1) pulse feels deep or faint and (2) reversal cold of the extremities. We recruited SDP patients with the 2 sets of symptoms. In this study, 12 SDP patients and 12 non-SDP individuals were recruited from Taipei City Hospital, Linsen Chinese Medicine Branch, Taiwan. All of the participants provided informed consent.Figure 2
Frequency with which survey symptoms were assigned a score of 5 in determining the SDP profile.


### 3.2. Identifying Biomarker Genes from the Training Set Microarray Experiment

On the basis of the results for the 4951 genes of Loop 1 (the training set), 192 genes were identified (Table1) using an F test with null hypothesis H
0: (
λ
1
^
+
λ
3
^
+
λ
5
^
+
λ
7
^
+
λ
9
^
+
λ
11
^
)
-
(
λ
2
^
+
λ
4
^
+
λ
6
^
+
λ
8
^
+
λ
10
^
+
λ
12
^
)
=
0, at a Bonferroni-adjusted significance level of 0.05/4951. By applying hierarchical clustering with Spearman distance and average linkage, the 12 training samples were appropriately separated into 2 groups, as shown in Figure 3. In Figure 3, the gene expression value used to do the clustering is the relative expressions of the labeled sample.Table 1
The 192 biomarker genes identified using anF test with null hypothesis H
0: (
λ
1
^
+
λ
3
^
+
λ
5
^
+
λ
7
^
+
λ
9
^
+
λ
11
^
)
-
(
λ
2
^
+
λ
4
^
+
λ
6
^
+
λ
8
^
+
λ
10
^
+
λ
12
^
)
=
0, at a Bonferroni-adjusted significance level of 0.05/4951.
Gene name
ABHD3
ABI2
ACSM2B
ACTB
AGFG1
AKAP9
ARG1
ARID5A
ATP1B3
ATP6AP2
ATP6V0B
BLZF1
BPGM
C11orf79
C12orf24
C12orf32
C17orf91
C6orf106
CAPN1
CBARA1
CCL18
CCL5
CD1D
CD24
CD40
CD44
CD83
CD8A
CD93
CDC42
CDC42SE2
CLDN1
COMT
COX10
CP
CREM
CST7
CXCL9
DEFA1
DEFA3
DEFA4
DENND2D
DGCR6L
DLG2
DUSP6
EEF2
EIF3E
EIF3I
EPAS1
EXOSC9
FAM113B
FAM175B
FCER1A
FLNB
FMOD
FOLR3
FOS
FOSB
FYN
GCA
GIMAP1
GLRX
GNG11
GNG2
GNLY
GPR56
GSTA2
GTF2E1
GZMA
H3F3B
HERPUD1
HIST1H2AC
HIST1H2BD
HIST1H2BG
HIST1H2BJ
HIST2H2AA3
HLA-A
HLA-DMA
HLA-DQB1
HLA-DRB1
HLA-G
Hs.282050
Hs.467411
Hs.534061
Hs.586682
Hs.591807
Hs.599359
Hs.655903
Hs.674438
Hs.88605
HSP90B1
HSPA1A
HTATIP2
HTR7P1
ID2
IFIT1
IL10RA
IL15RA
IL1B
IL32
JARID1C
KIF3C
KLF10
KLF6
LCP1
LNP1
LOC400027
LOC729776
LTB
LYRM1
MAPK13
MBP
MCM3AP
MEST
MLH3
MNDA
MRPL3
MRPL49
MTMR6
MYL6
NCRNA00084
NHEDC2
NPY
NUCB1
NUDT16
NUP93
ORM1
PACS1
PCDH8
PDE4B
PEA15
PF4
PLAUR
PNLIPRP2
POLR1C
POLR2D
POU2F2
PPBP
PPP1R12B
PPP1R7
PPP2R2A
PPRC1
PPT1
PRL
PSMC6
RAB4A
RARRES3
RGS18
RGS2
RHOQ
RNF40
RNPEP
RPL41
RPS11
RPS24
RPS27
RPS4X
RPS4Y1
RUNX3
S1PR5
SEC22B
SELL
SERPINE2
SFRS1
SH2B3
SH3D19
SHCBP1
SLC2A3P1
SOD2
SPTY2D1
SRGN
STX8
SYTL2
TADA3L
TARP
THBS1
TIMP1
TMEM176A
TNFAIP8L2
TPCN1
TRA@
TRAPPC1
TRMT11
VIM
VPS13D
WDR42A
WIZ
XPNPEP1
YARS
Unknown∗
Unknown∗
ZNF286A
∗The sequence of the clone was not yet identified in the BLAST.Figure 3
192 genes, identified using theF test at a Bonferroni-adjusted significance level of 0.05/4951, were clustered into 2 groups. By applying hierarchical clustering with Spearman distance and average linkage, the 12 training samples were appropriately separated into 2 groups according to the 192 marker genes. Samples L1C2, L1A3, L1C3, L1A1, L1C1, and L1A2 comprised the SDP group, and samples L1D2, L1B3, L1D3, L1D1, L1B1, and L1B2 comprised the non-SDP group. In this figure, the gene expression value used to do the clustering is the relative expressions of the labeled sample.


### 3.3. Distinguishing Power of the Biomarker Genes to the Test Set Microarray Experiment

The 192 biomarker genes identified in the training set experiments were then used to test the ability to distinguish the SDP patients and non-SDP individuals in the test set. To include additional available genes to be tested, we relaxed the selection criteria of the spot-screening rules for the data of the test set (Loop 2). The logarithm of the ratios for all valid spots on each array was normalized under the same normalization conditions as those for the training set. We then used hierarchical clustering to verify whether the gene expression profiles of the 192 biomarker genes could classify the 12 test samples into 2 groups and how well they could distinguish the SDP patients and non-SDP individuals. Under the same clustering conditions (i.e., hierarchical clustering with Spearman distance and average linkage), Figure4 showed that the 192 biomarker genes can classify the 12 test samples into 2 groups. One group comprised L2A2, L2A1, L2A3, L2C3, L2C2, and L2D1 samples and the other comprised L2D2, L2B3, L2D3, L2B2, L2C1, and L2B1 samples. Except for L2D1 and L2C1 samples, all of the samples were correctly grouped into SDP and non-SDP groups. In Figure 4, the gene expression value used to do the clustering is the relative expressions of the labeled sample. The microarray data of 4 biomarker genes from 7 samples were confirmed by quantitative PCR (qPCR). The correlation coefficients r between microarray data and qPCR for each gene are 0.73, 0.83, 0.93, and 0.76 (Appendix A).Figure 4
Hierarchical clustering results for the 12 test samples. Under the same clustering conditions (i.e., hierarchical clustering with Spearman distance and average linkage), the gene expression profiles of the 192 marker genes classified the 12 test samples into 2 groups. One group consisted of the L2A2, L2A1, L2A3, L2C3, L2C2, and L2D1 samples and the other consisted of the L2D2, L2B3, L2D3, L2B2, L2C1, and L2B1 samples. Except for the L2D1 and L2C1 samples, all of the samples were correctly grouped into SDP and non-SDP groups. In this figure, the gene expression values that used to do the clustering are the relative expressions of the labeled sample.


### 3.4. Functional Annotation and Enrichment Analysis Results for the 192 Biomarker Genes

The functional annotation and enrichment analysis were accessed using DAVID. Seventy-eight GO terms were enriched with aP value less than 0.05. All enriched GO terms are listed in Appendix B. 78 GO terms were enriched with a P value of less than 0.05. We further applied 2 stringent criteria to filter the high-ranked GO terms: (1) GO terms containing more than 10 genes and (2) a P value of less than 0.01. Consequently, 11 GO terms were identified (Table 2).Table 2
The 11 enriched GO terms for 192 biomarker genes (sorted according toP value).
GO term
Description
# gene
P value
GO:0006952
Defense response
24
0.000000
GO:0006955
Immune response
24
0.000002
GO:0009617
Response to bacterium
11
0.000088
GO:0001817
Regulation of cytokine production
11
0.000172
GO:0048534
Hemopoietic or lymphoid organ development
13
0.000562
GO:0030097
Hemopoiesis
12
0.000999
GO:0002520
Immune system development
13
0.001183
GO:0007610
Behavior
14
0.004162
GO:0006412
Translation
12
0.004210
GO:0001775
Cell activation
11
0.008090
GO:0019725
Cellular homeostasis
13
0.009574Table2 showed that the 11 GO terms are substantially correlated with immune, inflammation, and hemopoiesis functions. In detail, the immune function includes defense response, immune response, response to bacterium, and immune system development. The inflammation function includes cytokine production regulation, behavior, cell activation, and cellular homeostasis. The hematopoiesis function includes hemopoietic or lymphoid organ development and hemopoiesis.


## 3.1. Questionnaires to Obtain a Common Profile for SDP

We distributed questionnaires to 300 senior clinical TCM practitioners, and 43 questionnaires were returned. Our questionnaire listed 41 symptoms related to SDP. In the survey, each senior clinical TCM practitioner assigned a score to each symptom, indicating the correlation between the symptom and SDP, from high to low (5, 4, 3, 2, 1, and 0). Based on the frequency with which participating TCM practitioners assigned a score of 5 to the symptoms, we identified symptoms that collectively constituted a common SDP profile (Figure2). The frequency threshold was set at higher than 35% of the 43 returned surveys (i.e., symptoms assigned a score of 5 by at least 16 surveys). This process revealed 7 symptoms to form our symptom profile of SDP patients. By grouping related symptoms, we obtained a common profile for SDP as 2 sets of symptoms: (1) pulse feels deep or faint and (2) reversal cold of the extremities. We recruited SDP patients with the 2 sets of symptoms. In this study, 12 SDP patients and 12 non-SDP individuals were recruited from Taipei City Hospital, Linsen Chinese Medicine Branch, Taiwan. All of the participants provided informed consent.Figure 2
Frequency with which survey symptoms were assigned a score of 5 in determining the SDP profile.


## 3.2. Identifying Biomarker Genes from the Training Set Microarray Experiment

On the basis of the results for the 4951 genes of Loop 1 (the training set), 192 genes were identified (Table1) using an F test with null hypothesis H
0: (
λ
1
^
+
λ
3
^
+
λ
5
^
+
λ
7
^
+
λ
9
^
+
λ
11
^
)
-
(
λ
2
^
+
λ
4
^
+
λ
6
^
+
λ
8
^
+
λ
10
^
+
λ
12
^
)
=
0, at a Bonferroni-adjusted significance level of 0.05/4951. By applying hierarchical clustering with Spearman distance and average linkage, the 12 training samples were appropriately separated into 2 groups, as shown in Figure 3. In Figure 3, the gene expression value used to do the clustering is the relative expressions of the labeled sample.Table 1
The 192 biomarker genes identified using anF test with null hypothesis H
0: (
λ
1
^
+
λ
3
^
+
λ
5
^
+
λ
7
^
+
λ
9
^
+
λ
11
^
)
-
(
λ
2
^
+
λ
4
^
+
λ
6
^
+
λ
8
^
+
λ
10
^
+
λ
12
^
)
=
0, at a Bonferroni-adjusted significance level of 0.05/4951.
Gene name
ABHD3
ABI2
ACSM2B
ACTB
AGFG1
AKAP9
ARG1
ARID5A
ATP1B3
ATP6AP2
ATP6V0B
BLZF1
BPGM
C11orf79
C12orf24
C12orf32
C17orf91
C6orf106
CAPN1
CBARA1
CCL18
CCL5
CD1D
CD24
CD40
CD44
CD83
CD8A
CD93
CDC42
CDC42SE2
CLDN1
COMT
COX10
CP
CREM
CST7
CXCL9
DEFA1
DEFA3
DEFA4
DENND2D
DGCR6L
DLG2
DUSP6
EEF2
EIF3E
EIF3I
EPAS1
EXOSC9
FAM113B
FAM175B
FCER1A
FLNB
FMOD
FOLR3
FOS
FOSB
FYN
GCA
GIMAP1
GLRX
GNG11
GNG2
GNLY
GPR56
GSTA2
GTF2E1
GZMA
H3F3B
HERPUD1
HIST1H2AC
HIST1H2BD
HIST1H2BG
HIST1H2BJ
HIST2H2AA3
HLA-A
HLA-DMA
HLA-DQB1
HLA-DRB1
HLA-G
Hs.282050
Hs.467411
Hs.534061
Hs.586682
Hs.591807
Hs.599359
Hs.655903
Hs.674438
Hs.88605
HSP90B1
HSPA1A
HTATIP2
HTR7P1
ID2
IFIT1
IL10RA
IL15RA
IL1B
IL32
JARID1C
KIF3C
KLF10
KLF6
LCP1
LNP1
LOC400027
LOC729776
LTB
LYRM1
MAPK13
MBP
MCM3AP
MEST
MLH3
MNDA
MRPL3
MRPL49
MTMR6
MYL6
NCRNA00084
NHEDC2
NPY
NUCB1
NUDT16
NUP93
ORM1
PACS1
PCDH8
PDE4B
PEA15
PF4
PLAUR
PNLIPRP2
POLR1C
POLR2D
POU2F2
PPBP
PPP1R12B
PPP1R7
PPP2R2A
PPRC1
PPT1
PRL
PSMC6
RAB4A
RARRES3
RGS18
RGS2
RHOQ
RNF40
RNPEP
RPL41
RPS11
RPS24
RPS27
RPS4X
RPS4Y1
RUNX3
S1PR5
SEC22B
SELL
SERPINE2
SFRS1
SH2B3
SH3D19
SHCBP1
SLC2A3P1
SOD2
SPTY2D1
SRGN
STX8
SYTL2
TADA3L
TARP
THBS1
TIMP1
TMEM176A
TNFAIP8L2
TPCN1
TRA@
TRAPPC1
TRMT11
VIM
VPS13D
WDR42A
WIZ
XPNPEP1
YARS
Unknown∗
Unknown∗
ZNF286A
∗The sequence of the clone was not yet identified in the BLAST.Figure 3
192 genes, identified using theF test at a Bonferroni-adjusted significance level of 0.05/4951, were clustered into 2 groups. By applying hierarchical clustering with Spearman distance and average linkage, the 12 training samples were appropriately separated into 2 groups according to the 192 marker genes. Samples L1C2, L1A3, L1C3, L1A1, L1C1, and L1A2 comprised the SDP group, and samples L1D2, L1B3, L1D3, L1D1, L1B1, and L1B2 comprised the non-SDP group. In this figure, the gene expression value used to do the clustering is the relative expressions of the labeled sample.


## 3.3. Distinguishing Power of the Biomarker Genes to the Test Set Microarray Experiment

The 192 biomarker genes identified in the training set experiments were then used to test the ability to distinguish the SDP patients and non-SDP individuals in the test set. To include additional available genes to be tested, we relaxed the selection criteria of the spot-screening rules for the data of the test set (Loop 2). The logarithm of the ratios for all valid spots on each array was normalized under the same normalization conditions as those for the training set. We then used hierarchical clustering to verify whether the gene expression profiles of the 192 biomarker genes could classify the 12 test samples into 2 groups and how well they could distinguish the SDP patients and non-SDP individuals. Under the same clustering conditions (i.e., hierarchical clustering with Spearman distance and average linkage), Figure4 showed that the 192 biomarker genes can classify the 12 test samples into 2 groups. One group comprised L2A2, L2A1, L2A3, L2C3, L2C2, and L2D1 samples and the other comprised L2D2, L2B3, L2D3, L2B2, L2C1, and L2B1 samples. Except for L2D1 and L2C1 samples, all of the samples were correctly grouped into SDP and non-SDP groups. In Figure 4, the gene expression value used to do the clustering is the relative expressions of the labeled sample. The microarray data of 4 biomarker genes from 7 samples were confirmed by quantitative PCR (qPCR). The correlation coefficients r between microarray data and qPCR for each gene are 0.73, 0.83, 0.93, and 0.76 (Appendix A).Figure 4
Hierarchical clustering results for the 12 test samples. Under the same clustering conditions (i.e., hierarchical clustering with Spearman distance and average linkage), the gene expression profiles of the 192 marker genes classified the 12 test samples into 2 groups. One group consisted of the L2A2, L2A1, L2A3, L2C3, L2C2, and L2D1 samples and the other consisted of the L2D2, L2B3, L2D3, L2B2, L2C1, and L2B1 samples. Except for the L2D1 and L2C1 samples, all of the samples were correctly grouped into SDP and non-SDP groups. In this figure, the gene expression values that used to do the clustering are the relative expressions of the labeled sample.


## 3.4. Functional Annotation and Enrichment Analysis Results for the 192 Biomarker Genes

The functional annotation and enrichment analysis were accessed using DAVID. Seventy-eight GO terms were enriched with aP value less than 0.05. All enriched GO terms are listed in Appendix B. 78 GO terms were enriched with a P value of less than 0.05. We further applied 2 stringent criteria to filter the high-ranked GO terms: (1) GO terms containing more than 10 genes and (2) a P value of less than 0.01. Consequently, 11 GO terms were identified (Table 2).Table 2
The 11 enriched GO terms for 192 biomarker genes (sorted according toP value).
GO term
Description
# gene
P value
GO:0006952
Defense response
24
0.000000
GO:0006955
Immune response
24
0.000002
GO:0009617
Response to bacterium
11
0.000088
GO:0001817
Regulation of cytokine production
11
0.000172
GO:0048534
Hemopoietic or lymphoid organ development
13
0.000562
GO:0030097
Hemopoiesis
12
0.000999
GO:0002520
Immune system development
13
0.001183
GO:0007610
Behavior
14
0.004162
GO:0006412
Translation
12
0.004210
GO:0001775
Cell activation
11
0.008090
GO:0019725
Cellular homeostasis
13
0.009574Table2 showed that the 11 GO terms are substantially correlated with immune, inflammation, and hemopoiesis functions. In detail, the immune function includes defense response, immune response, response to bacterium, and immune system development. The inflammation function includes cytokine production regulation, behavior, cell activation, and cellular homeostasis. The hematopoiesis function includes hemopoietic or lymphoid organ development and hemopoiesis.


## 4. Discussion

Of the 4951 genes of the training set, 192 genes were identified using anF test with null hypothesis. By applying hierarchical clustering with Spearman distance and average linkage, the 12 training samples were appropriately separated into SDP and non-SDP groups. The 192 biomarker genes were then used to distinguish the SDP patients and non-SDP individuals in the test set. Under the same clustering conditions, the 12 test samples were clustered into 2 groups. One group comprised L2A2, L2A1, L2A3, L2C3, L2C2, and L2D1 samples and the other comprised L2D2, L2B3, L2D3, L2B2, L2C1, and L2B1 samples. Except for L2D1 and L2C1 samples, all of the samples were correctly grouped into the SDP and non-SDP groups. Thus, we reviewed the original information of these 2 samples. L2C1 sample (classified as SDP according to TCM but clustered into the non-SDP group by using biomarkers) had a deep, rough, and rapid pulse and reversal cold of the extremities. Among the symptoms of L2C1 sample, reversal cold of the extremities matched the selected symptom of SDP; however, a deep, rough, and rapid pulse was inconsistent with a deep or faint pulse, a selected symptom for SDP. L2D1 sample (categorized as non-SDP according to TCM but clustered into the SDP group by using biomarkers) had a deep and tight pulse, which was particularly deep at theguan site (the guan is immediately central to the radial styloid at the wrist, where the tip of the physician’s middle finger is placed), and reversal cold of the extremities. Among the symptoms of L2D1 sample, a deep and tight pulse, particularly at the guan site, did not meet the symptom of a deep or faint pulse, and the sample was thus classified as non-SDP according to TCM. However, reversal cold of the extremities matched the selected symptom of SDP. These findings indicate that, in classifying these 2 samples into SDP and non-SDP groups, uncertain areas exist. Nevertheless, determining biomarkers by using the mRNA profiling method can facilitate identifying the SDP group. Therefore, the traditional dialectical method of TCM is quite trustworthy.The results of the functional annotation and enrichment analysis of 192 biomarker genes revealed that the 11 GO terms were substantially correlated with immune, inflammation, and hemopoiesis functions. Among the symptoms of SDP, a deep or faint pulse implied the insufficiency ofQ
i (in the field of medicine, Q
i refers both to the refined nutritive substance that flows within the human body and to its functional activities). Insufficiency of Q
i indicates the body’s resistance against disease is weak. Pathogens can invade the human body easily, and viruses and bacteria could then cause inflammation easily.


## 5. Conclusion

Our results provided the biomarker genes of SDP, which was first mentioned inShang Han Za Bing Lun. The results provide molecular evidence for TCM in the ZHENGs differentiation of SDP. Using mRNA profiling experiments on patients’ PBMCs to discover biomarker genes to test an unknown group of individuals can effectively identify a specific ZHENG. We hope that the mRNA profiling method will become a model for TCM physicians in diagnosis and assist TCM practitioners in evaluating the curative effect of the medication.



---

*Source: 8239817-2016-04-20.xml*