SAS - 산점도와 히스토그램

descriptive analysis에서 변수별로 기본적으로 해보는 것이 바로 산점도와 히스토그램이라고 할 수 있습니다. 산점도는 두 변수의 관계를 파악할 때 쓰이며, 히스토그램은 한 변수의 분포를 알아보기 위해 쓰입니다. 이번 포스팅에서는 SAS로 산점도와 히스토그램을 그리는 법을 정리해보았습니다.

기본 산점도

proc gplot data=data;
    plot var1*var2;
run;

• 기본 산점도는 위와 같이 간단하게 구할 수 있습니다. 근데 분석을 하다보면, 산점도를 그려도 두 변수의 관계가 확실히 보이지 않을 때가 있습니다. 이 경우에는 선형 회귀 선을 함께 그림으로써 그 경향성을 더욱 잘 파악할 수 있는데요. 링크 에는 산점도와 함께 선형회귀선을 그리는 코드를 보여주고 있습니다.
  

산점도 with 선형회귀선

ods graphics off;
proc reg data=data;
   model y= x;
   ods output ParameterEstimates=PE;
run;

data _null_;
   set PE;
   if _n_ = 1 then call symput('Int', put(estimate, BEST6.));   
   else            call symput('Slope', put(estimate, BEST6.)); 
run;

title 'y*x';
proc sgplot data=data noautolegend;
   reg y=y x=x;
   inset "Intercept = &Int" "Slope = &Slope" /
         border title="Parameter Estimates" position=topleft;
run;

• 코드를 정리하면 위와 같은데요. 선형회귀분석의 결과로 나온 파라미터들을 임시 데이터셋에 저장해놓고 이를 매크로변수에 저장하여 산점도를 그릴 때 불러오는 방식입니다. 만약 추정된 파라미터 값이 굳이 필요없으면 마지막 proc sgplot 만 실행하면 되고, reg y=y x=x; 라인만 추가하면 선형회귀선이 산점도에 나타나게 됩니다.
  

히스토그램 및 밀도함수(density plot)

• 히스토그램과 밀도함수는 아래의 코드로 그릴 수 있습니다. histogram statement와 함께 density statement를 추가하며 적당한 밀도함수를 그려주게 되며, 정규분포 밀도함수를 함께 그려주어 해당 데이터가 정규분포와 얼마나 유사한지를 보여줄 수 있습니다.
  

proc sgplot data=data;
  histogram var_name ;
  density var_name / type=normal legendlabel='Normal' lineattrs=(pattern=solid);
  density var_name   / type=kernel legendlabel='Kernel' lineattrs=(pattern=solid);
  keylegend / location=inside position=topright across=1;
  xaxis display=(nolabel);
 run;

SAS Proc report 옵션들

/* 2017.9.20 by. 3months */

SAS Base programming 자격증 공부하다가 proc report 관련 정리를 하였습니다. proc report는 데이터셋을 통해 무언가를 report할 때 쓰이는 SAS의 프로시져입니다.

데이터 정의

data a;
 input a$ b;
 cards;
 a 1
 a 2
 a 3
 a 4
 c 5
 c 6
 b 7
 b 8
 b 9
 b 10
 ;
run;

proc report 프로시져

proc report data=a;
   column a; * 이 컬럼에 대한 리포트를 하겠다;
   define a/group; * 중복값을 제외하고 group으로 묶겠다 ;
run;

proc report data=a;
 column a;
 define a/across; * across 옵션은 갯수를 센다;
run;

proc report data=a;
 column b;
 define b/analysis; * b 컬럼 분석 sum을 낸다. ;
run;

proc report data=a;
 column b;
 define b/display; * 그냥 그대로 보여준다.;
run;

차례 대로 각각의 proc report에 대한 결과

출처

http://blog.naver.com/sonobono88/220271045853

강화학습 Q-learning on non-deterministic worlds

출처

• https://github.com/fabulousjeong/OPENML-Deep-Reinforcement-Learning/blob/master/ipython/Lab5_Q-learning%20in%20non-deterministic%20world.ipynb
• http://hunkim.github.io/ml/
  

Frozen Lake 게임이 non-deterministic인 경우 Q를 어떻게 학습하는지 알아보겠습니다. 이 포스팅은 Frozen Lake 게임과 강화학습의 기본 개념을 안다는 전제로 김성훈 교수님의 강화학습 강의를 제 나름대로 정리한 포스팅입니다. 우선 non-deterministic이 무슨 의미인지부터 살펴보면 non-deterministic 이라는 것은 원하는대로 되지 않는다는 의미입니다. (stochastic이라고도 합니다.) 이런 경우에는 무조건 Q가 가라는 대로 가는 것보다는 "기존의 믿음" 을 차근 차근 업데이트 해나가는 것이 Q 학습이 더 잘 된다고 합니다. 이는 현실 세계와도 많이 닮아있습니다. 이를 멘토에 비교해볼 수 있는데, 무조건 멘토가 하라는 대로 하기보다는 결국 무언가를 하기 위해서는 타인의 조언 보다는 자기 자신의 신념이 더 중요하다고 할 수 있습니다.

이를 식으로 표현하면 기존 알고리즘은  

Q(state, action) = R + λmax(Q(new state)) (λ < 1)

이었는데

