I.

주성분 분석을 이용한 유전자 집합 기반의 마이크로어레이 자료 분석

Gene Set-Wise Analysis of Microarray Data Using Principal Component Analysis

Sung Bum Cho M.D. M.S.1, Tae Su Chung PhD.1,2, Ju Han Kim M.D., PhD.1,2*

1Seoul National University Biomedical Informatics (SNUBI), 2Human Genome Research Institute, College of Medicine, Seoul National University, Seoul 110-799, Korea, csb1749@snu.ac.kr, *correspondence to juhan@snu.ac.kr

Abstract: Gene set-wise analysis provided many researchers a powerful way to understand microarray data. Correlation between gene set members, however, can confuse the gene set-wise testing procedure. In this analysis, the authors performed principal component analysis(PCA) of each gene set to avoid possible disturbing effect of highly correlated gene set members. The resulting orthogonal property of principal components can remedy the high correlation between gene set members. Global test was used for gene set-wise analysis. The analysis of breast cancer data set(GSE 3494 from GeneExpress Ominibus web site) with 387 pathways revealed that 48 pathways showed significance for lymph node metastasis consistently both in raw and transformed data. To measure the degree of correlation between gene set members, we used correlation index(CI). Among 387 pathways, 38 pathways with higher CI than mean CI value showed non-significant result in global test. PCA transformation of original expression matrix of the pathways revealed that these pathways were significant. These result indicated that PCA might be helpful for concerning the correlation structure in the gene set-wise analysis of microarray

Keyword: gene set-wise microarray data analysis, principal component analysis, correlation index

배 경

마이크로어레이 자료는 한 번의 실험으로 유전체의 전사체(transcriptome) 정보를 제공한다는 커다란 이점을 가지고 있다. 유전자들의 발현량을 한꺼번에 측정하여, 그들 상호간의 관계에 대한 추론이 가능하다는 이점 때문에 많은 생물학 연구자들이 새로운 유전자의 기능을 연구하거나, 특정 질환이나 질병의 기전을 연구하는데 마이크로어레이 실험을 사용하고 있다.

마이크로어레이가 등장한 초기에는 서로 다른 조건하에서 차별 발현하는 유전자(differentially expressed genes, DEGs)를 발굴하는 알고리즘들이 많이 연구되었다.[1] 그러나 이러한 접근법들의 한계가 지적되면서 생물학적인 의미를 지니는 유전자 집합 기반(gene set-wise)의 분석이 개발되었다.[2],[3] 그러나 이러한 분석에 있어서 규정된 유전자 집합 내의 유전자끼리의 발현량이 서로 높은 연관성을 지니는 경우에 잘못된 결론을 낼 수 있다.[4]

주성분 분석(Principal component analysis, PCA)은 대표적인 차원 축소 알고리즘으로 변수들 간의 공분산 구조를 반영하여 각 변수의 일차 결합으로 자료를 설명할 수 있는 성분(component)을 추출할 수 있다. 이 때, 각 성분간에는 직교성(orthogonality)이 유지되므로 연관 구조(correalation structure)를 분해하는 효과를 기대할 수 있다.

본 연구에서는 유전자 집합기반의 마이크로어레이 분석에서 변수간의 연관성을 PCA를 통해 제거할 때 나타나는 분석 결과의 차이에 대해서 살펴보았다.

II. 방 법

1.  마이크로어레이 자료

GeneExpress Omnibus에서 유방암 환자를 대상으로 한 마이크로어레이 자료(GSE3494)를 분석에 사용하였다. 252명의 환자에 대한 정보가 제공되었고, 임파선 전이 여부를 알 수 있는 242명의 환자를 대상으로 분석을 하였다. 실험에 사용된 어레이는 Affymetrix 사의 HU-133A GeneChip이 사용되었다. 실험 시에 발생 가능한 오류를 제거하기 위해 RMA(robust multichip array) 방법과 quantile normalization을 이용하여 정규화(normalization)을 시행하였다.[5]

2.  유전자 집합의 추출

분석에 사용할 유전자 집합은 사용 가능한 생물학적인 경로(biological pathway)를 기반으로 수집하였다. ArrayXpath[6] 데이터베이스를 기반으로  모두 14751개의 프로브들을 467개의 경로에 연결하였다. 이 경로들 중에서 집합 원소의 수가 10 개 이하인 것은 분석에서 제외하여 모두 387개의 유전자 집합을 대상으로 분석을 하였다.

