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

필자가 설치한 환경은 운영체제 우분투 14.04 LTS Server 버전이고 GPU는 1080TI를 사용한다.

1. CUDA 설치

https://developer.nvidia.com/cuda-downloads

이렇게 선택후 아래 나오는 대로 한다. (deb파일 다운로드후 입력하라는 커맨드 입력하면 설치됨)

1. `sudo dpkg -i cuda-repo-ubuntu1404_8.0.61-1_amd64.deb`
2. `sudo apt-get update`
3. `sudo apt-get install cuda`

2. cuDNN설치

https://developer.nvidia.com/rdp/cudnn-download

해당 링크에서 로그인이 필요하다. 로그인 후에,,

Download cuDNN v7.0 (August 3, 2017), for CUDA 8.0 <= 이거 설치

이 때, CUDA 버전에 맞는 cuDNN을 설치하여야 한다. 저는 CUDA 8.0 이므로 위 링크 cuDNN을 설치.

* wget 으로 해당링크 다운받으려하니 404 forbidden 오류가 뜬다.

* FTP로 그냥 압축파일을 서버로 옮겼다.

tar xvzf cudnn-8.0-linux-x64-v7.tgz

sudo cp -P cuda/include/cudnn.h /usr/local/cuda/include

sudo cp -P cuda/lib64/libcudnn* /usr/local/cuda/lib64

sudo chmod a+r /usr/local/cuda/include/cudnn.h /usr/local/cuda/lib64/libcudnn*

sudo apt-get install libcupti-dev

또한 ~/.bashrc에 cuda, cudnn 라이브러리 경로를 적어두어야한다.

vi ~/.bashrc

로 파일을 열어서 맨 밑에 아래 텍스트 추가.

export PATH="/home/jyahn/anaconda2/bin:$PATH"
export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/usr/local/cuda/lib64:/usr/local/cuda/extras/CUPTI/lib64"
export CUDA_HOME=/usr/local/cuda

tar xvzf cudnn-8.0-linux-x64-v5.1-ga.tgz
sudo cp -P cuda/include/cudnn.h /usr/local/cuda/include
sudo cp -P cuda/lib64/libcudnn* /usr/local/cuda/lib64
sudo chmod a+r /usr/local/cuda/include/cudnn.h /usr/local/cuda/lib64/libcudnn*

tar xvzf cudnn-8.0-linux-x64-v5.1-ga.tgz
sudo cp -P cuda/include/cudnn.h /usr/local/cuda/include
sudo cp -P cuda/lib64/libcudnn* /usr/local/cuda/lib64
sudo chmod a+r /usr/local/cuda/include/cudnn.h /usr/local/cuda/lib64/libcudnn*

tar xvzf cudnn-8.0-linux-x64-v5.1-ga.tgz
sudo cp -P cuda/include/cudnn.h /usr/local/cuda/include
sudo cp -P cuda/lib64/libcudnn* /usr/local/cuda/lib64
sudo chmod a+r /usr/local/cuda/include/cudnn.h /usr/local/cuda/lib64/libcudnn*

tar xvzf cudnn-8.0-linux-x64-v5.1-ga.tgz
sudo cp -P cuda/include/cudnn.h /usr/local/cuda/include
sudo cp -P cuda/lib64/libcudnn* /usr/local/cuda/lib64
sudo chmod a+r /usr/local/cuda/include/cudnn.h /usr/local/cuda/lib64/libcudnn*

$ sudo apt-get install libcupti-dev

3. tensorflow GPU 설치

아래 링크에 따라 virtual env를 구성한다.

https://github.com/Dataweekends/zero_to_deep_learning_udemy

이후,

source activate ztdl

conda install tensorflow-gpu

4. 설치되었나 확인

bash에서

python

>> import tensorflow as tf

>> sess = tf.Session()

아래와 같이 GPU 이름이 나오면 성공.

