안녕하세요 
오랜만에 스티밋에 글을 씁니다. 저는 제 소개에서 언급했듯 암 시스템 생물학 분야에서 생명정보학적 방법으로 암 바이오마커 및 치료제 관련 연구와 함께 암을 통해 생명 현상에 대한 이해를 도모하고 있습니다. 제 분야가 아직은 생물학에 있어 new comer다보니 많은 분들이 프로그래머, 통계학자, 혹은 데이터 분석가 등으로 인식하지만, 컴퓨터를 실험 도구로 사용하는 암 생물학자 입니다.
오늘은 Python으로 데이터분석을 수행할 때, 가장 널리 사용하는 Pandas 모듈로 행렬 유사 데이터 구조를 다루는 법을 살펴 보려 합니다. 데이터는 대중에 공개된 The Cancer Genome Atlas (TCGA) 데이터 중 대장암 환자들의 임상데이터를 이용하려 합니다. Pandas는 통계학자겸 사업가인 Wes McKinney가 데이터 분석을 위한 소프트웨어 라이브러리로 개발 한 것으로 시계열 데이터나 cross-sectional data (행렬 유사 형태 데이터)를 다루는데 큰 도움을 줍니다.
참고로 금융 데이터 분석과 시스템 트레이딩에도 Pandas를 활용하곤 합니다. 참고:Python For Finance: Algorithmic Trading
오늘 다룰 내용은 데이터 분석에 앞서 어떤 특징이 있는지 살펴보는 단계로 저희는 보통 데이터 탐색 단계라 합니다. 보통 데이터 분석을 수행하는 경우, 위 순서를 통해 데이터를 살펴본 후 적정한 분석 기법을 리스트업하고 분석 업무노가다를 시작합니다.
이 포스팅을 작성하는데 있어 저는 Python 3에 Jupyter notebook으로 작성했고, 사용한 라이브러리는
입니다. 오늘 포스트에서는 Pandas 모듈만 이용합니다.
생물학적 기초에 대해 언급하기엔 지면이 협소해 이 포스팅에서는 생물학적 배경을 다루기는 힘듭니다. 이 분야에 관심이 있지만 생물학 분야에 대해 어려움을 겪는다면 저는 보통 교양서로는 후쿠오카 신이치 4부작를 권하는 편이고, 기초 교재로는 방통대 교재로도 선정되었던 알기 쉽고 재미있는 분자생물학을 권하는 편입니다.
TCGA 데이터세트 다운로드는 cBio Portal이나 NCI GDC에서 다운로드 받을 수 있지만, 여기서는 UCSC Cancer Browser에서 다운로드 받도록 합니다. (신규 브라우저인 Xena는 좀 불편...)
UCSC Cancer Browser 접속 후 아래 그림과 같이 Add Datasets 버튼을 클릭
팝업된 Datasets 창에 Search 박스에 "colon"을 입력해 대장암 관련 데이터 세트만 필터링 합니다. 그 후, TCGA의 대장암 데이터 세트인 "TCGA colon adenocarcinoma"를 엽니다. 추후 messenger RNA (mRNA) 발현량 데이터 분석까지 계속 이용할 생각이기 때문에, 여러 데이터 세트 중, "COAD gene expression (IlluminaHiSeq)"을 다운로드 합니다. 이 파일의 의미는 'colon adenocarcinoma (COAD) 암 환자들의 mRNA 발현량을 Illumina 사의 HiSeq 이란 장비로 측정한 것' 입니다.
압축을 풀면 아래와 같이 파일이 있습니다. 환자들의 임상 정보 (성별, 연령, 진단시 키 및 몸무게 등 기초 정보와 생존 기간과 생존 여부, TNM 스테이지, 분자진단학적 분류 등)는 clinical_data 라는 파일에 들어있습니다. 이 파일을 python의 pandas 모듈로 불러드릴 겁니다.
Jupyter notebook에서 그래프를 바로 표시토록 하는 설정과 함께 라이브러리 import, 그리고 작업 폴더 위치를 상수로 미리 선언해놓습니다. 작업할 폴더 위치에 문자열 formatting slot ({0})을 포함 하는걸 개인적으로 선호합니다.
%matplotlib inline
import pandas as pd
import numpy as np
import seaborn as sns
_workingdir_slot = "./TCGA_COAD_exp_HiSeqV2-2015-02-24/{0}"
clinical_data를 read_csv 함수로 읽어들입니다. 컬럼간 구분자를 Tab (\t)으로 설정하고, key로 사용할 index 컬럼은 첫 번째 (0번) 컬럼으로 합니다. encoding은 기본으로 utf-8을 사용하지만, TCGA의 clinical_data의 경우 기본 인코딩을 Western으로 해놔서 기본으로 읽어들이면, 오류가 나곤 합니다.
df_clinic = pd.read_csv(_workingdir_slot.format("clinical_data"), sep="\t", index_col=0, encoding="ISO 8859-1")
df_clinic.head()
읽어들인 파일의 내용을 살펴보기 위해 head() (Unix 및 Linux의 head와 동일합니다.)로 출력합니다.
sampleID 컬럼이 굵게 표시되는 것은, CSV 파일을 읽어들일 때, index로 설정했기 때문입니다. Row 선택 시 index (또는 key가 됩니다.)
위 이미지 맨 하단에 표시된 것처럼, 이 데이터는 155개의 column을 가지고 있습니다. column 목록만 출력해보겠습니다. (TCGA의 clinical_data의 column 관련 설명은 웹에서 확인 할 수 있습니다.)
df_clinic.columns.tolist()
column 중에서, 성별 (gender), 신장 (height), 체중 (weight; 암 확진시 체중), 연령 (age_at_initial_pathologic_diagnosis; 암 확진시 연령), 샘플 종류 (종양 근처 peri tumor, tumor 등), 샘플 종류 id, 생존기간 (_OS; overall survival), 생존여부 (_OS_IND; 0-생존, 1-사망), microsatellite instability (MSI) 분류 (CDE_ID_3226963) 만 선택하기 위해 컬럼명을 list 타입으로 생성 후, 데이터프레임에 list 데이터를 column key로 줍니다.
target_columns = [
"gender", "height", "weight",
"age_at_initial_pathologic_diagnosis",
"sample_type", "sample_type_id",
"_OS", "_OS_IND", "CDE_ID_3226963"]
df_clinic[target_columns]
원하는 column들만 선택하면 출력은 아래와 같습니다. 캡쳐에서 누락됐지만, 총 550 명의 환자 데이터가 있습니다.
이제 데이터내 포함된 샘플들 특성을 좀 살펴봐야 겠습니다. 위 이미지를 보면 gender, sample_type, CDE_ID_3226963 columns는 통계학에선 categorical variable라 부르는 데이터 타입입니다. 또 cheat sheet을 좀 빌리자면, sample_type_id와 _OS_IND도 categorical variable 입니다. TCGA-5M-AAT4-01부터 4개 샘플들은 값이 누락 (NaN)되어 있습니다. 샘플 획득이 어려워, 데이터 개수가 매우 중요할 경우 보간법 (통계학의 회귀식이라든지, 기타 machine learning 기법 등)을 이용해 값을 예측 값으로 채워넣을 수도 있지만, 적정한 데이터 수가 확보된다면, 불완전한 데이터는 버리는게 좋습니다.
이제 성별로 각 몇 샘플이 있는지를 보겠습니다.
df_clinic[target_columns].groupby("gender").size()
전체 샘플 수가 550이지만, 위 결과가 539명인 걸로 봐서 위에서처럼 성별 기입이 누락되었거나 잘못 기입된 환자 수가 11명 있나 봅니다.
df_clinic[target_columns][(df_clinic["gender"]!="FEMALE") & (df_clinic["gender"]!="MALE")]
잘못 기입된 샘플은 없고 모두 누락되어 있습니다. 이제, 데이터 세트 내 종양 샘플 (primary tumor) 개수를 살펴 봅니다.
df_clinic[target_columns][(df_clinic["sample_type_id"]==1)]
총 461 샘플이 있습니다. 이제, 종양 샘플만 가져와서 새로운 데이터프레임을 만든 후, 성별, 신장, 체중 및 그 외 모든 값이 누락되지 않은 샘플만 선택합니다.
df_clinic_tumor = df_clinic[target_columns][(df_clinic["sample_type_id"]==1)]
df_clinic_tumor.dropna(inplace=True)
df_clinic_tumor
224 샘플이 남았습니다. 종양 샘플만 추출한 데이터 내 성비가 어떻게 되는지 봅니다.
df_clinic_tumor.groupby("gender").size().plot(kind="pie", autopct="%d", fontsize=14, figsize=(6, 6))
Pie chart는 정보 전달력이 그다지 좋지도 않을 뿐더러 왜곡 시킬 수 있어서 그리 권하는 양식이 아닙니다. 체중, 신장, 생존 기간, 최초 암진단시 연령 등은 비범주형이므로 데이터가 어떻게 분포되어있는지 산포도가 있으면 파악하는데 도움이 되니 하나 그려봅니다.
pd.plotting.scatter_matrix(df_clinic_tumor[["height", "weight", "age_at_initial_pathologic_diagnosis", "_OS"]], alpha=0.2, figsize=(10, 10), diagonal='hist')
scatter plot matrix라 부르는 형태로 좌측 상단에서 우측 하단 가운대 대각선을 가로지르는 분포는 각 변수가 어떻게 분포되어 있는지 보여줍니다. 신장 (height)의 분포를 보면 대부분이 150cm ~ 190cm에서 종모양 (정규분포형태)를 유지하고 있습니다만, 80cm 근처에 데이터가 하나 있습니다. 이건 이상치 값 (outlier)인지 기입 실수인지 데이터 분석 수행 전에 잊지말고 살펴봐야 합니다.
분포도를 기준으로 가로방향이나 세로방향의 산포도는 동일합니다. (파란선은 손으로 대충 그려놓은 것입니다만, Pearson correlation coefficient 를 계산해보면 대충 rho 값이 0.5 이상일 겁니다.) 다만 X, Y 축이 바뀌어서 전치 (transpose)되었기 때문에 대칭으로 나타납니다.
scatter matrix는 개별 변수의 분포와 함께 변수간 상관관계를 한눈에 볼 수 있어, 데이터 탐색 단계에서 꼭 찍어보는게 좋습니다.
분량 관계로 이번 포스팅은 여기서 멈추고 다음번에는 column 값을 대상으로 계산을 하는 것과 시각화 방법 등에 대해서 다뤄보겠습니다.
덧글: 데이터 분석과 관련한 글을 틈나는데로 준비하다보니 알게된 것이 있습니다. 스티밋 마크다운 지원 참 좋긴한데, 코딩 syntax highligh나 레이텍 (LaTex) 미지원으로 코드나 수식 기입이 불편하거나 불가능하고, 또 몇 몇 markdown extension에서 지원하는 [TOC] keyword나 inner link를 활용한 참조 문헌 등록 같은 것들이 없다는 점은 앞으로 개선되면 좋겠습니다.