[scikit-learn] RandomForest Classifier 기계학습 모델로 유전자 발현 분류하기

RandomForest Classifier 기계학습 모델로 유전자 발현 분류하기

---

파이썬(Python)을 이용한 기계학습(Machine Learning, ML)은 다양한 분야에서 사용되고 있는데, 특히 생물정보학(Bioinformatics, BI)에서는 유전자를 이용한 연구에 많이 사용되고 있습니다.

 

 

그래서, 이번에는 BI 분야에서 흔히 사용하는 유전자(여기서는 RNA) 발현 분류를 테마로 잡아 기계학습 실습을 해 보겠습니다. 여기에서 사용할 모델은 랜덤포레스트 분류기(RandomForest Classifier, RFC) 라는 모델인데, 다양한 의사결정나무(DicisionTree, DT)를 여러개 만들어둔 다음 평균 예측을 통해 분류를 하는 지도학습 모델입니다.

 

 

이 분석에 사용한 데이터는 아래 정보를 참고해 주시기 바랍니다.

 

• Transcription profiling by high throughput sequencing of Arabidopsis roots, leaves, flowers and siliques 데이터
• Source: https://www.ebi.ac.uk/gxa/experiments/E-GEOD-38612/Downloads (TPM file)

 

랜덤포레스트 분류기는 파이썬에서 사이킷런(scikit-learn) 이라는 패키지의 RandomForestClassifier 라는 모듈을 이용하면 손쉽게 구현할 수 있습니다. 아래에 파이썬과 패키지 버전을 적어놓으니, 코드를 따라해 보실 분들은 참고해 주세요.

 

• Pyhton 3.10.9
• Pandas 1.5.3
• sklearn 1.2.1
• seaborn 0.12.2

 

 

파일 전처리하기

---

# 패키지 불러오기

import pandas as pd
import seaborn as sns
from sklearn.ensemble import RandomForestClassifier

 

# 파일 불러오기
arabidopsis = pd.read_csv('./Data/Arabidopsis_3/E-GEOD-38612-query-results.tpms.tsv',
                    sep='\t', comment='#', index_col=1)

# leaf, flower, fruit, root만 사용하도록 자르기
arabidopsis = arabidopsis.iloc[:,1:]
arabidopsis.head()

	leaf	flower	fruit	root
Gene Name
NAC001	4.0	10.0	3.0	38.0
ARV1	16.0	16.0	7.0	18.0
NGA3	5.0	7.0	1.0	2.0
DCL1	16.0	21.0	16.0	19.0
MIR838A	6.0	7.0	5.0	8.0

 

 

위 파일은 애기장대(Arabidopsis thaliana) 라는 식물의 RNA 정보입니다. index가 RNA 이름, column name이 각각 잎, 꽃, 종자, 뿌리를 나타내고, 각 값들이 발현량이라고 보시면 됩니다. 하지만, 이 표에는 NaN 이라고 하는 결측치(missing value)가 있는데요, 이 값은 분석 결과에 영향을 끼치므로 이번에는 삭제해 주도록 하겠습니다.

 

 

# NaN 값 개수 확인
arabidopsis.isnull().sum()

leaf      5750
flower    2121
fruit     3375
root      4616
dtype: int64

 

# 계산 편의상 missing value를 0으로 치환한다
arabidopsis = arabidopsis.fillna(0)

 

 

사이킷런을 이용해서 기계학습을 하면 대체로 넘파이(Numpy) array 자료형으로 입력을 하게 됩니다. 이 때, 입력할 데이터를 X로 두고, 라벨을 y로 두는 것이 일종의 불문율이랍니다.

 

 

X = arabidopsis.values.T  # RNA를 가로, 기관 이름을 세로 방향으로 배열
y = [1, 2, 3, 4]  # 숫자로 기관 이름을 분류

 

 

모델 학습시키기

---

사이킷런에서 제공하는 대부분의 기계학습 모델은 클래스로 구현되어 있습니다. 그래서 하이퍼파라미터(Hyperparameter)를 적절하게 조절한 뒤 fit() 메소드를 이용해 학습을 시키면 됩니다. 이 때, 아래 정보를 이용해서 하이퍼파라미터를 조절하면 됩니다. 자세한 내용은 사이킷런 홈페이지에 설명되어 있으니 그 쪽을 참고하셔도 좋습니다.

 

• n_estimator: 모델이 전개할 의사결정나무의 개수
• random_state: 랜덤 시드 (똑같은 시드로 코드를 실행하면 분석을 재현할 수 있습니다)
• max_depth: 의사결정나무의 최대 깊이

 

rfc = RandomForestClassifier(n_estimators=2000, random_state=12345)
rfc.fit(X, y)

RandomForestClassifier(n_estimators=2000, random_state=12345)

 

 

 

Feature Selection

---

Feature Selection은 기계학습 모델의 학습 결과를 토대로 각 feature(여기서는 RNA)의 중요도(importance)를 평가한 뒤, Feature Impotance가 높은 대상만 선별하는 과정입니다.

 

 

우리도 공부를 할 때 관련 지식을 갖추면 쉬워지고, 쓸모없는 지식을 같이 배우면 어려워지는것처럼 기계학습 모델 역시 이러한 특성을 가지고 있습니다. 여기서는 각 부위를 분류할 수 있는 RNA는 어떤게 있을지를 알아보면서 코드를 구현해 보도록 하겠습니다.

 

 

