IT Note/Data&AI

COVID-19의 기원을 Data Science로 예측해보자!

GRuuuuu 2020. 7. 28. 15:08

원문 : 호롤리한하루/COVID-19의 기원을 Data Science로 예측해보자!

Overview

1.해당 포스팅은 다음 문서를 바탕으로 재가공한 문서입니다.
-> Using Data Science to Predict the Origin of SARS-CoV-2 Coronavirus from Public Genome Data
원 포스팅이 20년 2월에 작성되었으니 사실과 다른 내용이 포함되어 있을 수 있는 점 양해부탁드립니다.

2.본문에서 사용한 코드는 여기 있습니다.

SARS-CoV-2 코로나 바이러스는 어디에서 시작되었을까요?

전문가들은 박쥐가 바이러스의 근원지라고 추측하고 있습니다. 어째서 그런 결론이 나왔을까요?

전통적으로 의학자들은 BLAST같은 생물정보학(bioinformatic) 툴을 사용합니다. 이번 포스팅에서는 코로나 바이러스의 게놈을 데이터사이언스 관점에서 분석하는 방법을 말씀드리겠습니다.

BLAST: 단백질의 아미노산 서열 또는 DNA 및 / 또는 RNA 서열의 뉴클레오티드와 같은 1 차 생물학적 서열 정보를 비교하기위한 알고리즘 및 프로그램

코로나 바이러스에 감염된 다른 생물들 ( 닭, 소, 오리, 박쥐 )의 게놈 순서를 비교하여 COVID-19의 기원을 예측해 볼 것입니다.

Datasets

코로나바이러스의 게놈정보는 U.S National Library of Medicine website에 공개되어있습니다.

위 링크를 클릭하면 다음과 같은 사진을 보실 수 있습니다.
image

사이트에서는 두가지 타입의 바이러스 게놈을 확인할 수 있습니다.
nucleotide와 protein타입입니다. 해당 포스팅에서는 닭, 박쥐, 소, 오리 4가지 동물의 nucleotide 타입을 분석해보겠습니다. 물론 protein타입도 아래 nucleotide분석방법과 동일하게 진행하실 수 있습니다.

Data Collection

분석의 목표는 바이러스의 근원이 대략 어디서 왔는가를 찾는 것입니다.
제대로 하려면 모든 동물의 바이러스를 스캔하여 비교해야하지만 컴퓨팅리소스의 한계로 4가지 동물만 가지고 분석을 진행하도록 하겠습니다.

첫째로, 코로나바이러스의 nucleotide타입을 다운로드 받겠습니다.
image
링크로 들어가게 되면 자동으로 중증 급성 호흡기 증후군 코로나 바이러스(Severe acute respiratory syndrome coronavirus 2)가 선택되어 결과가 출력됩니다.

좌측 메뉴에 "Nucleotide Completeness"에서 Complete를 선택한 후, 우측상단의 Download 버튼을 눌러서 게놈순서를 다운로드 받습니다. (.fasta확장자)

2020.07.23 기준 총 7194개의 염기서열을 다운로드 받을 수 있습니다. 200개정도만 선택해서 다운로드 받겠습니다.

image

Nucleotide > Download Selected Records > use default 순으로 선택
image

메모장으로 열어보면 '>'뒤에 바이러스의 정보가 나오고 뒤에 염기서열이 나와있습니다.
image
이건 사람을 감염시킨 코로나바이러스의 염기서열이니 이제 각 동물들의 코로나바이러스 염기서열을 받아보겠습니다.

Virus탭에서 COVID-19 coronavirus 옵션을 체크해제하고 Coronaviridae를 검색해주세요. Coronaviridae는 동물을 포함한 모든 코로나바이러스를 지칭하는 학명입니다.
image

먼저 닭부터 시작하겠습니다. 닭의 학명은 Gallus gallus입니다.
Host탭에 gallus gallus를 검색해서 위와 동일하게 데이터를 다운로드 받아줍니다.
image

받아야 할 동물들의 학명은 다음과 같습니다. 각 데이터가 전부 400개 이하라 얼마 안되니 전부 다운로드 받아줍니다.
닭 : gallus gallus
소 : bos taurus
오리 : anatidae
박쥐 : chiroptera

image

Code

1. Numeric Representation of the Sequences

염기서열 스트링으로부터 feature를 생성하는 함수입니다.
염기서열은 알파벳의 나열로 나타내지기 때문에 일정길이로 잘라서 숫자로 변경한 뒤 분석을 진행할 것입니다.
이렇게 자르게 되면 단어간의 순서나 각 염기서열토큰의 출현빈도등의 정보를 유지할 수 있습니다.