Q(state, action) = (1-a)Q(state, action) + a(R+λmax(Q(new state)) (a = learning rate)

을 사용하게 됩니다. learning rate 개념을 통해 Q를 한 번에 업데이트 하지 않고, 기존에 믿음을 토대로 조금씩 업데이트 해나가는 것입니다. 이를 구현하는 코드를 보겠습니다. 기존 코드에서 Q 업데이트 하는 부분만 달라지기 때문에 크게 다르지는 않습니다.

라이브러리 로드

import gym
import numpy as np
import matplotlib.pyplot as plt 
from gym.envs.registration import register
import random as pr

environment 생성 및 하이퍼 파라미터 초기화

env = gym.make('FrozenLake-v0')
# Q-table 초기화
Q = np.zeros([env.observation_space.n, env.action_space.n])

# Learning Rate
learning_rate = 0.85

# discount factor
dis = .99
    
# learning 횟수
num_episodes = 2000

# 학습 시 reward 저장
rList = []

Q-learning 알고리즘

for i in range(num_episodes):
    
    # env 리셋
    state = env.reset()
    rAll = 0
    done = False
    
    #decaying E-greedy
    e = 1./((i//100)+1)
    
    # Q-table Learning algorithm
    while not done:
        
        # Choose an action by e greedy
        '''
        if np.random.rand(1) < e:
            action = env.action_space.sample()
        else:
            action = np.argmax(Q[state, :])
        '''
        
        # add Random noise
        action = np.argmax(Q[state, :]+ np.random.randn(1, env.action_space.n)/(i + 1))
        
        # new_state, reward 업데이트 
        new_state, reward, done, _ = env.step(action)
        
        # update Q-table
        Q[state, action] = (1-learning_rate) * Q[state, action] + learning_rate*(reward + dis*np.max(Q[new_state, :]))
        
        rAll += reward
        state = new_state

    rList.append(rAll)

action을 선택하는 부분에서 e-greedy 방법을 사용하는 방법과 random noise를 사용하는 방법이 있는데 둘 다 exploit & exploration을 구현한 것으로 둘 중 하나를 사용하면 됩니다. 밑에 #update Q-table 부분에 Q를 업데이트하는 코드를 볼 수 있습니다. learning rate를 통해 기존의 믿음을 얼마나 받아들일지, Q의 판단을 얼마나 받아들일지를 조정합니다.

print('success rate: ', str(sum(rList)/num_episodes))
print('Final Q-table values')
print(Q)
plt.bar(range(len(rList)), rList, color = 'blue')
plt.show()

('success rate: ', '0.4135')
Final Q-table values
[[  6.05216575e-01   4.77802714e-03   7.66890982e-03   1.14402756e-02]
 [  1.89070559e-04   2.71412272e-05   2.94641678e-03   6.06539650e-01]
 [  6.20782244e-03   3.96037865e-03   2.02160738e-03   5.88158391e-01]
 [  1.38851811e-04   1.81853506e-04   1.08700028e-03   3.36662917e-01]
 [  5.63057417e-01   1.35053120e-03   4.64530097e-04   1.81633578e-03]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  8.64124213e-05   2.51760080e-04   7.21585211e-02   3.41902605e-05]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  3.70706358e-05   1.37044981e-03   2.40643636e-04   7.95889822e-01]
 [  1.09162680e-04   5.60870413e-01   9.16914837e-04   2.64997619e-04]
 [  9.20012820e-01   5.54630225e-05   7.32392792e-05   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  4.37606936e-04   4.03103644e-04   9.16837079e-01   3.65822303e-04]
 [  0.00000000e+00   0.00000000e+00   2.34205607e-03   9.98477546e-01]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]]

success rate는 0.4135인데, 이는 0.4135의 확률로 Goal에 도착한다는 것입니다. non-deterministic world 에서는 상당히 높은 성공률입니다. 또 이 부분에서 중요한 점은 이 non-deterministic world에서의 강화학습 알고리즘의 성공률은 사람의 성공률 보다 높을 수 있다는 것입니다. 왜냐하면 deterministic world에서는 사람이 한 눈에 게임의 rule을 알 수 있고, 매우 쉽게 Goal을 찾을 수 있는 반면, non-deterministic world에서는 사람에게는 Goal을 찾아가는 것이 쉽지 않습니다. 하지만 non-deterministic world에서도 확률적으로 패턴이 존재하고, 이 사람이 이해하기 힘든 패턴속에서 강화학습 알고리즘은 에러를 최소화할 최적의 경로를 찾는다는 것입니다. 이 부분이 진정한 강화학습의 강점이 아닐까 생각해보았습니다.

강화학습 Q-Learning 기초 실습