3.  유전자 집합 기반의 분석

한 유전자 집합의 원소들 간의 상호 연관성을 측정하가 위하여 연관성 지표(correlation index, CI)를 각 집합마다 계산하였다. 다음의 식처럼 집합 원소 유전자들의 모든 원소 쌍의 연관 계수(correlation coefficient)의 총합을 원소의 개수로 나누어 연관성 지표를 계산하였다.

CI(correlation index) =

N은 집합의 유전자 개수를 의미하고 CCij는 각 유전자 쌍의 연관 계수를 의미한다.

유전자 집합을 대상으로 한 분석은 global test를 사용하였다. global test는 GLM(generalized linear model)을 이용하여 유전자 집합이 임상적인 표현형과 연관이 있는지 검사하는 방법이다. 이 때 사용하는 구체적인 모델은 다음과 같다.

여기서 α는 절편을 의미하고 β는 회귀 계수를, h는 연결 함수(linker function)를 의미한다. 그리고 y는 임상적인 표현형(phenotype)을 가리킨다. 본 연구에서는 임파선 전이 유무를 임상 표현형으로 정하였다. 이 모델에서 다음의 귀무가설을 검정한다. (p<0.05)

유전자 집합 원소간의 상관관계를 제거하기 위해서 각 유전자 집합에 해당하는 유전자 발현량의 행렬을 PCA를 적용하여 전체 분산에 대한 설명력이 큰 상위 4개의 주성분을 추출하여 global test에 적용하였다.

한 번에 많은 수의 유전자 집합에 대한 검사를 시행하기 때문에 다중 검사 보정(multiple testing correction)이 필요하였다. 이 연구에서는 Benjamini와 Hochberg가 제안한 FDR(false discovery rate)의 개념을 사용하였다. 이 때 q 값은 0.05로 고정하였다.

III. 결 과

1.  유전자 집합의 연관성 지수의 분포

각 유전자 집합의 연관성 지표(CI)를 계산 했을 때 그 분포는 그림1과 같다. 평균은 0.123(분산 0.083)이었고 가장 높은 CI를 보이는 경로는 Lck and Fyn tyrosine kinases in initiation of TCR Activation (CI=0.659)였고 가장 낮은 CI를 보이는 경로는 MAPKinase Signaling(CI=0.024)이었다. 모두 231개의 경로가 평균보다 낮은 CI값을 보였다.

2.  임파선 전이와 유의하게 연관이 있는 경로

q 값이 0.05수준에서 원래 값으로 검정 했을 때는 모두 134개의 경로가 유의하게 나왔고 PCA를 적용한 발현량을 사용했을 때는 137개의 경로가 유의하다고 판정 되었다. 이 두개의 결과에서 공통적으로 나타나는 48개의 경로를 표1에 나열하였다.

Figure 1. Histogram of correlation indexes(CI) of  389 pathways

3.  연관성 지수와 global test p 값과의 관계

전체 387개의 경로에서 원래 값으로는 유의성이 없었으나 PCA로 변환한 값에서는 유의성이 나타난 경로는 78개였다. 이 78개의 경로들에서 연관성 지수가 평균(0.123)보다  높게 나타난 경로는 모두 40개였다.

IV. 결론

본 연구에서는 유전자 집합 기반의 마이크로어레이 분석에서 유전자간의 발현량의 상관관계가 결과에 영향을 미칠 수 있는지를 조사하였다.

Mootha등이 제시한 바와 같이 유전자 집합 기반의 문석은 유전자 단위 분석에서는 발견하기 힘든 현상을 찾아낼 수 있는 기회를 제공한다.[2] 유전자 집합을 정할 때는 생물학적인 지식을 바탕으로 사용하며