name: GeForce GTX 1080 Ti
major: 6 minor: 1 memoryClockRate (GHz) 1.582
pciBusID 0000:02:00.0
Total memory: 10.91GiB
Free memory: 10.75GiB
2017-08-30 17:41:13.373074: W tensorflow/stream_executor/cuda/cuda_driver.cc:523] A non-primary context 0x2d63250 exists before initializing the StreamExecutor. We haven't verified StreamExecutor works with that.
2017-08-30 17:41:13.374815: I tensorflow/core/common_runtime/gpu/gpu_device.cc:940] Found device 1 with properties:
name: GeForce GTX 1080 Ti
major: 6 minor: 1 memoryClockRate (GHz) 1.582
pciBusID 0000:03:00.0
Total memory: 10.91GiB
Free memory: 10.75GiB
2017-08-30 17:41:13.695077: W tensorflow/stream_executor/cuda/cuda_driver.cc:523] A non-primary context 0x2d67600 exists before initializing the StreamExecutor. We haven't verified StreamExecutor works with that.
2017-08-30 17:41:13.696140: I tensorflow/core/common_runtime/gpu/gpu_device.cc:940] Found device 2 with properties:
name: GeForce GTX 1080 Ti
major: 6 minor: 1 memoryClockRate (GHz) 1.582
pciBusID 0000:82:00.0
Total memory: 10.91GiB
Free memory: 10.75GiB
2017-08-30 17:41:14.016589: W tensorflow/stream_executor/cuda/cuda_driver.cc:523] A non-primary context 0x2d6b9b0 exists before initializing the StreamExecutor. We haven't verified StreamExecutor works with that.
2017-08-30 17:41:14.017710: I tensorflow/core/common_runtime/gpu/gpu_device.cc:940] Found device 3 with properties:
name: GeForce GTX 1080 Ti
major: 6 minor: 1 memoryClockRate (GHz) 1.582
pciBusID 0000:83:00.0
Total memory: 10.91GiB
Free memory: 10.75GiB

참고

http://goodtogreate.tistory.com/entry/TensorFlow-GPU-%EB%B2%84%EC%A0%84-%EC%9A%B0%EB%B6%84%ED%88%AC-1604%EC%97%90-%EC%84%A4%EC%B9%98-%ED%95%98%EA%B8%B0

CIGAR Format

출처 - https://genome.sph.umich.edu/wiki/SAM

/* 2017.8.13 by 3months. */

CIAGR Format은 어떠한 Sequence A가 reference genome에 align되어있을 때 이를 나타내 주는 포맷이다.

예를 들어, 아래와 같은 reference sequence와 Read가 있다고 하자. 이 Read를 reference sequence에 align 해보면

RefPos:     1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T  G  A  A  C  T  G  A  C  T  A  A  C
Read: ACTAGAATGGCT

이와 같이 align 된다.

RefPos:     1  2  3  4  5  6  7     8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T     G  A  A  C  T  G  A  C  T  A  A  C
Read:                   A  C  T  A  G  A  A     T  G  G  C  T

이 때 이것을 나타내는 CIGAR Format은

POS: 5
CIGAR: 3M1I3M1D5M

이렇게 되는데 reference의 5번 position부터 시작해서 3개의 Matching, 1개의 insertion, 다시 3개 matching, 1개 deletion, 5개 matching 이라는 뜻이다. ref의 14번 position에 read와 시퀀스가 다르지만, align이 이렇게 되기 때문에 matching으로 쳐준다.

Lastz를 통한 CIGAR Format 생성

lastz_test.fa

lastz_test2.fa

의도적으로 lastz_test fasta 파일의 sequence에다가 4개의 sequence를 insert하고 2개의 sequence를 delete해서 lastz_test2.fa 파일을 생성하였다.

그 다음

lastz lastz_test.fa lastz_test2.fa --format=general:cigar > cigar.txt

를 실행 후, cigar.txt 파일을 열어보면