Q-Learning을 통해 Frozen Lake 문제를 해결하는 방법에 대한 포스팅입니다. 실습 코드의 내용 및 출처는 김성훈 교수님의 강화학습강의( https://www.youtube.com/watch?v=VYOq-He90bE&feature=youtu.be) 입니다.

이 코드를 이해하기 위해서는 강화학습에 대한 기본적인 내용과 Frozen Lake 문제에 대하여 알아야합니다. Frozen Lake 문제는 간단히 말하면 위와 같은 상황에서 state를 이동하여 H(Hall)에 빠지지 않고 G(Goal)로 이동하는 문제입니다. 이 포스팅에서 강화학습에 대한 기본 내용들은 다루지 않겠고, 코드만 정리하여 포스팅 하려고 합니다.

환경 세팅

실습에서 사용하는 환경은 gym 패키지 환경인데, 이 환경에서 강화학습으로 해결할만한 여러 toy example들을 로드할 수 있습니다.

gym이 설치되지 않은 경우 pip install gym 으로 설치하시면 됩니다.

import gym
import numpy as np
import matplotlib.pyplot as plt
from gym.envs.registration import register

register(
    id='FrozenLake-v3',
    entry_point = 'gym.envs.toy_text:FrozenLakeEnv',
    kwargs={'map_name':'4x4',
           'is_slippery':False}
)

env = gym.make('FrozenLake-v3')

위와 같이 FrozenLake 게임을 로드할 수 있습니다.

Q-Table 초기화 및 하이퍼 파라미터 초기화

# Q Table을 모두 0으로 초기화 한다. : 2차원 (number of state, action space) = (16,4)
Q = np.zeros([env.observation_space.n, env.action_space.n])

# discount 정의 => 미래의 reward를 현재의 reward 보다 조금 낮게 본다.
dis = 0.99

# 몇 번 시도를 할 것인가 (에피소드)
num_episodes = 2000

# 에피소드마다 총 리워드의 합을 저장하는 리스트
rList = []

Q-learning은 결국 Q 를 배우는 것입니다. Q는 주어진 state에서 어떤 action을 취할 것인가에 대한 길잡이입니다. 이 문제에서 Q라는 2차원 배열에는 현재 state에서 action을 취할 때 얻을 수 있는 reward를 저장하고 있습니다. 이 Q 2차원 배열에서 argmax 함수를 이용하면 어떤 action을 취할지를 얻어낼 수 있는 것이죠. 이 문제에서 state는 16, action은 4입니다. (4x4 Frozen Lake 그리고 action은 위, 아래, 왼쪽, 오른쪽 4개)

Q-learning 알고리즘에서 Q를 업데이트하는 것

Q(state, action) = R + max(Q(new state)) 로 업데이트하게 됩니다. R은 reward로 게임 내부에서 지정되는 값입니다. 현재 state에서의 어떤 action을 취할 때의 Q 값은 그 action을 통해 얻어지는 reward와 action으로 변화된 state에서 얻을 수 있는 reward의 최댓값을 더한 것입니다. 즉 의미는 현재 리워드와 미래에 가능한 리워드의 최대치를 더하는 것입니다.

근데 이 때, Q(state, action) = R + discount * max(Q(new state)) 로 미래 가능한 리워드에 1 미만의 discount factor를 곱해주어(이 실습에서는 0.99) 미래 리워드에 약간의 패널티를 주기도 하는데 이런 방식을 통해 Q가 조금 더 optimal한 방법으로 learning 될 수 있습니다. 이 때 discount는 하이퍼 파라미터(hyperparameter)로 여러번 시도하면서 좋은 값을 찾을 수 있습니다.

Q-learning 알고리즘에서 Action을 고르는 것

강화학습에서 action 을 고르는 것은 이슈중 하나입니다. 이 실습에서는 두 가지 방법을 쓰는데 첫 번째, random noise 방식, 두 번째, E-Greedy 방식입니다. 둘 다 exploit & exporation 방법을 구현한 것으로 볼 수 있는데 약간의 차이가 있습니다. 둘의 공통적인 아이디어는 '무조건 Q 가 시키는대로만 가지말고 새로운 길로도 가보자는 것' 입니다. 왜냐하면 Q가 시키는 것이 optimal한 것이 아닐 수 있기 때문입니다.

1. random noise 방식

random noise 방식이란 현재 state에서 가능한 action에 따른 Q값(총 4가지)에 random noise를 주어서, 이 것이 최대값이 되는 action을 action으로 선택하게 됩니다. 그러면 무조건 최대 Q 값만 따르지 않고 가끔은 다른 action을 취하기도 합니다. random noise를 (i+1)로 나누는 것은 여러번 해보면서 어느 정도 Q 값이 optimal하게 될 것이기 때문에 exploration을 줄이고 exploit 위주로 한다는 것입니다.

for i in range(num_episodes) : 
    state = env.reset()
    rAll = 0
    done = False
    
    # Q learning 알고리즘
    while not done : 
        # Action 중에 가장 R(Reward)이 큰 Action을 고른다. 
        # 이 때, random noise 방식으로 decaying Exploit & Exploration 구현 
        action = np.argmax(Q[state, :] + np.random.randn(1, env.action_space.n) / (i+1))
        
        # 해당 Action을 했을 때 environment가 변하고, 새로운 state, reward, done 여부를 반환 받음
        new_state, reward, done, _ = env.step(action)
        
        # Q = R + Q 
        Q[state, action] = reward + dis * np.max(Q[new_state, :])
        
        rAll += reward
        state = new_state
        
    rList.append(rAll)

2. E-greedy 방식

E-Greedy 방식은 어떠한 확률값 e를 주어, e의 확률로 exploration한다는 것입니다. 예를 들어 e=0.99 이면 99%의 확률로 exploration 하고, 1%의 확률로 exploit해서 새로운 길을 찾게 됩니다. 이 값을 (i/100) 으로 나눈 것은 위 1번과 마찬가지로 여러버 해보면서 Q 값이 optimal 해질 것이기 때문에 exploration을 줄이면서 exploit 위주로 한다는 뜻입니다.

for i in range(num_episodes) : 
    state = env.reset()
    rAll = 0
    done = False
    
    # exploration의 확률 (decaying)
    e = 1./((i / 100) + 1)
    
    # Q learning 알고리즘
    while not done : 
        
        # E-Greedy 알고리즘으로 action 고르기
        if np.random.rand(1) < e :
            action = env.action_space.sample()
        else : 
            action = np.argmax(Q[state, :])
        
        # 해당 Action을 했을 때 environment가 변하고, 새로운 state, reward, done 여부를 반환 받음
        new_state, reward, done, _ = env.step(action)
        
        # Q = R + Q 
        Q[state, action] = reward + dis * np.max(Q[new_state, :])
        
        rAll += reward
        state = new_state
        
    rList.append(rAll)

최종 결과

Success rate는 Goal까지 실패 안하고 갈 확률입니다. 그래프를 보면 초반에만 실패하고 나중에는 다 성공하도록 Q 가 학습되어 나간다는 것을 볼 수 있습니다.

print("Success rate : "+str(sum(rList) / num_episodes))
print("Final Q-Table Values")
print(Q)

plt.bar(range(len(rList)), rList, color="blue")
plt.show()

Success rate : 0.807
Final Q-Table Values
[[ 0.94148015  0.95099005  0.95099005  0.94148015]
 [ 0.94148015  0.          0.96059601  0.95099005]
 [ 0.95099005  0.970299    0.95099005  0.96059601]
 [ 0.96059601  0.          0.95099005  0.        ]
 [ 0.95099005  0.96059601  0.          0.94148015]
 [ 0.          0.          0.          0.        ]
 [ 0.          0.9801      0.          0.96059601]
 [ 0.          0.          0.          0.        ]
 [ 0.96059601  0.          0.970299    0.95099005]
 [ 0.96059601  0.9801      0.9801      0.        ]
 [ 0.970299    0.99        0.          0.970299  ]
 [ 0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.        ]
 [ 0.          0.9801      0.99        0.970299  ]
 [ 0.9801      0.99        1.          0.9801    ]
 [ 0.          0.          0.          0.        ]]

Introduction to Next Generation Sequencing data analysis

아래 글을 참고하였습니다.

http://genestack-user-tutorials.readthedocs.io/guide/intro-to-ngs.html

NGS 기술은 차세대 시퀀싱 기술이라고 말하며, DNA의 시퀀스 정보를 빠르고 값 싸게 읽어내는 기술이다. NGS 기술이란 광범위한 의미인데, 예를 들어, WGS, RNA-Seq, WES, WGBS, ChIP-Seq 등을 NGS 기술이라고 말한다. WES(Whole Exon Sequencing) 데이터의 예를 들어서 NGS 데이터 분석하는 과정에 대해서 알아보자. WES이란 DNA 상의 Exon 부위의 염기서열을 읽어내는 기술을 말한다.

우선 큰 그림을 살펴보면, 일반적인 NGS 데이터 분석은 raw reads quality control, preprocessing, mapping, post-alignment processing, variant calling, annotation, and prioritization 으로 이루어져 있다.

이제 차례대로 WES 데이터를 처리할 때 어떤 스텝을 거쳐야 하는지 알아보자. 시퀀싱 데이터를 가지고 첫 번째로 해야할일은 바로 시퀀스의 quality를 체크하는 일이다. 이를 quality control(QC)이라고 한다. 일반적으로 read의 길이와 수를 체크하게 되고 contaminating 시퀀스나 낮은 quality의 시퀀스가 있는지 찾아야한다. 이후에는 preprocessing 과정을 거치는데 이는 시퀀스의 퀄리티를 증가시키기 위한 과정이다. QC와 preprocessing 과정은 매우 중요하며, 이 과정이 제대로 되어야지 이후의 분석 결과를 신뢰할 수 있게 된다.

다음 스텝은 mapping이다. 또는 aligning이라고도 불린다. reads를 reference genome이나 reference transcriptome에 정렬하는 것을 뜻한다.  (reference genome은 표준게놈 또는 참조게놈 또는 참조 서열 등으로 불리며, HG19, HG38 등을 예로 들 수 있다.) 이렇게 함으로써 데이터의 뉴클레오타이드 서열을 de novo assembly 없이 연구할 수 있게 된다. 예를 들어 WES 데이터의 경우 WES 데이터의 read를 reference genome에 mapping 한다면, reference와 WES 데이터의 시퀀스 간의 다른 부분(variant)을 알아낼 수 있게 되고, 이 variant의 정확도는 mapping accuracy에 의존적이다.

따라서 그 다음에 해야할 일은 바로 mapping quality를 체크하는 일이다. 데이터의 특정 종류의 bias는 mapping step 이후 나타난다. mapping quality가 만족스럽지 않으면 이를 processing하는 과정을 거쳐야하는데 이를 post-alignment processing 이라고 부른다. 예를 들어, 중복된 mapped read를 제거하는 것을 의미한다. (이는 PCR의 인공물이다.) 이는 매우 중요하며 앞으로의 분석에 큰 영향을 미친다.

시퀀스가 참조 서열에 align 되었다면, 데이터를 experiment-specific한 방법으로 분석을 해야할 필요가 있다. 우리는 여기서 WES 데이터를 참조 서열에 mapping한 후에, variant 분석을 할 것이다. 즉, variant calling과 그 variant가 gene에 미치는 영향 (단백질의 변화, frame shift 등)을 알아볼 것이다. 이 과정에서 시퀀스를 참조 서열과 비교하고, 차이를 확인하여 이 차이가 유전자에 얼마나 큰 영향을 미칠지에 대해 분석하여야한다. 예를 들어, 이것이 snynonymous variant(mRNA가 생성하는 아미노산에 변화가 없는 변이)인 경우에는 이 영향이 미미할 것이다. 하지만 그 variant가 사이즈가 큰 deletion 이라면, 해당 시퀀스를 포함하는 유전자에 큰 영향을 줄 것이라고 예측해볼 수 있다. 분석과는 별개로 데이터를 visualization 해볼 수도 있다. mapped read 데이터를 visualising 함에 있어 가장 표준적인 tool 중 하나는 Genome Browser이다.

Gradient Boosting 의 기초

Gradient Boosting이란 도대체 무엇인지, 약한 Learner를 강화 시켜서 강한 학습기를 만드는 것이 무슨의미인지 직관적으로 이해할 수 있던 자료. https://www.youtube.com/watch?v=sRktKszFmSk Gradient Boosting = Residual Fitting 이다. 매우 Simple한 모델을 구축한 후, Residual에 Fitting한 모델을 만들고, 이 두 개를 결합한다. 그리고 결합된 모델에서 다시 Residual이 나오면 다시 이 Residual에 Fitting하는 모델을 만들어나가고 이를 계속 반복하여 최종 모델을 만들게 된다.

그리고 Stackoverflow에 누군가 Gradient Boosting에 대해 설명한 내용이다. GB는 residual fitting method이다. 아래 그림은 Regression에 GB 방법을 적용시킨 것이지만, Classification에도 GB 를 적용할 수 있다. 적절한 residual만 설정하면 똑같이 적용할 수 있다. 또한 GB 에 가장 인기 많은 베이스모델은 바로 Shallow한 Decision tree이다. decision tree는 다른 통계모형과는 다르게 가정이 적고, 범주형이든 연속형이든 제약없이 쉽게 만들 수 있는 모델이라고 설명하고 있다.

ClinVar DB 를 통한 질병 연관 변이 찾기

심근증(cardiomyopathy)과 관련있는 gene은 무엇이고, 그 gene의 어떤 위치의 variant가 pathogenic(질병의 원인이 될 수 있는)한지 알아보고자 한다. 이를 알아보기 위해서는 genotype-phenotype 데이터베이스를 이용하면되는데 이번 포스팅에서는 NCBI의 ClinVar 이라는 데이터베이스를 이용해보도록 하자. 우선 NCBI 홈페이지에 들어가 cardiomyopathy를 검색해본다.

그러면 아래와 같은 화면이 나오는데 ClinVar 을 클릭하여 들어간다.

이곳에서 왼쪽 메뉴바에서 Pathogenic을 클릭한다. 그러면 cardiomyopathy에 pathogenic한 gene의 variant를 찾을 수 있다.

또 많은 연구자들에 의해 보고된 pathgenic한 variant를 알아보고 싶다면, 아래와 같이 Multiple submitters를 체크한 후에 살펴보면 된다.

또 그 variant의 인해 분자적인 특성이 어떻게 변하는지 궁금할 수 있다. 이 때는 Molecular consequence에서 관심있는 분자적 특성에 체크를 하면된다. 예를 들어, missense(variant로 인해 mRNA의 코돈이 변해서 정상적인 경우와 다른 아미노산을 생성하는 분자적 변이)를 클릭하면 아래와 같은 화면을 볼 수 있다. 아래와 같이 다양한 pathogenic한 variant와 해당 gene을 알 수 있다. 해당 gene의 variant로 인해 missense가 일어나고 이로 인해 cardiomyopathy 질병을 일으키는 것이다.

크게보면 아래와 같이 variant를 나타내고 있다. 이렇게 variant를 표기하는 방법을 HGVS notation이라고 한다.

이를 통해 검색할 수도 있는데 예를 들어,

NM_213653.3(HFE2):c.959G>T (p.Gly320Val)

를 ClinVar 검색창에 검색해보자.

그러면 위와 같은 화면을 볼 수 있다. Links를 보면 해당 Variant가 UniProtKB, OMIM, dbSNP, 100 Genome Browser에 등록되어 있다는 것을 알 수 있다. 이를 통해 이 variant가 Multiple Submitters에 submit 되었다는 것을 알 수 있다. (Number of Submissions에 4가 표시되어 있다.)

ClinVar 데이터 다운로드 받기

https://www.ncbi.nlm.nih.gov/clinvar/에 접속 후 FTP => Go to the FTP site

여기서 vcf_GRCh37, vcf_GRCh38이 variant를 vcf 파일 형태로 보여주는 것이고, xml은 xml 파일 형태로 보여주는 것이다.

들어가면 다양한 버전이 존재하는데 가장 위에 있는 것을 받으면 된다. 

XGBoost 모델에서의 변수선택

/* 2017.9.3 by 3months */

원본글

https://machinelearningmastery.com/feature-importance-and-feature-selection-with-xgboost-in-python/

Gradient Boosting에서의 변수 중요도 (Feature Importance)

gradient boosting의 장점은 boosted tree가 완성된 후, feature의 중요도 스코어를 내는 것이 상대적으로 쉽다는 것이다.  어떠한 변수가 decision tree에서 "중요한 결정" 을 내리는데 사용된다면 그 변수의 중요도는 높을 것이다. 중요도는 변수별로 계산되며, 각각의 변수는 중요도 랭킹에 따라 정렬될 수 있다. 어떠한 변수의 중요도는 하나의 decision tree에서 그 변수로 인해 performance measure를 증가하는 양으로 계산된다. 이 때 performance measure는 purity(Gini index)가 사용된다. 최종적인 변수의 중요도는 모든 decision tree의 중요도의 평균값이 된다. boosted decision tree에서 변수의 중요도를 구하는 기술적인 방법은 https://www.amazon.com/dp/0387848576?tag=inspiredalgor-20 이 책의 page 367 에서 볼 수 있다고 한다. 또한 https://stats.stackexchange.com/questions/162162/relative-variable-importance-for-boosting 이 링크도 참고할만 하다.

Python XGBoost에서 변수의 중요도 구하기

get_xgb_imp 함수를 이용해 XGBoost의 변수 중요도를 구할 수 있다. 해당 함수의 파라미터로 트레이닝된 XGBoost의 모델과 feature의 이름을 리스트 데이터 타입으로 넣어주면 된다.

출처 - https://stackoverflow.com/questions/38212649/feature-importance-with-xgbclassifier

def get_xgb_imp(xgb, feat_names):
    from numpy import array
    imp_vals = xgb.booster().get_fscore()
    imp_dict = {feat_names[i]:float(imp_vals.get('f'+str(i),0.)) for i in range(len(feat_names))}
    total = array(imp_dict.values()).sum()
    return {k:v/total for k,v in imp_dict.items()}

예)