Table 1. result of global test
Pathway Names	Correlation Index				Global test p-value	PCA p-value
BioCarta//Hs_Acetylation and Deacetylation of RelA in The Nucleus		0.079	0.008			0.001
BioCarta//Hs_Actions of Nitric Oxide in the Heart		0.051	0.000			0.002
BioCarta//Hs_Cadmium induces DNA synthesis and proliferation in macrophages		0.080	0.000			0.004
BioCarta//Hs_CD40L Signaling Pathway		0.126	0.006			0.002
BioCarta//Hs_Cell Cycle: G1/S Check Point		0.053	0.000			0.001
BioCarta//Hs_Control of Gene Expression by Vitamin D Receptor		0.058	0.002			0.005
BioCarta//Hs_Corticosteroids and cardioprotection		0.108	0.002			0.004
BioCarta//Hs_Cyclins and Cell Cycle Regulation		0.095	0.000			0.001
BioCarta//Hs_Estrogen-responsive protein Efp controls cell cycle and breast tumors growth		0.105	0.002			0.003
BioCarta//Hs_Internal Ribosome entry pathway		0.216	0.009			0.001
BioCarta//Hs_Lissencephaly gene (LIS1) in neuronal migration and development		0.074	0.014			0.001
BioCarta//Hs_Low-density lipoprotein (LDL) pathway during atherogenesis		0.163	0.014			0.001
BioCarta//Hs_Mechanism of Protein Import into the Nucleus		0.202	0.000			0.001
BioCarta//Hs_Melanocyte Development and Pigmentation Pathway		0.152	0.000			0.012
BioCarta//Hs_NFAT and Hypertrophy of the heart (Transcription in the broken heart)		0.030	0.000			0.010
BioCarta//Hs_Nuclear receptors coordinate the activities of chromatin remodeling complexes and coactivators to fa		0.067	0.000			0.003
BioCarta//Hs_p53 Signaling Pathway		0.047	0.001			0.002
BioCarta//Hs_Pertussis toxin-insensitive CCR5 Signaling in Macrophage		0.098	0.000			0.009
BioCarta//Hs_RB Tumor Suppressor/Checkpoint Signaling in response to DNA damage		0.100	0.000			0.006
BioCarta//Hs_Regulation of eIF4e and p70 S6 Kinase		0.029	0.008			0.008
BioCarta//Hs_Rho cell motility signaling pathway		0.046	0.005			0.003
BioCarta//Hs_Role of BRCA1, BRCA2 and ATR in Cancer Susceptibility		0.092	0.009			0.000
BioCarta//Hs_Role of EGF Receptor Transactivation by GPCRs in Cardiac Hypertrophy		0.108	0.002			0.006
Pathway Names	Correlation Index			Global test p-value		PCA p-value
BioCarta//Hs_Role of Mitochondria in Apoptotic Signaling		0.063	0.000			0.000
BioCarta//Hs_Role of PPAR-gamma Coactivators in Obesity and Thermogenesis		0.091	0.004			0.006
BioCarta//Hs_Segmentation Clock		0.087	0.008			0.005
BioCarta//Hs_Signal transduction through IL1R		0.057	0.007			0.004
BioCarta//Hs_Signaling Pathway from G-Protein Families		0.039	0.000			0.016
BioCarta//Hs_Spliceosomal Assembly		0.341	0.001			0.006
BioCarta//Hs_Stathmin and breast cancer resistance to antimicrotubule agents		0.060	0.000			0.006
BioCarta//Hs_T Cell Receptor Signaling Pathway		0.059	0.001			0.004
BioCarta//Hs_Telomeres, Telomerase, Cellular Aging, and Immortality		0.054	0.000			0.001
BioCarta//Hs_The 4-1BB-dependent immune response		0.048	0.015			0.010
BioCarta//Hs_The IGF-1 Receptor and Longevity		0.066	0.003			0.002
BioCarta//Hs_The PRC2 Complex sets Long-term gene Silencing The Modification of Histone Tails		0.081	0.002			0.001
BioCarta//Hs_TNFR2 Signaling Pathway		0.109	0.007			0.001
BioCarta//Hs_TSP-1 Induced Apoptosis in Microvascular Endothelial Cell		0.179	0.000			0.002
GenMAPP//Hs_GPCRs_Class_A_Rhodopsin-like		0.219	0.000			0.001
GenMAPP//Hs_MAPK_Cascade		0.048	0.005			0.011
GenMAPP//Hs_Ovarian_Infertility_Genes		0.076	0.013			0.000
GenMAPP//Hs_S1P_Signaling		0.055	0.003			0.006
KEGG//Hs_Aminoacyl-tRNA_biosynthesis		0.176	0.001			0.001
KEGG//Hs_Aminosugars_metabolism		0.074	0.003			0.006
KEGG//Hs_Arginine_and_proline_metabolism		0.037	0.002			0.015
KEGG//Hs_Fatty_acid_metabolism		0.072	0.000			0.007
KEGG//Hs_One_carbon_pool_by_folate		0.153	0.001			0.002
KEGG//Hs_Riboflavin_metabolism		0.132	0.001			0.001
PharmGKB//Hs_Irinotecan_pathway_liver_cell		0.150	0.001			0.000