#cigar
208M4I163M2D138M

위와 같은 결과가 나왔음을 볼 수 있다. 결과를 해석하면 208개 matching 4개 insertion, 163개 matching 2개 deletion, 138개 matchin 이다.

2017.10.8 추가

lastz 사용해서 cigar을 얻을 때 다른 옵션을 주면 substitution도 고려가 가능하다.

lastz_test.fa

lastz_test2.fa

Lastz 실행

/* 2017.8.13 by. 3months */

Lastz를 설치하였으면 간단하게 결과를 뽑아보자.

두 개의 fasta 파일을 폴더에 붙여넣고,

해당 폴더에서 아래 명령어를 입력해보자.

lastz lastz_test.fa lastz_test2.fa

결과

그러면 콘솔창에 이런 결과가 나온다.

또는 아래와 같이 옵션을 주어서 실행하면 maf 파일로 내보낼 수 있다.

lastz lastz_test.fa lastz_test2.fa --format=maf > lastz_result.maf

lastz_result.maf 파일을 열어보면 아래와 같이 되어있다.  (이것이 무엇을 뜻하는지는 다음에 알아보자..)

##maf version=1 scoring=lastz.v1.04.00
# lastz.v1.04.00 --format=maf
#
# hsp_threshold      = 3000
# gapped_threshold   = 3000
# x_drop             = 910
# y_drop             = 9400
# gap_open_penalty   = 400
# gap_extend_penalty = 30
#        A    C    G    T
#   A   91 -114  -31 -123
#   C -114  100 -125  -31
#   G  -31 -125  100 -114
#   T -123  -31 -114   91

a score=47093
s ref 0 511 + 511 GAAGAGCAGCACACTCAGATCTTCTTACACCATCTCTGGGGGCCGACTGGCAGAGAAGGCTAGGGACTTGGCCAATATCATGCAGCCTTCATGCCTTATTCGGGAAAGGTGCTCCAAATAGGAGGAATATTTGGAAGAAATCCCAATAGGGATTTGGAGACATGCCATTTGGGAATGCAAACCCAGCCCCTCCT--------AGTGCTCCACCCTG----CCCATGTCCCCATTGAACCGACAACCTGCAGAGAGAGCAGGATGCCCACAGCAAGCCCCAAGGACAAAGGCCACAGAGTCAGGGCACCCAAATCCTAACCTAACCTTGCCTGTGCTTCCTGTGTACTGTGGCCATTTTACTTTCCCCCTGTTCCCTCCTCTGTAAATGAGGCCCTTGGACTGCATCAGTGGTTTTTACACTGAGCTCCCTGGAGCTCTGGGGGCCTCAGGGACTATTGTGGAGGCCTTGCTGCATTTGAGAGTGGTGTCTTGTACCCTACTCAACCCATTTTTTCTGTTTTAC

s ref 0 521 + 521 GAAGAGCAGCACACTCAGATCTTCTTACACCATCTCTGGGGGCCGACTGGCAGAGAAGGCTAGGGACTTGGCCAATATCATGCAGCCTTCATGCCTTATTCGGGAAAGGTGCTCCAAATAGGAGGAATATTTGGAAGAAATCCCAATAGGGATTTGGAGACATGCCATTTGGGAATGCAAACCCAGCCCCTCCTAGAGAGAGAGTGCTCCACCCTGCCCCCCCATGTCCCCATTGAACCGACAACCTGCAGAGAGAGCAGGATGCCCACAGCAAGCCCCAAGGACAAAGGCCACAGAGTCAGGGCACCCAAATCCTAACCTAACCTTGCCTGTGCTTCCTGTGTACTGTGGCCATTTTACTTTCCCCCTGTTCCCTCCTCTGT--ATGAGGCCCTTGGACTGCATCAGTGGTTTTTACACTGAGCTCCCTGGAGCTCTGGGGGCCTCAGGGACTATTGTGGAGGCCTTGCTGCATTTGAGAGTGGTGTCTTGTACCCTACTCAACCCATTTTTTCTGTTTTAC