from xgboost import XGBClassifier

model = XGBClassifier()

# train_x와 train_y는 학습하고자 하는 트레이닝 셋의 설명변수와 종속변수의 ndarray 타입이다.
model.fit(train_x, train_y)

feat_names = [i for i in range(0,100)] # 변수가 100개인 경우 가정
feature_importance = get_xgb_imp(model,feat_names)

{0: 0.021963824289405683, 1: 0.020671834625323, 2: 0.041343669250646, 3: 0.023255813953488372, 4: 0.05297157622739018, 
5: 0.03488372093023256, 6: 0.025839793281653745, 7: 0.021963824289405683, 8: 0.04780361757105943, 9: 0.04392764857881137, 
10: 0.06847545219638243, 11: 0.023255813953488372, 12: 0.006459948320413436, 13: 0.012919896640826873, 
14: 0.012919896640826873, 15: 0.015503875968992248, 16: 0.0012919896640826874, 17: 0.003875968992248062, 
18: 0.006459948320413436, 19: 0.003875968992248062, 20: 0.002583979328165375, 21: 0.00904392764857881, 
22: 0.0012919896640826874, 23: 0.0012919896640826874, 24: 0.00516795865633075, 25: 0.0, 26: 0.0, 
27: 0.014211886304909561, 28: 0.002583979328165375, 29: 0.0012919896640826874, 30: 0.00516795865633075, 
31: 0.002583979328165375, 32: 0.007751937984496124, 33: 0.003875968992248062, 34: 0.00904392764857881, 
35: 0.0012919896640826874, 36: 0.012919896640826873, 37: 0.0103359173126615, 38: 0.015503875968992248, 
39: 0.003875968992248062, 40: 0.006459948320413436, 41: 0.006459948320413436, 42: 0.012919896640826873, 
43: 0.006459948320413436, 44: 0.0103359173126615, 45: 0.02454780361757106, 46: 0.025839793281653745, 
47: 0.021963824289405683, 48: 0.006459948320413436, 49: 0.0103359173126615, 50: 0.01808785529715762, 
51: 0.020671834625323 .....

* 결과는 위와 같이 dict 타입으로 반환된다.

변수 중요도 그래프 그리기

아래와 같이 matplotlib를 통해 변수 중요도의 barchart를 그릴 수 있다.

from matplotlib import pyplot
pyplot.bar(range(len(feature_importance)), feature_importance.values())
pyplot.show()

Keras를 통한 LSTM 구현

/* 2017.8.29 by. 3months */

- LSTM에 대한 개념 학습

http://blog.naver.com/kmkim1222/221069000196

https://www.youtube.com/watch?v=6niqTuYFZLQ&list=PL3FW7Lu3i5JvHM8ljYj-zLfQRF3EO8sYv

- 코드 및 데이터 출처  -https://github.com/Dataweekends/zero_to_deep_learning_udemy

LSTM 개념

LSTM은 기존 vanilla RNN을 개선한 모델로써 데이터의 long-term dependency를 학습하는데 효과적인 모델입니다. 또한 내부적으로 구성요소간에 additive, elementwise interaction을 통해 backpropagation 할 때, gradient explosion, vanishing 문제를 기존 RNN에 비해 매우 줄인 것이 중요한 개선점입니다.

코드 실습을 통해 LSTM을 학습하기 위해서는 어떠한 데이터 구조가 필요한지, 즉, INPUT은 어떤 모양이어야 하고, OUTPUT이 어떻게 생성되며, 네트워크 구조를 코드 상으로 어떻게 정의하는지 알아보겠습니다.

사용할 데이터

cansim-0800020-eng-6674700030567901031.csv

데이터 임포트

데이터를 임포트합니다. 데이터는 날짜별로 특정한 값이 있는 데이터입니다. Unadjusted이 Raw 값이고, Seaonally adjusted는 Season의 효과를 보정한 값입니다. 보통 주식과 같은 시계열 데이터는 Season에 의한 효과가 있는데, Seasonally adjusted는 이 Season에 의한 효과를 보정하여 전체적인 추세만 보고자할 때 사용하는 변수입니다. 우리가 학습할 변수는 Season에 의한 효과를 보정하지 않은 Unadjusted입니다.

import pandas as pd
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt

df = pd.read_csv('data/cansim-0800020-eng-6674700030567901031.csv',
                 skiprows=6, skipfooter=9,
                 engine='python')
df.head()

전처리

from pandas.tseries.offsets import MonthEnd

df['Adjustments'] = pd.to_datetime(df['Adjustments']) + MonthEnd(1)
df = df.set_index('Adjustments')
print(df.head())
df.plot()

우선 문자열로된 Adjustments 변수를 datetime 객체로 만들어주고, MonthEnd(1)을 더하면 해당 월의 맨 끝 날짜가 지정되게 됩니다. 그리고 이 날짜를 index로 만들어줍니다. 이후에 plot을 하면 날짜를 index로 해서 Unadjusted와 Seasonally adjusted 변수의 그래프를 그릴 수 있습니다.

인덱스가 0,1,2,3 ... 에서 날짜로 변했음을 알 수 있습니다. (가장 왼쪽열. 또한 index의 이름이 Adjustements가 되었음을 알 수 있음.)

우리가 학습할 변수는 파란색 선 Unadjusted 입니다.

트레이닝셋 테스트셋 SPLIT

2011/1/1 을 기준으로 트레이닝셋과 테스트셋으로 구분합니다. 그리고 학습과 테스트에 사용할 Unadjusted 변수만 남겨놓습니다.

split_date = pd.Timestamp('01-01-2011')
# 2011/1/1 까지의 데이터를 트레이닝셋.
# 그 이후 데이터를 테스트셋으로 한다.

train = df.loc[:split_date, ['Unadjusted']]
test = df.loc[split_date:, ['Unadjusted']]
# Feature는 Unadjusted 한 개

ax = train.plot()
test.plot(ax=ax)
plt.legend(['train', 'test'])

변수 Scaling

MinMax Scaling을 통해 변수를 Scaling합니다. MinMax Scaling : http://3months.tistory.com/167

그러면 아래와 같이 Scaling 된 것을 볼 수 있고, scaling의 결과는 2차원의 numpy ndarray 타입으로 변환이 되게 됩니다.

from sklearn.preprocessing import MinMaxScaler

sc = MinMaxScaler()

train_sc = sc.fit_transform(train)
test_sc = sc.transform(test)

train_sc

(240, 1)
[[ 0.01402033]
 [ 0.        ]
 [ 0.0704258 ]
 [ 0.09531795]
 [ 0.16362761]
 [ 0.13514108]
 [ 0.12395846]
 [ 0.12617398] 
......

Pandas Dataframe으로 변환

이를 다시 pandas dataframe 데이터 타입으로 변환해봅시다. Dataframe에서 작업을 해야지 Window를 만들기 유용하기 때문입니다. Window는 RNN을 트레이닝 하기위한 단위로 Window size는 Timestep과 같다고 생각하시면 됩니다.

train_sc_df = pd.DataFrame(train_sc, columns=['Scaled'], index=train.index)
test_sc_df = pd.DataFrame(test_sc, columns=['Scaled'], index=test.index)
train_sc_df.head()

pandas shift를 통해 Window 만들기

shift는 이전 정보 다음 row에서 다시 쓰기 위한 pandas의 함수입니다. 이를 통해 아래와 같이 과거의 값들을 shift_s 와 같은 형태로 저장할 수 있습니다. 과거값은 총 12개를 저장하며, timestep은 12개가 됩니다. 우리의 목적은 과거값 shift1~12를 통해 현재값 Scaled를 예측하는 것입니다.

for s in range(1, 13):
    train_sc_df['shift_{}'.format(s)] = train_sc_df['Scaled'].shift(s)
    test_sc_df['shift_{}'.format(s)] = test_sc_df['Scaled'].shift(s)

train_sc_df.head(13)

트레이닝셋과 테스트셋 만들기

dropna를 통해 NaN이 있는 데이터를 제거하고, shift_1 ~ shift_12는 X로 Scaled는 Y로 지정을 합니다.

X_train = train_sc_df.dropna().drop('Scaled', axis=1)
y_train = train_sc_df.dropna()[['Scaled']]

X_test = test_sc_df.dropna().drop('Scaled', axis=1)
y_test = test_sc_df.dropna()[['Scaled']]

최종 트레이닝셋과 테스트셋의 모양은 아래와 같습니다.

X_train.head()

y_train.head()

다시 ndarray로 변환하기

데이터구조를 다시 ndarray로 변환시켜줍니다. 실제 deep learning 모델의 트레이닝과 테스트로 사용되는 데이터는 일반적으로 numpy의 ndarray입니다. 아래와 같이 .values를 통해 ndarray 값을 얻을 수 있습니다. 트레이닝 데이터는 timestep=12, size = 228 입니다.

X_train = X_train.values
X_test= X_test.values

y_train = y_train.values
y_test = y_test.values

print(X_train.shape)
print(X_train)
print(y_train_shape)
print(y_train)

(228, 12)
[[ 0.20091289  0.13173822  0.11139526 ...,  0.0704258   0.          0.01402033]
 [ 0.03002688  0.20091289  0.13173822 ...,  0.09531795  0.0704258   0.        ]
 [ 0.01999285  0.03002688  0.20091289 ...,  0.16362761  0.09531795
   0.0704258 ]
 ..., 
 [ 0.79916654  0.81439355  0.86398323 ...,  0.92972161  0.71629034
   0.77368724]
 [ 0.80210057  0.79916654  0.81439355 ...,  0.59734863  0.92972161
   0.71629034]
 [ 0.81482896  0.80210057  0.79916654 ...,  0.53166512  0.59734863
   0.92972161]]
(228, 1)
[[ 0.03002688]
 [ 0.01999285]
 [ 0.06596369]
 [ 0.10983126]
 [ 0.14912986]
....

최종 트레이닝셋과 테스트셋의 X 만들기

이부분이 중요한데, keras에서는 RNN 계열의 모델을 트레이닝할 대 요구되는 데이터의 형식이 있습니다. 바로 3차원 데이터여야하며 각각의 차원은 (size, timestep, feature) 을 순서대로 나타내주어야하는 것입니다. 따라서 이 형태로 데이터를 reshape 해주어야합니다. 일반적인 MLP 모델에서는 size와 feature만 있기 때문에 2차원이지만, RNN에서는 "시간" 이라는 개념이 있기 때문에 차원이 한 차원 늘어나게 된 것입니다. 합리적인 데이터 구조라고 볼 수 있습니다.

아래 코드를 통해 Training set과 Test set의 X를 RNN 학습에 맞는 형태로 reshape 해줍니다.

X_train_t = X_train.reshape(X_train.shape[0], 12, 1)
X_test_t = X_test.reshape(X_test.shape[0], 12, 1)

print("최종 DATA")
print(X_train_t.shape)
print(X_train_t)
print(y_train)

최종 DATA
(228, 12, 1)
[[[ 0.20091289]
  [ 0.13173822]
  [ 0.11139526]
  ..., 
  [ 0.0704258 ]
  [ 0.        ]
  [ 0.01402033]]

 [[ 0.03002688]
  [ 0.20091289]
  [ 0.13173822]
  ..., 
  [ 0.09531795]
  [ 0.0704258 ]
  [ 0.        ]]

LSTM 모델 만들기

아래와 같이 keras를 통해 LSTM 모델을 만들 수 있습니다. input_shape=(timestep, feature)으로 만들어줍니다. size는 모델 설계시에는 중요하지 않으므로, feature, timestep만 모델에 알려주면 됩니다. 또 예측하고자하는 target의 갯수가 1이므로 마지막에 Dense(1)을 하나 추가해줍니다. 또한 실제 연속적인 값을 예측하는 것이기 때문에 loss function은 mean squared error가 됩니다. 또한 일반적으로 optimizer는 adam을 자주 사용합니다.

from keras.layers import LSTM
from keras.models import Sequential
from keras.layers import Dense
import keras.backend as K
from keras.callbacks import EarlyStopping

K.clear_session()
model = Sequential() # Sequeatial Model
model.add(LSTM(20, input_shape=(12, 1))) # (timestep, feature)
model.add(Dense(1)) # output = 1
model.compile(loss='mean_squared_error', optimizer='adam')

model.summary()

모델 Fitting

모델을 Fitting한다는 것은 Training data set으로 optimization 과정을 통해 모델의 weight를 찾는 것입니다. early stopping 객체를 이용해 epoch마다 early stopping을 체크합니다. 아래와 같은 결과를 내면서 트레이닝이 잘 되는 것을 확인할 수 있습니다.

early_stop = EarlyStopping(monitor='loss', patience=1, verbose=1)

model.fit(X_train_t, y_train, epochs=100,
          batch_size=30, verbose=1, callbacks=[early_stop])

Epoch 1/100
228/228 [==============================] - 0s - loss: 0.3396     
Epoch 2/100
228/228 [==============================] - 0s - loss: 0.2958     
Epoch 3/100
228/228 [==============================] - 0s - loss: 0.2551     
Epoch 4/100
228/228 [==============================] - 0s - loss: 0.2175   
...

학습된 모델을 통해 테스트셋 Test 하기

테스트셋은 2011/1/1 이후의 데이터를 나타냅니다. 트레이닝에는 1991년~ 2010년 까지 데이터가 이용되었기 때문에, 이 실습은 1991년~2010년 데이터를 통해 2011년 이후의 데이터를 예측하는 것으로 볼 수 있습니다.

테스트셋의 모양은 아래와 같이 생겼습니다. 당연한 이야기지만 트레이닝셋과 구조가 같습니다.

print(X_test_t)

[[[ 1.06265011]
  [ 0.87180554]
  [ 0.84048091]
  [ 0.86220767]
  [ 0.88363094]
  [ 0.89302107]
  [ 0.92552046]
  [ 0.89993326]
  [ 0.83505683]
  [ 0.77259579]
  [ 0.56926634]
  [ 0.61423187]]


....

model.predict를 통해 테스트셋의 X에 대한 예측값 y_hat을 얻을 수 있습니다.

y_pred = model.predict(X_test_t)
print(y_pred)

[[ 0.82410765]
 [ 0.85221326]
 [ 0.88795143]
 [ 0.90008551]
 [ 0.90281236]
 [ 0.90063977]
 [ 0.89379823]
 [ 0.88957471]
 [ 0.88779473]

...

Python Scaling 방법 정리

/* 2017.8.27 by. 3months */

머신러닝에서 데이터를 모델에 트레이닝할 때, 일반적으로 데이터를 Scaling 해서 사용한다. Scaling을 통해 데이터의 Scale을 맞추면 Weight의 scale도 일관성 있도록 맞출 수 있을 것이다.

MinMax Scaling은 최댓값 = 1, 최솟값 = 0으로 하여, 그 사에 값들이 있도록 하는 방법이다.

즉 데이터 X의 scaling된 값 Xscale은 아래와 같다.

Xscale = (X-Xmin) / (Xmax-Xmin)

from sklearn.preprocessing import MinMaxScaler
a = [[10], [9], [8], [6], [2]]
scaler = MinMaxScaler(feature_range=(0,1))
a = scaler.fit_transform(a)
print(a)

결과

[[ 1.   ]
 [ 0.875]
 [ 0.75 ]
 [ 0.5  ]
 [ 0.   ]]

StandardScaling은 평균=0과 표준편차=1이 되도록 Scaling 하는 방법이다. 이는 표준화라고도 부른다.

Xscale = (X-X_mean) / Xstd

from sklearn.preprocessing import StandardScaler
a = [[10], [9], [8], [6], [2]]
scaler = StandardScaler()
a = scaler.fit_transform(a)
print(a)
print(a.mean())
print(a.std())

결과

[[ 1.06066017]
 [ 0.70710678]
 [ 0.35355339]
 [-0.35355339]
 [-1.76776695]]
0.0
1.0