CountVectorizer와 scikit-learn 라이브러리를 사용하여 단어들의 벡터를 생성해주겠습니다. 

#염기서열을 토큰별로 잘라서 수치화한 다음, 리스트로 반환
def generate_ngrams(s1):
    count_vect = CountVectorizer(lowercase=False, ngram_range=(2,4),analyzer='char')
    X1 = count_vect.fit_transform(s1)

    lcount = list()
    lcount = []
    for i in s1:
        count = len(i)
        #print(count)
        lcount.append(count)

    count_vect_df = pd.DataFrame(X1.todense(), columns=count_vect.get_feature_names())
    count_vect_df=count_vect_df.apply(lambda x: x / lcount[x.name] ,axis=1)

    return count_vect_df

2. Data Transformation

데이터 가공파트입니다. 얻은 5가지 데이터들을 fasta파일에서 pandas dataframe형식으로 변경해줘야 합니다.
'>'로부터 시작하는 description은 제거하고 ATTAAG같은 염기서열데이터만 뽑아주겠습니다.

염기 서열 데이터만 뽑으면 위의 generate_ngrams함수를 통해 염기서열 스트링을 토큰 리스트로 변경시킬 것입니다. 

#파일로부터 스트링읽어서 각염기서열의 description삭제
def process_file(filename,target_val):
    f = open(filename) #'datasets\\corona-nucleo-chicken-complete.fasta')
    lines = ""
    s1 = list()
    step = 0
    term = 0
    for line in f:
        line = ''.join(line.split())
        #print('step: ',step,' ',line)
        if line.startswith(">") and step==0:
            line = line.split('>',1)[0].strip()
            step = step + 1
        if line.startswith(">") and step>=1:
            line = line.split('>',1)[0].strip()
            s1.append(lines)
            lines = ""
            step = step + 1
            term = 0
        lines = lines + line

    count_vect_df = generate_ngrams(s1) 
    count_vect_df['target'] = target_val
    return count_vect_df

3. Exploratory Analysis

각 dataframe을 합쳐서 %로 출력해보면 다음과 같은 그래프를 얻을 수 있습니다.
현재 데이터양은 닭이 63%를 차지하고 박쥐가 32% 소 오리순으로 있습니다. 

# 각 데이터셋의 % 그래프로 출력
import matplotlib.pyplot as plt
plot_size = plt.rcParams["figure.figsize"]
plot_size[0] = 8
plot_size[1] = 6
plt.rcParams["figure.figsize"] = plot_size

df=pd.concat([df1,df2,df3,df4])

# NaN값을 지닌 열은 삭제
df=df.dropna(axis=1)
df['target'].value_counts().plot(kind='pie', autopct='%1.0f%%')

image

4. Data Preprocessing

동물 데이터들은 df라는 이름의 dataframe으로 합쳐놨습니다.
이제 사람 데이터를 로드하여 동물 데이터와 동일한 열 정보를 가질 수 있게 데이터 전처리를 진행해 줄 것입니다.

사람 데이터를 cov라는 이름으로 로드해주고 

cov = process_file('genome/homo_sapiens.fasta',"COVID-19")

#모델로 사용할 것이 아니기때문에 target은 drop
cov = cov.drop('target', axis=1)

두 데이터의 열을 보면
image

image

df는 348열, cov는 923개의 열을 가지고 있습니다.
두 데이터의 포맷이 동일하지 않으니 일정하게 맞춰주어야 합니다. 두 데이터가 동시에 가지고 있는 열만 남도록 나머지는 drop시켜주도록 하겠습니다.

동물에겐 있고 사람에겐 없는 열 삭제하기

# 동물에겐 있고 사람에겐 없는 열찾기
y=df.pop('target')
mc = df.columns.difference(cov.columns)
mc
#해당 열 삭제
df = df.drop(mc, axis=1)

y=df.pop('target')는 모델 트레이닝 할 때의 라벨을 위해 미리 추출해놓겠습니다.
사람 데이터에서 target부분을 삭제했기 때문에 아래 과정에서 동물데이터의 target도 삭제되기 때문입니다.

사람에겐 있고 동물에겐 없는 열 삭제하기

#사람에겐 있고 동물에겐 없는 열찾기
rf = cov.columns.difference(df.columns)
rf
#해당 열 삭제
cov = cov.drop(rf, axis=1)

5. Building the Predictive Model