m = rfc.feature_importances_ > 0  # Feature Importance > 0 만을 선택하기
fi = arabidopsis.index[m]  # 주요 feature를 마스킹으로 선별하기
arabidopsis_fi = arabidopsis.loc[fi]

# 부위별 RNA 발현 수준을 clustermap으로 표현
sns.clustermap(arabidopsis_fi, standard_scale=0)

<seaborn.matrix.ClusterGrid at 0x1d7f09347c0>

 

 

위의 클러스터맵을 확인해보니, 음영이 뚜렷하게 구분되는 부분이 몇군데 보입니다. 여기서 밝은 음영 부분을 공유하고 있는 RNA 집단은 해당 기관을 분류할 때 주요한 영향을 끼치는(Feature Importance가 높은) RNA라고 볼 수 있겠습니다.

 

 

규모가 큰 Sample을 이용한 발현 수준 분류

---

이번에는 애기장대 17개에서 RNA 정보를 표 형태로 정리해 놓은 파일을 이용해 보겠습니다. 아래 파일을 다운로드한 뒤, 위에서 진행한 분석과 같이 실행해 보겠습니다. 데이터 정보는 아래와 같습니다.

 

 

• Transcription profiling by high throughput sequencing from seedling, root, and floral buds of 17 Arabidopsis thaliana accessions 데이터
• Source: https://www.ebi.ac.uk/gxa/experiments/E-GEOD-53197/Downloads (TPM file)

 

 

ara_17 = pd.read_csv('./Data/Arabidopsis_17/E-GEOD-53197-query-results.tpms.tsv',
                     sep='\t', comment='#', index_col=1)

ara_17 = ara_17.iloc[:, 1:]
# 계산 편의상 missing value를 0으로 치환한다
ara_17 = ara_17.fillna(0)

ara_17.head()

	Bur-0, aerial part	Bur-0, floral bud	Bur-0, root	Ct-1, aerial part	Ct-1, floral bud	Ct-1, root	Edi-0, aerial part	Edi-0, floral bud	...	Ws-0, floral bud	Ws-0, root	Wu-0, aerial part	Wu-0, floral bud	Wu-0, root	Zu-0, aerial part	Zu-0, floral bud	Zu-0, root
Gene Name
NAC001	5.0	13.0	43.0	4.0	4.0	40.0	4.0	12.0	...	14.0	40.0	4.0	10.0	38.0	3.0	7.0	51.0
ARV1	7.0	9.0	10.0	4.0	7.0	5.0	5.0	7.0	...	7.0	9.0	5.0	5.0	8.0	6.0	6.0	8.0
NGA3	4.0	6.0	1.0	4.0	5.0	0.2	5.0	5.0	...	5.0	1.0	6.0	5.0	0.9	6.0	6.0	1.0
DCL1	11.0	24.0	15.0	17.0	26.0	17.0	14.0	30.0	...	15.0	11.0	11.0	15.0	10.0	10.0	22.0	14.0
MIR838A	2.0	8.0	2.0	0.7	7.0	3.0	4.0	6.0	...	5.0	2.0	1.0	3.0	0.7	2.0	16.0	3.0

5 rows × 51 columns

 

 

X = ara_17.values.T

 

 

이번 데이터는 column 개수가 많아져서 y 값을 설정할 때 분류를 다르게 할 필요가 있습니다. 다만, 각 column이 접두사는 다르지만, 각각 aerial part, floral bud, root의 셋 중 하나인 것으로 확인되어 이를 기준으로 분류하도록 하겠습니다.

 

 

각 column의 이름은, Bur-0,aerial part 와 같이, 부위의 이름과 접두사의 사이에 콤마와 띄어쓰기(, )로 구분자가 있습니다. 이것을 기준으로 문자를 분리한다면, 각각의 부위별 이름으로 column을 통합할 수 있을겁니다.

 

 

part = [x.split(', ')[1] for x in ara_17.columns]
part = list(set(part))  # 중복되는 요소를 합친 후 리스트화

part

['aerial part', 'floral bud', 'root']

 

y = [part.index(x.split(', ')[1]) for x in ara_17.columns]

 

 

이렇게, floral bud 관련 부위는 모두 0으로, aerial part 관련 부위는 모두 1로, root 관련 부위는 모두 2로 지정하도록 y를 만들어 놓았습니다. 그러면, 각 부위를 구분할 수 있도록 RFC를 이용해 기계학습을 하고, Feature Importance를 확인해 보도록 하겠습니다.

 

 

rfc = RandomForestClassifier(n_estimators=2000, random_state=0)
rfc.fit(X, y)

 

RandomForestClassifier(n_estimators=2000, random_state=0)

 

 

m = rfc.feature_importances_ > 0
ara_17_fi = ara_17.index[m]

sns.clustermap(ara_17.loc[ara_17_fi], standard_scale=0)

<seaborn.matrix.ClusterGrid at 0x1d7901736a0>

 

 

분석 결과를 살펴보면, 애기장대의 종류(x축 이름)에 관계없이, 전체적인 RNA 발현(Global Pattern)은 비슷한 경향을 보이고 있습니다. 따라서, 이 결과를 근거로, 식물 실험을 할 때는 개체별로 조직에 대한 큰 차이는 없다고 가정하고 실험을 진행할 수 있다는 사실을 알 수 있습니다.