이것은 분석 후의 해석을 용이하게 하는 장점도 가지고 있다. 예를 들어, 한 유전자 집합이 임파선 전이와 연관이 있다고 판명됐고 그 집합이 특정 경로를 나타낸다면 그 경로가 임파선 전이에 연관이 있는 것으로 쉽게 해석할 수 있다. 그러나 이러한 분석에서는, 유전자 집합에 속한 원소 유전자들의 연관 관계가 높을 경우 분석 결과에 영향을 미칠 수 있다.

이 연구에서도 연관 관계가 결과에 영향을 미칠 수 있음을 확인하였다. 모두 40개의 평균보다 높은 연관성 지수를 가지며 검사에서 유의하지 않은 경로들이 PCA를 적용했을 때 유의한 결과를 보였다. 회귀 분석에서 변수간의 다중 공선성(multi-collinearity)은 그 회귀식이 의미 없는 것으로 판정될 가능성을 높이는 효과를 지닌다고 알려져 있다. 위의 결과에서 나타난 40개의 경로들은 연관성 지수가 다른 경로들보다 상대적으로 높은 상태에서 PCA를 통해서 연관성을 제거했을 때 global test의 결과가 유의했기 때문에 이런 경로들은 연관 구조를 제거한 후에 검사를 시행하는 것이 적절하다고 판단된다.

표1에서 estrogen responsive finger protein(Efp) 경로(Estrogen-responsive protein Efp controls cell cycle and breast tumors growth)는 특히 흥미 있는 결과였다. Efp는 유방암의 성장에 관여하는 것으로 알려져 있고 임파선 전이와도 연관이 있는 것으로 나타났다.[7] 따라서 Efp 경로는 임파선 전이에 연관된 정보를 제공함과 동시에 임파선 전이 예측이나  진단에 쓸 수 있는 예후 인자(prognostic marker)로서의 가능성을 지닌다.

이 연구에서 저자들은 유전자 집합 기반의 마이크로어레이 자료를 PCA를 통해서 분석하였다. 그리고 이러한 과정이 유전자 집합 단위의 분석에서  집합 원소간의 연관 관계를 제거하여 좀 더 유의한 결과를 도출하는데 도움이 될 수 있음을 확인하였다.

감사의 글

This study was supported by a grant from Korea Health 21 R&D Project, Ministry of Health & Welfare, Republic of Korea (03-PJ1-PG3-21000- 0009, KPGRN-H-10).

참고문헌

1. Tusher et al. Significance analysis of microarrays applied to theionizing radiation response, PNAS Vol 98 April 2001, 5116-5121

2. Mootha et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet 34, 2003, 267-273.

3. J. Goeman, S. A. van de Geer, F. de Kort and J. C. van Houwelingen, 2004, A global test for groups of genes: testing association with a clinical outcome, Bioinformatics 20(1),  2004, 93–9

4. Tian L, Greenberg SA, Kong SW, Altschuler J, Kohane IS, Park PJ. Discovering statistically significant pathways in expression profling studies, Proc Natl Acad Sci USA 102:1, 2005, 3544-9

5. Irizarry, RA, Hobbs, B, Collin, F, Beazer-Barclay, YD, Antonellis, KJ, Scherf, U, Speed, TP  Exploration, Normalization, and Summaries of High Density Oligonucleotide Array Probe Level Data. Biostatistics. Vol. 4, 2003, 2: 249-264

6. Hee-Joon Chung, Chan Hee Park, Mi Ryung Han, Seoko Lee, Jung Hun Ohn, Jihoon Kim, Jihoon Kim,  Ju Han Kim. ArrayXPath II: mapping and visualizing microarray gene expression data with biomedical ontologies and integrated pathway resources using Scalable Vector Graphics. NAR  1;33, 2005, 621-626

7. Suzuki T, Urano T, Tsukui T, Horie-Inoue K, Moriya T, Ishida T, Muramatsu M, Ouchi Y, Sasano H, Inoue S. Estrogen-responsive finger protein as a new potential biomarker for breast cancer. Clin Cancer Res. 1;11(17), 2005, 6148-54.