이제 예측 모델을 만들 준비가 다끝났습니다. xgboost(Extream Gradient Boosting) 알고리즘으로 모델을 생성하겠습니다.

XGBoost?
머신러닝 앙상블 boosting 알고리즘.
병렬처리가 가능하기 때문에, 학습과 분류가 빠르다. 

#모델 생성
from sklearn.model_selection import train_test_split 
from xgboost import XGBClassifier
from xgboost import plot_importance
import xgboost

#y=df.pop('target')
X=df.values

# create a train/test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=7, shuffle=True)

model = XGBClassifier()
model.fit(X_train, y_train)

6. Prediction

XGBoost모델에 dataframe형식의 데이터를 input으로 넣으면 feature_names mismatch 에러가 발생합니다.

참고 : ValueError: feature_names mismatch: in xgboost in the predict() function

그래서 모델 생성할 때나 예측할 때 전부 데이터를 ndarray형식으로 변환해주어야 합니다. 

# dataframe -> numpy expression
c=cov.values

동물 데이터셋인 df는 위에 모델 생성하기 전에 ndarray형식으로 바꿔주었습니다.

그리고 생성했던 모델에 사람데이터를 input으로 넣으면 다음과 같이 결과가 출력됩니다. 

model.predict(c)
array(['bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat', 'bat',
       'bat'], dtype=object)

사람의 데이터 199개에 대한 예측 결과가 출력됩니다.
더 분석해보기 위해 predict_proba API를 사용해서 각 예측에 대한 확률정보를 얻어보겠습니다. 

import numpy as np
print(model.classes_)
similarities = model.predict_proba(c)
np.round(similarities, 3)
['bat' 'cattle' 'chicken' 'duck']
array([[0.987, 0.001, 0.011, 0.001],
       [0.987, 0.001, 0.011, 0.001],
       [0.987, 0.001, 0.011, 0.001],
       [0.987, 0.001, 0.011, 0.001],
       [0.987, 0.001, 0.011, 0.001],

...

       [0.985, 0.001, 0.012, 0.001],
       [0.987, 0.001, 0.011, 0.001],
       [0.987, 0.001, 0.011, 0.001],
       [0.987, 0.001, 0.011, 0.001],
       [0.987, 0.001, 0.011, 0.001],
       [0.987, 0.001, 0.011, 0.001]], dtype=float32)

예측 결과의 확률정보를 살펴보니 박쥐의 확률이 다른 세마리동물보다 월등히 높다는 것을 확인할 수 있습니다.

이 결과 대로라면 박쥐가 맨 처음 코로나바이러스의 진원지로 지목당한 것도 어느정도 이유가 있는 선택이었겠다고 생각이 듭니다.

Appendix -데이터의 NaN값 처리에 대하여

사실 데이터 전처리에 대한 부분은 원문과는 약간 내용과 순서가 다릅니다.
원문에서의 데이터 처리는

  1. 동물데이터 로드 -> 각 데이터 concat -> NaN값 그대로 놔두고 모델생성
  2. 사람데이터 로드 -> (동물o 사람x)인 데이터는 -999, (동물x 사람o)인 데이터는 drop시킴

위와같은 과정으로 진행을 하게 되어있습니다.

비교하기 위해 본문 과정도 나열하자면 다음과 같습니다.

  1. 동물데이터 로드 -> 각 데이터 concat -> NaN값이 1개라도 있는 열은 drop
  2. 사람데이터 로드 -> (동물o 사람x)인 데이터는 동물데이터에서 drop, (동물x 사람o)인 데이터는 사람데이터에서 drop시킴

두 과정의 차이는 크게 동물데이터의 NaN값을 그대로 두느냐로 갈립니다.

일단 동물데이터에 NaN값이 생기는 이유는, 4마리 동물 데이터를 합치는 과정에서 index가 일치하지 않기 때문입니다.
이러한 NaN값은 Pandas 내에서 계산을 할 때 그 자체적으로 수행이 안됩니다. 내부적으로 np.nan으로 처리가 되는데, 이 값에 어떤 값을 더하거나 빼도 NaN값이 반환되어 정상적인 결과를 얻을 수 없습니다. 그래서 반드시 NaN값을 어떤식으로도 처리를 해야합니다.

여러 방식으로 NaN값을 처리하는 방법이 인터넷이 많이 올라와 있습니다.
NaN을 배제하고 계산한다던지, NaN을 다른 수로 변환해서 계산한다던지 하는 방법이 올라와 있는데 제가 최종적으로 이번 포스팅에서 선택한 방법은 NaN을 완전히 배제하고 연산을 수행하는 것 입니다.