아래 링크에서 lastz 실행과 관련된 다양한 예제를 볼 수 있다.

http://www.bx.psu.edu/miller_lab/dist/README.lastz-1.02.00/README.lastz-1.02.00a.html#examples

DNA 결합 단백질

Regulatory Genomics 관련 논문을 읽던 중 DNA binding protein을 보고 몰라서 공부하였다.

DNA => RNA => Protein이 되는 과정에서 DNA=>RNA로 변하는 transcription 과정 중 DNA가 RNA로 변하는 것의 시작은 전사개시복합체(transcription initiation complex)가 Promoter 부분에 붙으면서 시작된다. 이 복합체 단백질 중 하나인 RNA중합효소2가 유전자를 전사헤 pre-mRNA를 만든다. 그 다음에 RNA Processing(5'캡 생성, 폴리A 꼬리 생성, 인트론 splicing)을 통해 mature mRNA가 생성된다.

이 과정에서 RNA 중합효소가 전사를 개시하려면 전사인자(transcription factor)라는 단백질이 필요하다. 이들은 보편전사인자(general transcription factor)라고도 불린다. 하지만 이것만 가지고는 전사개시속도가 매우 낮으며 특수전사인자라는 조절 단백질의 작용이 필요하다.

promoter 주변에는 DNA 시퀀스에는 근거리 조절요소가 있다. 그리고 유전자에 멀리 떨어져 있는 원거리 조절요소인 인핸서가 있다. (한 개의 인핸서는 하나의 유전자에만 작용하며, 한 유전자는 여러개의 인핸서를 갖을 수 있다. 인핸서에 붙는 단백질들은 특수전사인자라고 불리며 전사인자 종류중 하나라고 볼 수 있다 )

진핵생물에서 유전자 발현율은 인핸서의 조절요소에 활성자(activator) 또는 억제자(repressor)가 결합함에 따라 증가되거나 감소된다. 결합되는 부분은 DNA Binding site라고 불리는듯하다. 즉 인핸서의 조절 요소에 특정한 단백질이 붙어서 그 인핸서의 표적 유전자의 발현을 조절하는 것이다.  활성자 단백질은 RNA 중합효소를 promoter로 recruting한다. 그래서 전사속도가 빨라진다. 반면 억제자 단백질은 RNA 중합효소가 이동하지 못하게 하여 전사가 될수 없게 만든다.

이렇게 gene expression 조절을 수행하는 여러가지 형태의 단백질(RNA 중합효소, 활성자, 억제자 단백질)을 DNA 결합 단백질 (DNA Binding protein)이라고 하며 이들은 sequence specific하다. 즉, DNA Binding protein을 코딩하는 시퀀스에 문제가 있을 경우 해당 protein의 표적 유전자의 발현에 이상이 생길 수 있다. 그래서 DNA Binding protein의 sequence specificity를 안다면 표적 유전자를 알 수 있게 되고, DNA Binding protein 시퀀스의 disease causing variant를 파악하는데 도움이 된다.

이러한 Transcription factor 등과 같은 단백질에 결합하는 DNA 시퀀스를 알아내는 방법으로 CHIP-seq 방법이 있다.

참고 - 캠밸 생명과학

유전체학을 위한 신경망 모델 (Neural Networks for Genomics)

원본글

http://blog.qure.ai/notes/neural-networks-for-genomics

위 글을 번역한 포스팅입니다.

Introduction

딥러닝을 처음 유전체학에 적용하려고 한다면, 최신으로 이루어지고 있는 연구(state-of-the-art)가 무엇인지, 연구자들이 딥러닝을 통해 어떠한 문제를 해결하고자 하고, 이를 어떻게 접근하고 있는지 알아야할 필요가 있다. 이 포스팅은 유전체학에서 유명한 네트워크 구조를 소개하고, 이 네트워크를을 트레이닝할 때 쓰이는 데이터의 타입, 그리고 모델을 통해 최종적으로 예측한 결과가 어떻게 나오는지를 요약한다.

시퀀싱 기술의 발달, 그리고 1000-genome project, encode, geo와 같은 public dataset의 범람에도 불구하고 아직까지 genotype-phenotype 구분, 그리고 시퀀스를 통해 질병을 예측하는 것을 완벽하게 해내지 못하고 있다. Bredan Frey의 이 talk는 딥러닝과 유전체학의 컨텍스트상에서 왜 더 많은 genome을 시퀀싱 하는 것이 반드시 해답이 아니라는 것을 설명하고 있다. genome은 매우 복잡하고, 서로 상호작용하는 수많은 정보 레이어를 갖고 있다. 현재의 대부분의 접근법은 시퀀스를 통해 직접적으로 표현형을 추측하는 것이 아니라, 시퀀스의 일부를 해석하는 시스템을 구축하는 것이다. 아래에 딥러닝이 유전체학이나 전사체학에 적용된 예를 소개한다.

얕은 fully connected networks 를 이용한 초기 연구

초기 연구중 몇몇은 주성분분석을 통해 차원을 축소한 gene expression 데이터를 single layer fully connected network에 적용한 것이다. 이러한 초기 연구는 gene expression을 통해 tumor type을 구분하거나, tumor grade를 예측하거나, 환자의 생존 여부를 예측하는데에 사용되었다. 이러한 연구의 발달은 가장 예측력이 높은 gene의 subset을 찾거나, gene의 signature를 찾는 연구의 발달을 가져왔다. 이와 비슷하게 인공신경망은 microRNA의 발현 패턴을 통해 대장암의 등급을 예측하거나, 암 상태를 예측할 수 있는 microRNA를 찾아내는 연구에 사용되었다.

Autoencoder를 이용한 차원 축소와 특징값 추출

Autoencoder는 유전체학에서 gene expression 데이터로부터 feature space를 줄이거나 유용한 feature를 찾아내는데 사용되어왔다. 한 예는 autoencoder를 사용해 labling 되지 않은 gene expression data로부터 정교한 feature representation을 구현하고, 이를 통해 tumor를 예측하는 classifier를 구현한 이 페이퍼이다. 더 최근 연구는 denoising autoencoder를 이용해 유방암 gene expression data에서 feature extraction을 한 연구이다. 다음은 ADAGE(paper and repository)인데, 박테리아 gene expression 데이터에 비슷한 접근을 한 연구이다. autoencoder를 통해 추출한 유용한 feature들은 암 분류, 에스트로겐 수용체의 상태, 환자의 생존 예측에 활용될 수 있다. 그러면 하나의 autoencoder로 부터 밝혀진 gene expression profile들을 통해 어떠한 조직에도 일반적으로 적용할 수 있는 feature를 밝혀낼 수 있을까? gene expression 데이터 말고 DNA sequence에 autoencoder를 적용하여 유용한 feature를 추출할 수 있을까?

Deep learning을 통해 gene expression, transcript expression level을 예측하기

 유전자의 발현은 공동으로 이루어지는 경우가 많기 때문에, 서로 다른 유전자의 발현 정도는 상관관계가 매우 높다. 이는 특정 subset의 유전자 발현 데이터만 가지고도 나머지 gene의 발현 정도를 예측해볼 수 있다는 것을 뜻한다. 이는 잠재적으로 유전자 발현 profiling에서 비용과 복잡성을 줄일 수 있다. D-GEX라고 불리는 방법이 개발되었는데 이는 퍼블릭 데이터인 CMAP을 통해 트레이닝된 심층신경망 모델로써, 1000개의 gene expression 데이터를 통해 나머지 gene 의 expression 정도를 예측한다. 최근의 topcoder challenge도 이와 비슷한 문제가 올라와있다. 이와 비슷하지만 더 hard한 태스크는 exon이나 transcript의 expression level을 DNA sequence data로 부터 예측하는 것이다. expression level은 sequence로만 결정되는 것이아니라, cellular context 에도 의존한다. ‘Deep learning of the tissue-regulated splicing code’ 라는 타이틀의 연구는 exon 주변의 dna sequence를 통해 그 exon이 splice in 되는 확률 (PSI)를 예측한다. 직접 손으로 만든 genomics feature들(이 genomics feature들은 특정한 mouse의 cell type에서 splicing pattern을 예측할 수 있는 것들이다.) 이 모델을 트레이닝하는데 사용되었다. autoencoder를 통해 이 feature의 차원을 축소하고, cell type을 나타내는 추가적인 input을 통해 multilayer fully connected network를 트레이닝하였다. 이 방법을 통해 저자들은 tool을 개발하고 검증하였는데, 이것을 통해 SNV가 splicing애 미치는 effect를 scoring할 수 있다.

Epigenomics를 위한 Convolutional Networks

최근 연구의 많은 부분이 바로 전사인자 결합 위치나, 인핸서 지역, 시퀀스를 통한 크로마틴 accessibility 예측과 같은 epigenomic의 문제에 CNN을 이용하는 것이다. DeepBind는 주어진 sequence를 통해 DNA나 RNA 결한 proteins의 specificity를 예측할 수 있는 방법이다. 이것을 하기 위해 많은 수의 protein binding micro-array chip-seq 데이터를 통해 CNN이 training 되었다. convolutional stage에는 one-hot encoding된 시퀀스에서 'motif detector' matiireces의 set을 scan한다. learning된 filter는 유전쳏학에서 DNA seuqnece motif를 설명하고 이해하는데 사용된다. 네트워크는 알려진 motif와 알려지지 않은 motif까지도 learning할 수 있다. DeepMotif는 비슷한 계열의 더 깊은 모델로써, 시퀀스에 따라 transcription factor binding(yes or no)를 분류하고, motif를 추출하는데 더 중점을 둔 모델이다. DeepSea는 ENCODE의 epigenomics 데이터로부터 트레이닝되었고, 단일염기의 mutation이 transcription factor binding과 DNAse sensitivity에 미치는 영향을 예측한다. 이 모델의 주요 기능은 다양한 규모의 시퀀스 기능을 학습하는 계층적 아키텍처의 사용, 광범위한 시퀀스 컨텍스트 스캔 기능 및 예측 기능을 공유하는 다양한 염색질 요소의 멀티 태스크 공동 학습이다. 다른 epigenomics에 사용되는 CNN의 예는 Basset인데, DNA 시퀀스로부터 chromatin accessibility code를 예측한다. DeepCpG는 DNA metylation을 모델링하였고, DEEP은 enhancer나 transcription factor binding으로 인해 gene의 transcription이 증가하는 DNA의 region을 예측하는 ensemble framework이다. 또한 chip-seq 실험에서 noise를 감소시키거나 signal을 증폭하기 위한 CNN-based 방법도 있다.

언어처리로부터 영감을 얻은 모델들

언어처리와 genome해석은 어느정도 비슷한 면이있다. 이는 언어처리에 사용되는 방법이 genomics를 이해하는데 유용할 수 있다는 것을 시사한다.

Recurent Neural Networks

RNN은 긴 범위를 갖는 DNA 시퀀스의 상호작용을 발견하기 위해 사용되어왔다. DanQ는 CNN과 bi-directional LSTM(long-short term memory)의 하이브리드로써, convolution layer는 regulatory motif를 발견하고, recurrent layer는 regulatory 'grammer'를 이해하기 위해 motif들 간의 long-term dependency를 발견하는데 사용된다. 다른 예는 DNA-level splice junction prediction인데 RNN이 splice junction이나 intro과 exon의 경계를 발견하기위해 training이 된다. 또한 바이러스 gemoe에서 proten coding region을 발견하는 것도 있다.

Word embedding이나 word2vec 같은 모델

‘Gene2vec: Neural word embeddings of genetic data’ 라는 타이틀을 갖고있는 한 프로젝트는 구글의 오리지날 word2vec을 genome 시퀀스로 구현한 것이다. (genome을 27bp의 'word'로 splitting한다.) 실제로 이것이 인간 genom에 실용적으로 적용될 수 있는지는 흥미로운 주제일 것이다. 비슷한 word2vec 스타일의 모델은 gene expression data로도 트레이닝이 가능하다.

기타 링크

dna binding protein motif

http://2013.igem.org/wiki/index.php?title=Team:XMU_Software/Project/promoter&oldid=359507

Regulatory Genomics 개념

regulatory genomics 관련 자료를 찾아보다 아래 링크를 찾아 번역하여 정리하였습니다.

https://biology.stackexchange.com/questions/17810/what-does-regulatory-genomics-mean

Genome은 크게 두 파트로 나뉘어진다.

1. 단백질을 코딩하는 파트 (coding region)

2. 단백질을 코딩하지 않는 파트 (non-coding region)

이 때 단백질을 코딩하지 않는 두 번째 파트는 또 다시 두 개의 클래스로 나뉘어진다.

2-1. transcribe가 되어 생물학적인 활동을 하는 파트 (long noncoding RNA, miRNA, competing endogenous RNA 등)

2-2. transcribe가 안되는 부분 (unscribed non-coding region)

이 중에서 2-2번 unscribed noncoding region은 아마도 1번파트와 2-1번 파트에 regulation 작용을 할 것으로 생각된다. (transcription factor나 transcriptional coactivators/corepressors에 binding 함으로써)

※ transcriptional coactivators/corepressors는 gene의 시작부분 (promoter) 또는 멀리 떨어진 부분(enhancer)에서 해당 gene이 expression이 될지 안될지를 조절한다. enhancer은 gene과 멀리 떨어져있지만 3차원 상에서는 붙어있을 수 있기 때문에 gene expression을 조절할 수 있다. (genome은 세포 안에서 접혀져 있기 때문이다.)

Regulatory Genomics는 바로 genomics "features"라고도 불리는 이 unscribed non-coding region에 대한 연구이다. 어떻게 그들이 gene regulation을 하는지 알아내는 것이다. 이 분야의 연구의 대표적인 예는 ENCODE project의 파트로써 publish된 페이퍼들이다.

ENCODE 프로젝트 링크

http://www.nature.com/encode/#/threads

ENCODE 프로젝트 소개 (홈페이지 내용 번역)

ENCODE(Encyclopedia of DNA Elements) 프로젝트는 National Human Genome Research Institute가 후원하는 프로젝트로써 genome 상에서의 transcription, transcription factor association, chromatin structure, histone modification을 밝혀내기 위한 프로젝트이다. 이러한 genome 상에서의 기능적인 요소들을 식별함으로써 현재 인간 genome의 80%의 부분이 최소 1개 이상의 biochemical function 을 한다는 것 알아내었다. 이러한 functional annotation에 관한 광활한 자원들은 genome과 유전자의 regulation과 oragnization에 있어 새로운 통찰력을 제공하고 있다.

ENCODE 관련 한글로 정리된 블로그

http://blog.daum.net/kimuks/7535222

https://madscientist.wordpress.com/2012/09/18/%EC%98%A4%EB%8A%98%EC%9D%98-%EB%85%BC%EC%9D%BD%EB%82%A8-%EC%97%94%EC%BD%94%EB%93%9C%EB%A5%BC-%EB%94%94%EC%BD%94%EB%93%9C%ED%95%98%EA%B8%B0/

Regulatory Genomics 및 Epigenomics 관련 강의 자료

https://simons.berkeley.edu/talks/regulatory-genomics-epigenomics