이번 챕터에서는 제가 수행했던 여러 작업들을 간략하게 설명해드리려고 합니다.

용어정리
df : 4마리 동물을 concat한 dataframe 이름
cov : 사람의 dataframe 이름

1번째 시도 : 원문 그대로(creating model -> preprocessing)

원문에서 한 방법은 위에서도 잠깐 언급했듯이, 데이터안의 NaN값을 그대로 놔두고 모델을 생성한 뒤, 각 데이터셋의 인덱스를 맞춰주는 순서로 진행됩니다.

df의 NaN값은 그대로 두고, cov의 새로생성된 index에 -999를 집어넣어서 predict를 했습니다.
결과는 예상했던 박쥐가 아니라 닭이 나왔죠.

image

물론 원문과 결과가 동일하지 않다고 해서 무조건 틀린 결과라고는 할 수 없습니다.

본문의 데이터셋은 대략 동물:580 사람:200 인 반면에 원문의 데이터셋은 사람:37 (동물의 데이터 수는 정확하게 나와있진 않습니다.) 정도로 적습니다.

동일한 데이터를 사용하지 않은 만큼 별개의 결과가 나올 수 있지만, 원문에서 predict에 대한 확률과 제 결과를 비교해보니 너무나도 차이가 나더군요.
image
그래서 뭔가 이상하다고 느꼈고, 무엇이 문제인지 찾아보기 시작했습니다.

2번째 시도 : df의 NaN값을 다른값으로 변경

데이터에 NaN값이 있으면 좋지 않다는 사실을 이 테스트를 진행하면서 처음 알게되었습니다.

그래서 df에 있던 모든 NaN값을 fillna함수를 통해 다른 값으로 치환해봤습니다. 

# df의 NaN값 -999로 치환
df.fillna(-999)
...
# cov에 새로 포함되는 값 -999로 생성
for newcol in mc: cov[newcol]=-999

번외로 0.0도 해봤습니다.

결론은 -999로 바꾸면 박쥐가 나오긴 하지만 닭과 확률이 비슷하여 돌릴때마다 결과가 변했고, 0.0으로 바꾸면 명확하게 닭으로 결과가 나왔습니다.

그나마 -999로 바꾼게 원문의 결과와 비슷하지만 여전히 확률쪽은 전혀 비슷하지 않습니다.

3번째 시도 : 전처리를 하고 모델만들기

이번에는 df와 cov의 index맞춰주는 전처리를 모델을 만들기 전에 해봤습니다.
여태까지는 df의 index는 건들지 않고 모델을 생성한 뒤 cov의 인덱스를 늘리고 줄이는 방법을 사용했었습니다.

df와 cov의 index맞춰주기 후, 모델생성

  • df의 NaN값 그대로 유지 : 박쥐 (소와 오리는 0)
  • df의 NaN값을 -999로 : 박쥐 (소와 오리는 0)
  • df의 NaN값을 0.0로 : 박쥐

총 3가지 결과를 얻었습니다.
우선 NaN값을 0.0으로 치환한 결과는 원문과 비슷한 확률을 가진 "박쥐"가 선택되었지만 NaN값을 그대로 유지하는 것과 -999로 변환하는 것의 결과는 박쥐이지만 소와 오리의 확률이 0이었습니다.

확률이 0이라니..아무리 낮아도 0은 아닌것 같아서 뭔가 이상함을 느꼈습니다.

4번째 시도 : NaN값을 아예 없애보자

AI쪽 하는 친구에게 조언을 구했더니, NaN을 아예 없애보면 어떠냐고 조언을 해줬습니다.
그래서 이번엔 df의 NaN값을 아예 없애버렸습니다.(dropna) 

# NaN값을 지닌 열은 삭제
df=df.dropna(axis=1)

이렇게 하니 모델을 만들고나서 index를 맞추든 index를 맞추고나서 모델을 생성하던 결과가 특정 확률이 0이되는 일 없이 일정하게 "박쥐"로 나왔습니다.

결론

아직 공부가 부족해서 NaN값을 가진 데이터를 무조건 날려버리는게 옳은 선택인지는 잘 모르겠습니다.
후에 좀 더 경험이 쌓이면 어떤 방향으로 처리를 하는것이 옳은지 알 수 있을겁니다.

그 때가 되면 또 관련 포스팅을 올리도록 하겠습니다.

반응형
저작자표시 비영리 변경금지