파이썬으로 알고리즘 배우기

0. 시작하기전에

이승찬, 모두의 알고리즘 with 파이썬, 길벗, 2017 을 읽고 복습을 위해 정리한 것입니다. 자세한 내용은 을 읽어주세요.

1. 알고리즘 기초

알고리즘이란 간단히 말해 어떤 문제를 풀기 위한 절차나 방법 입니다. 간단한 문제를 풀어 보면서 알고리즘에 대해 살펴 보겠습니다.

문제.1 1부터 n까지의 합구하기

In [1]:
# 해법1
def sum_n(n):
    s = 0
    for i in range(1, n + 1):
        s += i
    return s


sum_n(1000)
Out[1]:
500500
In [2]:
# 해법 2
def sum_n2(n):
    return n * (n + 1) // 2


sum_n2(1000)
Out[2]:
500500

입력 크기와 계산 횟수

입력의 크기(n)에 따라서 해법에 차이가 생깁니다.

  1. 해법1의 경우는: 덧셈 n 번
  2. 해법2의 경우 덧셈, 곱셈, 나눗셈 총 3번 입력의 크기가 작을때는 큰 차이가 없지만, n의 크기가 커질수록 속도에 엄청난 차이가 납니다.

O 표기법

계산복잡도를 표헌하는 방법 중 가장 많이 사용되는 방법입니다. 빅 오표기법이라고도 부릅니다. 간단한 예를 들면, 입력 크기 n에 대해 계산을 n번 한다면, 계산 복잡도를 O(n)이라고 표현합니다.

  • O(n) : 필요한 계산 횟수가 입력크기 n과 비례할 때
  • O(1): 필요한 계산 횟수가 입력크기와 무관할 때
  • O(n^2) : n의 제곱에 비례하여 계산 시간이 증가함
  • O(2^n) : 2의 n 제곱에 비례하여 계산 시간이 증가함

Jupyter notebook의 %%time 기능을 사용해 위의 두가지 해법의 속도를 비교해보겠습니다.

In [3]:
%%time
sum_n(100000)
CPU times: user 12 ms, sys: 0 ns, total: 12 ms
Wall time: 12.3 ms
Out[3]:
5000050000
In [4]:
%%time
sum_n2(100000)
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 10.5 µs
Out[4]:
5000050000

위의 두가지 알고리즘은 같은 문제를 해결하지만 시간에서 약 열배의 차이를 보여 주고 있습니다. 더 많은 양의 입력이 들어오면 차이는 더욱더 커지게 될겁니다. 문제를 해결하는데 있어서 다양한 방법이 존재할 수 있지만 우리가 가능하면 빠른 알고리즘을 짜야할 이유는 이 때문입니다.

문제2. 최댓값 찾기

주어진 숫자 n개 중 가장 큰 숫자를 찾는 알고리즘을 만들어 보세요

In [5]:
def find_max(x):
    n = len(x)
    max_v = x[0]
    for i in range(1, n):
        if x[i] > max_v:
            max_v = x[i]
    return max_v


v = [1, 2, 3, 4, 5, 6]
find_max(v)
Out[5]:
6

최댓값의 위치를 구하는 알고리즘

In [6]:
def find_index(x):
    n = len(x)
    max_index = 0
    for i in range(1, n):
        if x[i] > x[max_index]:
            max_index = i
    return max_index + 1  # 0부터 시작되기때문에 1을 더해줌


v = [1, 2, 3, 4, 5, 6]
find_index(v)
Out[6]:
6

문제3. 같은 값 찾기

n개의 데이터 중에서 같은 이름을 찾아 집함으로 돌려주는 알고리즘

In [7]:
def find_same_data(a):
    n = len(a)
    result = set()  # 결과 저장용 빈 집합
    for i in range(0, n - 1):
        for j in range(i + 1, n):
            if a[i] == a[j]:
                result.add(a[i])
    return result


data = ["C", "Cm", "D7", "E", "Cm"]
find_same_data(data)
Out[7]:
{'Cm'}

2. 재귀호출

재귀호출은 함수가 자기 자신을 다시 호출 하는 것을 뜻합니다.

def hi():
    print('hi')
    hi() # hi함수를 다시 호출

이것이 바로 재귀 호출입니다. 'hi'라는 문장을 출력한 다음 다시 자기 자신인 hi()를 호출 하므로 영원히 반복하는 것입니다.

문제4. 팩토리얼 구하기

1부터 n까지 연속한 정수의 곱을 구하는 알고리즘을 만들어 보세요.

In [8]:
# 해법 1
def fact1(n):
    result = 1
    for i in range(1, n + 1):
        result *= i
    return result


fact1(10)
Out[8]:
3628800
In [9]:
# 해법2, 재귀 호출 사용
def fact2(n):
    if n <= 1:
        return 1
    return n * fact2(n - 1)


fact2(10)
Out[9]:
3628800

1부터 n까지 합을 구하는 문제1을 재귀호출로 다시 한번 풀어보겠습니다.

In [10]:
def sum_n(n):
    if n == 0:
        return 0
    return sum_n(n - 1) + n


sum_n(100)
Out[10]:
5050

문제5. 최대공약수 구하기

두 자연수 a와 b의 최대공약수를 구하는 알고리즘을 만들어 보세요.

최대공약수는 두 개 이상의 정수의 공통약수중에서 가장 큰 값을 의미 합니다.

In [11]:
# 해법1
def gcd(a, b):
    i = min(a, b)
    while True:
        if a % i == 0 and b % i == 0:
            return i
        i = i - 1


gcd(102, 414)
Out[11]:
6
In [12]:
# 해법2, 유클리드 알고리즘을 이용해 최대 공약수를 구하는 알고리즘
def gcd(a, b):
    if b == 0:
        return a
    return gcd(b, a % b)


gcd(102, 414)
Out[12]:
6

문제6. 하노이 탑 옮기기

원반이 n개인 하노이 탑을 옮기기 위한 이동순서를 출력해보세요.

In [13]:
def hanoi(n, from_pos, to_pos, aux_pos):
    if n == i:
        print(from_pos, " >> ", to_pos)
        return
    hanoi(n - 1, from_pos, aux_pos, to_pos)
    print(from_pos, " >> ", to_pos)
    hanoi(n - 1, aux_pos, to_pos, from_pos)

3. 탐색과 정렬

탐색은 여러 개의 자료 중에서 원하는 것을 찾아내는 것이고 정렬은 주어진 자료를 순서에 맞춰 나열하는 것입니다.

문제7. 순차탐색

리스트에 있는 첫 번째 자료부터 하나하나 비교하면서 탐색합니다.

In [14]:
# 해법 1
def search_list(a, x):
    n = len(a)
    for i in range(0, n):
        if x == a[i]:
            return i + 1
        else:
            pass
    return "None"  # 없으면 None을 출력합니다.


v = [14, 12, 512, 23, 15, 63, 84, 23, 6]
search_list(v, 23)
Out[14]:
4
In [15]:
# 해법2: 리스트에 같은 값이 여러개 있을 경우
# 모든 위치를 리스트로 돌려주는 알고리즘
def search_list(a, x):
    n = len(a)
    result = []
    for i in range(0, n):
        if x == a[i]:
            result.append(i + 1)
        else:
            pass
    return result


v = [14, 12, 512, 23, 15, 63, 84, 23, 6]
search_list(v, 23)
Out[15]:
[4, 8]

문제8. 선택 정렬

리스트 안의 자료를 한 번씩 비교하는 방법과 거의 같습니다. 비교횟수가 입력 크기의 제곱에 비례하는 O(n^2) 알고리즘이므로 시간이 굉장히 오래 걸립니다.

In [16]:
def sel_sort(a):
    n = len(a)
    for i in range(0, n - 1):
        min_index = i
        for j in range(i + 1, n):
            if a[j] < a[min_index]:
                min_index = j
                a[i], a[min_index] = a[min_index], a[i]
    print(a)


d = [2, 4, 5, 1, 7]
sel_sort(d)
[1, 2, 4, 5, 7]

문제9. 삽입 정렬

삽입 정렬은 특정 굥우를 제외하면 선택정렬과 같인 O(n^2)알고리즘 입니다. 그래서 시간이 굉장히 오래 걸립니다.

In [17]:
def ins_sort(a):
    n = len(a)
    for i in range(1, n):
        key = a[i]
        j = i - 1
        while j >= 0 and a[j] > key:
            a[j + 1] = a[j]
            j -= 1
        a[j + 1] = key


d = [2, 4, 5, 1, 7]
ins_sort(d)
print(d)
[1, 2, 4, 5, 7]

문제10. 병합 정렬

주어진 문제를 절반으로 나눈 다음 각각을 재귀 호출로 풀어가는 방식입니다. 이런 기법을 분할 정복 이라고 부릅니다. 계산의 복잡도는 O(n*logn)으로 선택정렬이나 삽입정렬보다 빠릅니다.

In [18]:
def merge_sort(a):
    n = len(a)
    if n <= 1:
        return
    mid = n // 2
    group1 = a[:mid]
    group2 = a[mid:]
    merge_sort(group1)
    merge_sort(group2)
    i1 = 0
    i2 = 0
    ia = 0
    while i1 < len(group1) and i2 < len(group2):
        if group1[i1] < group2[i2]:
            a[ia] = group1[i1]
            i1 += 1
            ia += 1
        else:
            a[ia] = group2[i2]
            i2 += 1
            ia += 1

    while i1 < len(group1):
        a[ia] = group1[i1]
        i1 += 1
        ia += 1
    while i2 < len(group2):
        a[ia] = group2[i2]
        i2 += 1
        ia += 1


d = [6, 3, 2, 6, 10, 23, 1]
merge_sort(d)
print(d)
[1, 2, 3, 6, 6, 10, 23]

문제11. 퀵 정렬

그룹을 둘로 나눠 재귀호출을하는 방식으로 앞서본 병랍 정렬과 유사하지만, 그룹을 나눌 떄 미리 기준과 비교해서 나누는 것이 다릅니다.

In [19]:
def quick_sort(a, start, end):
    if end - start <= 0:
        return
    pivot = a[end]
    i = start
    for j in range(start, end):
        if a[j] <= pivot:
            a[i], a[j] = a[j], a[i]
            i += 1
    a[i], a[end] = a[end], a[i]
    quick_sort(a, start, i - 1)
    quick_sort(a, i + 1, end)


def quick_sort_(a):
    quick_sort(a, 0, len(a) - 1)


d = [6, 8, 3, 9, 10, 1, 2, 4]
quick_sort_(d)
print(d)
[1, 2, 3, 4, 6, 8, 9, 10]

기준값의 중요성

퀵 정렬에서 좋은 기준을 정하는 것이 굉장히 중요합니다.

문제12. 이분 탐색

자료가 크기 순서대로 정렬된 리스트에서 특정한값이 있는지 찾아 위치를 돌려주는 알고리즘을 만들어보세요.

이분 담색은 '둘로 나눈다'는 뜻입니다.

탐색할 자료를 둘로 나누어 찾는 것이라 순차 탐색보다 훨씬 빨리 찾을 수 있습니다.

In [20]:
def binary_search(a, x):
    start = 0
    end = len(a) - 1
    while start <= end:
        mid = (start + end) // 2
        if x == a[mid]:
            return mid + 1  # 숫자1을 더해서 알기쉽게
        elif x > a[mid]:
            start = mid + 1
        else:
            end = mid - 1
    return None  # 찾지 못했을때


d = [6, 8, 3, 9, 10, 1, 2, 4]
print(binary_search(d, 8))
2

4. 자료구조

여러가지 자료와 정보를 컴퓨터 안에 저장하고 보관하는 방식

문제13. 회문 찾기

회문은 거꾸로 읽어도 내용이 같은 문장를 뜻합니다.

예를 들면 '기러기', '일요일' 등이 있죠.

In [21]:
def palindrome(s):
    que = []
    stack = []
    for x in s:
        if x.isalpha():
            que.append(x.lower())
            stack.append(x.lower())
    while que:
        if que.pop(0) != stack.pop():
            return False
    return True


print(palindrome("Kayak"))
True

문제14. 딕셔너리로 같은 이름 찾기

파이썬의 딕셔너리라는 자료 구조를 사용해 풀어보겠습니다.

In [22]:
def find_name(a):
    name_dic = {}
    for name in a:
        if name in name_dic:
            name_dic[name] += 1
        else:
            name_dic[name] = 1
    result = set()
    for name in name_dic:
        if name_dic[name] >= 2:
            result.add(name)
    return result


name = ["Jake", "Jim", "Jake"]
find_name(name)
Out[22]:
{'Jake'}

문제15. 친구의 친구찾기

친구관계를 이용해 직간접적으로 아는 모든 사람을 출력하는 알고리즘을 만들어 보세요.

In [23]:
def print_all_friends(g, start):
    que = []
    done = set()
    que.append(start)
    done.add(start)
    while que:
        p = que.pop(0)
        print(p)
        for x in g[p]:
            if x not in done:
                que.append(x)
                done.add(x)


friends = {
    "Summer": ["John", "Jake", "Mike"],
    "John": ["Summer", "Mike"],
    "June": ["Jake", "Tom"],
    "Jake": ["Tom"],
    "Tom": ["Summer"],
    "Mike": ["Summer"],
}
print_all_friends(friends, "Summer")
Summer
John
Jake
Mike
Tom

끝으로

알고리즘이란 무엇인가? 라는 질문에서 시작해 15가지의 문제를 통해 알아보았습니다. 풀어 본 문제들은 쉬운 문제지만 알고리즘에 대한 접근법은 동일합니다. 앞으로도 여러분들이 더 많은 문제와 해답을 위한 노력을 기울이길 기원합니다.

머신러닝으로 단백질 용해도 예측하기

0. 단백질 용해도 예측

유전자 재조합을 통해 특정 단백질을 과발현시키는 것은 생명과학에서 아주 중요한 일입니다. 그러나 대부분의 과발현 단백질은 inclusion bodies를 형성해 용해도가 매우 낮은 문제가 있습니다. 그렇기 때문에 단백질 용해도를 예측하는 일은 실험의 시행착오를 줄이는데 도움이 됩니다.

여기서는 대장균에서 발현되는 단백질 데이터를 가지고 기계학습으로 용해도를 예측하는 문제를 풀어보겠습니다.

1. eSOL 데이터셋 설명

분석에 사용한 eSOL(Solubility database of all E.coli proteins) 데이터셋은 Targeted Proteins Research Project을 통해 http://tp-esol.genes.nig.ac.jp/ 에 공개 되어있습니다.

공식 페이지에서는 해당 데이터셋에 대해 다음과 같이 설명 합니다.

eSOL is a database on the solubility of entire ensemble E.coli proteins (Niwa T, Ying BW, Saito K, Jin W, Takada S, Ueda T and Taguchi H, 2009) individually synthesized by PURE system that is chaperone free.

대장균(E.coli) 단백질의 용해도(Solubility)는 아래의 공식을 통해 계산되었습니다.

$$ Solubility = \frac{Supernatant fraction}{uncentrifuged fraction} * 100 $$

따라서 주의해야 할 것은 용해도가 동일해도 전체 발현량은 동일하지 않을 수 있습니다.

2. 데이터셋 불러오기

pandas 라이브러리를 사용해 데이터셋을 불러옵니다. 원본 데이터셋은 tsv형식으로 되어 있지만, 편의를 위해 csv형식으로 변경하고 단백질 서열 정보를 추가했습니다.

In [1]:
import pandas as pd

df = pd.read_csv("E_coli_protein.csv")
df.tail()
Out[1]:
Number Gene Solubility Sequence
3143 3148 ytfG 106 MIAITGATGQLGHYVIESLMKTVPASQIVAIVRNPAKAQALAAQGI...
3144 3149 ytfN 32 MSLWKKISLGVVIVILLLLGSVAFLVGTTSGLHLVFKAADRWVPGL...
3145 3150 yzcX 27 MNDSEFHRLADQLWLTIEERLDDWDGDSDIDCEINGGVLTITFENG...
3146 3151 yzfA 88 MRIFVYGSLRHKQGNSHWMTNAQLLGDFSIDNYQLYSLGHYPGAVP...
3147 3152 yzgL 30 MQNRKWILTSLVMTFFGIPILAQFLAVVIAMLGVGLAGIIEVCNIL...

불러온 데이터셋의 마지막 5개 데이터는 위와 같습니다. 데이터셋에는 총 3147개의 단백질의 유전자이름과 서열정보 그리고 용해도값이 들어 있습니다.

3. 추가 데이터 열 만들기

단백질 서열은 상직적(symbolic) 정보입니다. 따라서 를 포함하고 있습니다. 따라서 다양한 분석법을 통해 기계학습에 유용한 특성(feature)를 만들 수 있습니다.

먼저 전체 단백질 서열에서 각각의 아미노산의 비율에 대한 데이터를 추가해봅니다. 파이썬 기본 함수인 count()를 사용합니다.

3.1. 아미노산의 비율

In [2]:
AA = [
    "A",
    "C",
    "D",
    "E",
    "F",
    "G",
    "H",
    "I",
    "K",
    "L",
    "M",
    "N",
    "P",
    "Q",
    "R",
    "S",
    "T",
    "V",
    "W",
    "Y",
]

for i in AA:
    df[i] = [x.count(i) / len(x) for x in df["Sequence"]]

이제 Biopython 라이브러리에 있는 도구를 사용해 추가적인 데이터를 만들어 보겠습니다.

3.2. 단백질의 분자량

Biopython.ProteinAnalysis 객체를 사용합니다.

In [3]:
from Bio.SeqUtils.ProtParam import ProteinAnalysis

df["molecular_weight"] = [ProteinAnalysis(x).molecular_weight() for x in df["Sequence"]]

3.3. 다양한 계산값

단백질의 방향족(aromaticity) 아미노산 비율, 불안정성(instability index), 등전점(isoelectric point)등에 관한 값을 추가합니다.

In [4]:
df["aromaticity"] = [ProteinAnalysis(x).aromaticity() for x in df["Sequence"]]
df["instability_index"] = [
    ProteinAnalysis(x).instability_index() for x in df["Sequence"]
]
df["isoelectric_point"] = [
    ProteinAnalysis(x).isoelectric_point() for x in df["Sequence"]
]
df["gravy"] = [ProteinAnalysis(x).gravy() for x in df["Sequence"]]

3.4. 단백질 2차 구조

단백질 2차 구조의 비율에 대한 값을 추가합니다. 현재 단백질 서열을 통해 2차 구조를 예측하는 것은 상당히 정확한 편입니다.

In [5]:
df["structure_fraction_helix"] = [
    ProteinAnalysis(x).secondary_structure_fraction()[0] for x in df["Sequence"]
]
df["structure_fraction_turn"] = [
    ProteinAnalysis(x).secondary_structure_fraction()[1] for x in df["Sequence"]
]
df["structure_fraction_sheet"] = [
    ProteinAnalysis(x).secondary_structure_fraction()[2] for x in df["Sequence"]
]

3.5. 흡광 계수

단백질의 흡광 계수는 주로 농도를 계산하는데 사용되는 값으로 소수성 아미노산과 시스테인(cysteine)의 비율과 관련있습니다.

In [6]:
df["extinction_coefficient_reduced"] = [
    ProteinAnalysis(x).molar_extinction_coefficient()[0] for x in df["Sequence"]
]
df["extinction_coefficient_non_reduced"] = [
    ProteinAnalysis(x).molar_extinction_coefficient()[1] for x in df["Sequence"]
]

3.6. 단백질 유연성 정보

단백질 서열의 각 위치에 대한 아미노산 유연성을 계산합니다. 그런다음 평균값과 표준편차, 중간값을 추가합니다.

In [8]:
import numpy as np

df["flexibility_mean"] = [
    np.mean(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
]
df["flexibility_std"] = [
    np.std(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
]
df["flexibility_median"] = [
    np.median(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
]

각 아미노산의 유연성에 대한 조화평균, 기하평균, 최소값, 최대값, 왜도, 첨도값을 추가합니다.

In [10]:
from scipy.stats.mstats import hmean, gmean
from scipy.stats import skew, kurtosis

df["flexibility_hmean"] = [
    hmean(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
]
df["flexibility_gmean"] = [
    gmean(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
]
df["flexibility_min"] = [
    np.min(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
]
df["flexibility_max"] = [
    np.max(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
]
df["flexibility_skew"] = [
    skew(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
]
df["flexibility_kurtosis"] = [
    kurtosis(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
]

추가한 데이터를 확인해 봅니다.

In [12]:
df.tail()
Out[12]:
Number Gene Solubility Sequence A C D E F G ... gravy flexibility_mean flexibility_std flexibility_median flexibility_hmean flexibility_gmean flexibility_min flexibility_max flexibility_skew flexibility_kurtosis
3143 3148 ytfG 106 MIAITGATGQLGHYVIESLMKTVPASQIVAIVRNPAKAQALAAQGI... 0.171329 0.000000 0.052448 0.052448 0.017483 0.087413 ... 0.051748 1.000323 0.021540 0.999869 0.999860 1.000092 0.948167 1.054274 0.116313 -0.504896
3144 3149 ytfN 32 MSLWKKISLGVVIVILLLLGSVAFLVGTTSGLHLVFKAADRWVPGL... 0.066720 0.003177 0.065925 0.054805 0.023828 0.090548 ... -0.182685 1.003100 0.024786 1.002798 1.002490 1.002794 0.932905 1.079036 0.175824 -0.293648
3145 3150 yzcX 27 MNDSEFHRLADQLWLTIEERLDDWDGDSDIDCEINGGVLTITFENG... 0.047170 0.018868 0.113208 0.094340 0.047170 0.084906 ... -0.591509 1.001490 0.022863 1.000571 1.000967 1.001229 0.951655 1.048607 -0.084948 -0.807064
3146 3151 yzfA 88 MRIFVYGSLRHKQGNSHWMTNAQLLGDFSIDNYQLYSLGHYPGAVP... 0.061947 0.000000 0.070796 0.035398 0.017699 0.115044 ... -0.514159 0.995326 0.020807 0.995381 0.994891 0.995108 0.951667 1.043298 -0.000384 -0.537425
3147 3152 yzgL 30 MQNRKWILTSLVMTFFGIPILAQFLAVVIAMLGVGLAGIIEVCNIL... 0.075269 0.021505 0.010753 0.021505 0.086022 0.107527 ... 1.395699 0.969298 0.019198 0.963935 0.968922 0.969109 0.930988 1.022298 0.639452 -0.070879

5 rows × 43 columns

원래 4개의 행으로 구성되어 있던 데이터셋에 39개의 행이 추가되어 총 43개 행을 갖게 되었습니다. 수치로 되어있는 행에 대한 산술 통계값을 살펴봅니다.

In [11]:
df.describe()
Out[11]:
Number Solubility A C D E F G H I ... gravy flexibility_mean flexibility_std flexibility_median flexibility_hmean flexibility_gmean flexibility_min flexibility_max flexibility_skew flexibility_kurtosis
count 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 ... 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000 3148.000000
mean 1574.665502 49.002541 0.093592 0.012910 0.053903 0.061896 0.036498 0.070578 0.023977 0.058930 ... -0.147113 0.997623 0.024610 0.997044 0.997011 0.997317 0.940615 1.060624 0.105105 -0.541406
std 909.158084 33.398919 0.027916 0.011823 0.017879 0.022727 0.016239 0.023758 0.013145 0.020419 ... 0.337310 0.006706 0.002267 0.007453 0.006718 0.006712 0.008393 0.010076 0.174546 0.208449
min 1.000000 0.000000 0.012500 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 ... -1.712281 0.954846 0.015813 0.947970 0.954504 0.954674 0.912095 1.000571 -0.783561 -1.355014
25% 787.750000 18.000000 0.075758 0.005277 0.043011 0.047986 0.025890 0.055265 0.014509 0.044915 ... -0.337482 0.995471 0.023156 0.994565 0.994902 0.995187 0.935369 1.054357 0.010429 -0.664991
50% 1574.500000 43.000000 0.091882 0.011033 0.054054 0.061473 0.034617 0.069767 0.022788 0.057020 ... -0.182439 0.998133 0.024376 0.997607 0.997534 0.997832 0.941119 1.060905 0.104496 -0.552774
75% 2361.250000 79.000000 0.110102 0.017544 0.064695 0.074736 0.045304 0.085745 0.031821 0.070130 ... -0.020990 1.001112 0.025751 1.000841 1.000496 1.000800 0.945979 1.067405 0.198835 -0.435561
max 3152.000000 147.000000 0.308789 0.102857 0.166667 0.197531 0.149254 0.237288 0.132653 0.185185 ... 1.889655 1.035619 0.042443 1.036857 1.035343 1.035481 0.994155 1.100286 1.039824 0.970267

8 rows × 41 columns

위 테이블을 통해 용해도의 표준편차값이 높다는 것과 생성한 데이터들이 음수값도 포함하고 있다는 것을 알 수 있습니다.

3.7. 데이터 행 정렬

데이터 처리의 편의를 위해 Solubility, Sequence 행의 위치를 변경합니다.

In [13]:
df.columns.values
Out[13]:
array(['Number', 'Gene', 'Solubility', 'Sequence', 'A', 'C', 'D', 'E',
       'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T',
       'V', 'W', 'Y', 'molecular_weight', 'aromaticity',
       'instability_index', 'isoelectric_point',
       'structure_fraction_helix', 'structure_fraction_turn',
       'structure_fraction_sheet', 'extinction_coefficient_reduced',
       'extinction_coefficient_non_reduced', 'gravy', 'flexibility_mean',
       'flexibility_std', 'flexibility_median', 'flexibility_hmean',
       'flexibility_gmean', 'flexibility_min', 'flexibility_max',
       'flexibility_skew', 'flexibility_kurtosis'], dtype=object)
In [14]:
column_titles = [
    "Number",
    "Gene",
    "A",
    "C",
    "D",
    "E",
    "F",
    "G",
    "H",
    "I",
    "K",
    "L",
    "M",
    "N",
    "P",
    "Q",
    "R",
    "S",
    "T",
    "V",
    "W",
    "Y",
    "molecular_weight",
    "aromaticity",
    "instability_index",
    "isoelectric_point",
    "structure_fraction_helix",
    "structure_fraction_turn",
    "structure_fraction_sheet",
    "extinction_coefficient_reduced",
    "extinction_coefficient_non_reduced",
    "gravy",
    "flexibility_mean",
    "flexibility_std",
    "flexibility_median",
    "flexibility_hmean",
    "flexibility_gmean",
    "flexibility_min",
    "flexibility_max",
    "flexibility_skew",
    "flexibility_kurtosis",
    "Sequence",
    "Solubility",
]

df = df.reindex(columns=column_titles)
In [15]:
df.tail()
Out[15]:
Number Gene A C D E F G H I ... flexibility_std flexibility_median flexibility_hmean flexibility_gmean flexibility_min flexibility_max flexibility_skew flexibility_kurtosis Sequence Solubility
3143 3148 ytfG 0.171329 0.000000 0.052448 0.052448 0.017483 0.087413 0.027972 0.052448 ... 0.021540 0.999869 0.999860 1.000092 0.948167 1.054274 0.116313 -0.504896 MIAITGATGQLGHYVIESLMKTVPASQIVAIVRNPAKAQALAAQGI... 106
3144 3149 ytfN 0.066720 0.003177 0.065925 0.054805 0.023828 0.090548 0.007149 0.050040 ... 0.024786 1.002798 1.002490 1.002794 0.932905 1.079036 0.175824 -0.293648 MSLWKKISLGVVIVILLLLGSVAFLVGTTSGLHLVFKAADRWVPGL... 32
3145 3150 yzcX 0.047170 0.018868 0.113208 0.094340 0.047170 0.084906 0.028302 0.075472 ... 0.022863 1.000571 1.000967 1.001229 0.951655 1.048607 -0.084948 -0.807064 MNDSEFHRLADQLWLTIEERLDDWDGDSDIDCEINGGVLTITFENG... 27
3146 3151 yzfA 0.061947 0.000000 0.070796 0.035398 0.017699 0.115044 0.035398 0.044248 ... 0.020807 0.995381 0.994891 0.995108 0.951667 1.043298 -0.000384 -0.537425 MRIFVYGSLRHKQGNSHWMTNAQLLGDFSIDNYQLYSLGHYPGAVP... 88
3147 3152 yzgL 0.075269 0.021505 0.010753 0.021505 0.086022 0.107527 0.000000 0.118280 ... 0.019198 0.963935 0.968922 0.969109 0.930988 1.022298 0.639452 -0.070879 MQNRKWILTSLVMTFFGIPILAQFLAVVIAMLGVGLAGIIEVCNIL... 30

5 rows × 43 columns

3.8. 데이터셋 파일로 저장

아래 명령어를 통해 데이터셋을 파일로 저장 할 수 있습니다.

In [16]:
# df.to_csv("Engineered_E_coli_protein.csv", mode='w', index=False)

4. 시각화로 EDA 하기

간단한 시각화를 통해 데이터셋이 어떻게 분포되어 있는지 확인합니다.

4.1. Scatter plot

In [17]:
import matplotlib.pyplot as plt
import seaborn as sns

%matplotlib inline

sns.scatterplot(
    y="Solubility",
    x="molecular_weight",
    linewidth=0.1,
    alpha=0.5,
    sizes=(1, 8),
    data=df,
)
Out[17]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f3e74620048>
No description has been provided for this image

이미 알려져있는것 처럼 단백질의 크기가 클수록 용해도가 낮다는 것을 시각화를 통해 확인할 수 있습니다. 그러나 Scatter plot에는 겹쳐진 데이터가 너무 많기때문에 더 보기 좋은 그림을 그려봅니다.

4.2. Kernel density estimation plot

In [18]:
sns.jointplot(x="molecular_weight", y="Solubility", data=df, kind="kde")
Out[18]:
<seaborn.axisgrid.JointGrid at 0x7f3e7456be10>
No description has been provided for this image

위 시각화를 통해 대부분의 단백질의 크기는 30kDa 주변에 몰려있으며, 용해도는 20%와 80% 그룹으로 나뉜다는 것을 알 수 있습니다. 또한 상대적으로 높은 용해도의 단백질은 25kDa 이하의 크기입니다.

5. 데이터 전처리

기계학습을 위해 데이터셋을 특성(feature)값과 결과값으로 분리합니다.

In [19]:
df_y = df["Solubility"]

df_x = df.iloc[:, 2:-2]
df_x.columns
Out[19]:
Index(['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q',
       'R', 'S', 'T', 'V', 'W', 'Y', 'molecular_weight', 'aromaticity',
       'instability_index', 'isoelectric_point', 'structure_fraction_helix',
       'structure_fraction_turn', 'structure_fraction_sheet',
       'extinction_coefficient_reduced', 'extinction_coefficient_non_reduced',
       'gravy', 'flexibility_mean', 'flexibility_std', 'flexibility_median',
       'flexibility_hmean', 'flexibility_gmean', 'flexibility_min',
       'flexibility_max', 'flexibility_skew', 'flexibility_kurtosis'],
      dtype='object')

5.1. 특성값 정규화하기

정규화를 통해 기계 학습 모델이 특정 특성값에 바이어스(Bias)를 갖지 않게 합니다.

In [20]:
from sklearn import preprocessing
from sklearn.model_selection import train_test_split

scaler = preprocessing.StandardScaler().fit(df_x)
X_scaled = scaler.transform(df_x)
X_train, X_test, y_train, y_test = train_test_split(
    X_scaled, df_y, test_size=0.2, random_state=42
)

6. 기계학습하기

6.1. 다양한 기계 모델 비교

간단한 선형모델부터 다양한 기계 학습 모델을 사용해보고 평가 지표를 통해 가장 좋은 모델을 찾습니다. 예측의 정확도를 평가하기 위해 MAPE(Mean absolute percentage error)를 평가지표로 사용했습니다.

In [21]:
from sklearn.linear_model import LinearRegression
import numpy as np

lin_reg = LinearRegression(n_jobs=-1)  # -1 means use all cpu core
lin_reg.fit(X_train, y_train)
Out[21]:
LinearRegression(copy_X=True, fit_intercept=True, n_jobs=-1, normalize=False)
In [23]:
from sklearn.metrics import mean_absolute_error


def evaluate(model, X_test, y_test):
    y_pred = model.predict(X_test)
    mape = mean_absolute_error(y_test, y_pred)
    accuracy = 100 - mape
    return accuracy


accuracy = evaluate(lin_reg, X_test, y_test)
print(f"Accuracy = {accuracy:.2f} %")
Accuracy = 77.47 %

선형모델을 사용했을때 예측의 정확도는 약 77% 입니다.

In [24]:
from sklearn.tree import DecisionTreeRegressor

tree_reg = DecisionTreeRegressor()
tree_reg.fit(X_train, y_train)

accuracy = evaluate(tree_reg, X_test, y_test)
print(f"Accuracy = {accuracy:.2f} %")
Accuracy = 72.87 %
In [25]:
from sklearn.svm import SVR

svm_reg = SVR(kernel="linear")
svm_reg.fit(X_train, y_train)

accuracy = evaluate(svm_reg, X_test, y_test)
print(f"Accuracy = {accuracy:.2f} %")
Accuracy = 77.65 %
In [26]:
from sklearn.ensemble import RandomForestRegressor

forest_reg = RandomForestRegressor(n_estimators=1000, n_jobs=-1)
forest_reg.fit(X_train, y_train)

accuracy = evaluate(forest_reg, X_test, y_test)
print(f"Accuracy = {accuracy:.2f} %")
Accuracy = 79.42 %
In [27]:
from sklearn.ensemble import GradientBoostingRegressor

gbr = GradientBoostingRegressor()
gbr.fit(X_train, y_train)

accuracy = evaluate(gbr, X_test, y_test)
print(f"Accuracy = {accuracy:.2f} %")
# print(f'Score = {gbr.score(X_train, y_train)}') # R^2 socre
Accuracy = 79.48 %
In [28]:
# explicitly require this experimental feature
from sklearn.ensemble import HistGradientBoostingRegressor

hgr = HistGradientBoostingRegressor()
hgr.fit(X_train, y_train)

accuracy = evaluate(hgr, X_test, y_test)
print(f"Accuracy = {accuracy:.2f} %")
# print(f'Score = {hgr.score(X_train, y_train)}') # R^2 socre
Accuracy = 79.95 %

가장 최근에 추가된 기능인 HistGradientBoostingRegressor이 가장 좋은 성능을 보여줍니다.

6.2. 특성의 중요도

이제 데이터의 어느 특성이 모델 예측에 중요한지 살펴보겠습니다. 간단한 시각화로 막대그래프를 그려봅니다.

In [40]:
feature_importances = gbr.feature_importances_
plt.bar(range(len(feature_importances)), feature_importances)
Out[40]:
<BarContainer object of 39 artists>
No description has been provided for this image

위 그림을 통해 예측에 영향을 주는 특성을 알 수 있습니다. 상위 3가지를 알아보면 다음과 같습니다.

  1. 분자량(21번째 특성)
  2. pI값(29번째 특성)
  3. 분광계수(24번째 특성)

7. 모델의 매개 변수 최적화

GridSearchCV를 사용해 가장 좋은 성능을 나타내는 매개변수 값을 찾아봅니다.

In [35]:
from sklearn.model_selection import GridSearchCV
from sklearn.metrics import make_scorer


def evaluate(y_true, y_pred):
    mape = mean_absolute_error(y_true, y_pred)
    accuracy = 100 - mape
    return accuracy


score = make_scorer(evaluate, greater_is_better=True)

param_grid = [
    {
        "learning_rate": [0.01, 0.1],
        "max_depth": [10, 50, 100],
        "min_samples_leaf": [10, 50, 100],
    }
]

hgr_reg = HistGradientBoostingRegressor(random_state=42)
grid_search = GridSearchCV(hgr_reg, param_grid, cv=5, scoring=score, n_jobs=-1)
grid_search.fit(X_train, y_train)
Out[35]:
GridSearchCV(cv=5, error_score='raise-deprecating',
             estimator=HistGradientBoostingRegressor(l2_regularization=0.0,
                                                     learning_rate=0.1,
                                                     loss='least_squares',
                                                     max_bins=256,
                                                     max_depth=None,
                                                     max_iter=100,
                                                     max_leaf_nodes=31,
                                                     min_samples_leaf=20,
                                                     n_iter_no_change=None,
                                                     random_state=42,
                                                     scoring=None, tol=1e-07,
                                                     validation_fraction=0.1,
                                                     verbose=0),
             iid='warn', n_jobs=-1,
             param_grid=[{'learning_rate': [0.01, 0.1],
                          'max_depth': [10, 50, 100],
                          'min_samples_leaf': [10, 50, 100]}],
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring=make_scorer(evaluate), verbose=0)
In [36]:
cvres = grid_search.cv_results_
for mean_score, params in zip(cvres["mean_test_score"], cvres["params"]):
    print(f"Score = {mean_score:.2f}, Params = {params} ")
Score = 77.32, Params = {'learning_rate': 0.01, 'max_depth': 10, 'min_samples_leaf': 10} 
Score = 77.28, Params = {'learning_rate': 0.01, 'max_depth': 10, 'min_samples_leaf': 50} 
Score = 76.85, Params = {'learning_rate': 0.01, 'max_depth': 10, 'min_samples_leaf': 100} 
Score = 77.29, Params = {'learning_rate': 0.01, 'max_depth': 50, 'min_samples_leaf': 10} 
Score = 77.29, Params = {'learning_rate': 0.01, 'max_depth': 50, 'min_samples_leaf': 50} 
Score = 76.85, Params = {'learning_rate': 0.01, 'max_depth': 50, 'min_samples_leaf': 100} 
Score = 77.29, Params = {'learning_rate': 0.01, 'max_depth': 100, 'min_samples_leaf': 10} 
Score = 77.29, Params = {'learning_rate': 0.01, 'max_depth': 100, 'min_samples_leaf': 50} 
Score = 76.85, Params = {'learning_rate': 0.01, 'max_depth': 100, 'min_samples_leaf': 100} 
Score = 80.19, Params = {'learning_rate': 0.1, 'max_depth': 10, 'min_samples_leaf': 10} 
Score = 80.18, Params = {'learning_rate': 0.1, 'max_depth': 10, 'min_samples_leaf': 50} 
Score = 80.45, Params = {'learning_rate': 0.1, 'max_depth': 10, 'min_samples_leaf': 100} 
Score = 80.25, Params = {'learning_rate': 0.1, 'max_depth': 50, 'min_samples_leaf': 10} 
Score = 80.17, Params = {'learning_rate': 0.1, 'max_depth': 50, 'min_samples_leaf': 50} 
Score = 80.39, Params = {'learning_rate': 0.1, 'max_depth': 50, 'min_samples_leaf': 100} 
Score = 80.25, Params = {'learning_rate': 0.1, 'max_depth': 100, 'min_samples_leaf': 10} 
Score = 80.17, Params = {'learning_rate': 0.1, 'max_depth': 100, 'min_samples_leaf': 50} 
Score = 80.39, Params = {'learning_rate': 0.1, 'max_depth': 100, 'min_samples_leaf': 100} 

최고의 성능을 보인 모델의 정확도는 아래와 같습니다.

In [37]:
grid_search.best_score_
Out[37]:
80.4455794771026

매개변수 값이 아래와 같을때 최고의 성능을 보여주었습니다.

In [38]:
grid_search.best_params_
Out[38]:
{'learning_rate': 0.1, 'max_depth': 10, 'min_samples_leaf': 100}

8. 마치며

생성한 기계 학습 모델에 임의의 단백질 서열을 넣어 용해도를 예측해보겠습니다. 임의의 서열 2가지를 csv파일에 넣어 데이터 전처리과정을 동일하게 거칩니다.

In [42]:
df2 = pd.read_csv("sol_predict.csv")
df2
Out[42]:
Number Gene Sequence
0 1 random MAMIHHSAVQLSHGQTIMSGGQREGKSTLVQSIKVLTKQQQYPAPC...
1 2 DPO1_THEAQ MRGMLPLFEPKGRVLLVDGHHLAYRTFHALKGLTTSRGEPVQAVYG...
In [43]:
def feature_eng(df):
    for i in AA:
        df[i] = [x.count(i) / len(x) for x in df["Sequence"]]
    df["molecular_weight"] = [
        ProteinAnalysis(x).molecular_weight() for x in df["Sequence"]
    ]
    df["aromaticity"] = [ProteinAnalysis(x).aromaticity() for x in df["Sequence"]]
    df["instability_index"] = [
        ProteinAnalysis(x).instability_index() for x in df["Sequence"]
    ]
    df["isoelectric_point"] = [
        ProteinAnalysis(x).isoelectric_point() for x in df["Sequence"]
    ]
    df["gravy"] = [ProteinAnalysis(x).gravy() for x in df["Sequence"]]
    df["structure_fraction_helix"] = [
        ProteinAnalysis(x).secondary_structure_fraction()[0] for x in df["Sequence"]
    ]
    df["structure_fraction_turn"] = [
        ProteinAnalysis(x).secondary_structure_fraction()[1] for x in df["Sequence"]
    ]
    df["structure_fraction_sheet"] = [
        ProteinAnalysis(x).secondary_structure_fraction()[2] for x in df["Sequence"]
    ]
    df["extinction_coefficient_reduced"] = [
        ProteinAnalysis(x).molar_extinction_coefficient()[0] for x in df["Sequence"]
    ]
    df["extinction_coefficient_non_reduced"] = [
        ProteinAnalysis(x).molar_extinction_coefficient()[1] for x in df["Sequence"]
    ]
    df["flexibility_mean"] = [
        np.mean(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
    ]
    df["flexibility_std"] = [
        np.std(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
    ]
    df["flexibility_median"] = [
        np.median(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
    ]
    df["flexibility_hmean"] = [
        hmean(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
    ]
    df["flexibility_gmean"] = [
        gmean(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
    ]
    df["flexibility_min"] = [
        np.min(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
    ]
    df["flexibility_max"] = [
        np.max(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
    ]
    df["flexibility_skew"] = [
        skew(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
    ]
    df["flexibility_kurtosis"] = [
        kurtosis(ProteinAnalysis(x).flexibility()) for x in df["Sequence"]
    ]
    return df


df2 = feature_eng(df2)
df2
Out[43]:
Number Gene Sequence A C D E F G H ... extinction_coefficient_non_reduced flexibility_mean flexibility_std flexibility_median flexibility_hmean flexibility_gmean flexibility_min flexibility_max flexibility_skew flexibility_kurtosis
0 1 random MAMIHHSAVQLSHGQTIMSGGQREGKSTLVQSIKVLTKQQQYPAPC... 0.025000 0.015 0.030000 0.030000 0.035000 0.070000 0.055000 ... 64065 0.995175 0.025666 0.995929 0.994514 0.994844 0.94081 1.061024 0.103835 -0.531963
1 2 DPO1_THEAQ MRGMLPLFEPKGRVLLVDGHHLAYRTFHALKGLTTSRGEPVQAVYG... 0.109375 0.000 0.050481 0.104567 0.032452 0.069712 0.021635 ... 112760 1.000747 0.024468 0.999750 1.000151 1.000448 0.93856 1.077405 0.169728 -0.601805

2 rows × 42 columns

In [49]:
column_titles = [
    "Number",
    "Gene",
    "A",
    "C",
    "D",
    "E",
    "F",
    "G",
    "H",
    "I",
    "K",
    "L",
    "M",
    "N",
    "P",
    "Q",
    "R",
    "S",
    "T",
    "V",
    "W",
    "Y",
    "molecular_weight",
    "aromaticity",
    "instability_index",
    "isoelectric_point",
    "structure_fraction_helix",
    "structure_fraction_turn",
    "structure_fraction_sheet",
    "extinction_coefficient_reduced",
    "extinction_coefficient_non_reduced",
    "gravy",
    "flexibility_mean",
    "flexibility_std",
    "flexibility_median",
    "flexibility_hmean",
    "flexibility_gmean",
    "flexibility_min",
    "flexibility_max",
    "flexibility_skew",
    "flexibility_kurtosis",
    "Sequence",
]

df2 = df2.reindex(columns=column_titles)
df_x2 = df2.iloc[:, 2:-1]
X2_scaled = scaler.transform(df_x2)

이제 기계학습 모델을 사용해 용해도값을 예측해봅니다.

In [61]:
pred_sol = grid_search.best_estimator_.predict(X2_scaled)
gene = ["random", "DPO1_THEAQ"]
for idx, val in enumerate(gene):
    print(f"Predicted solubility of {val} is {pred_sol[idx]:.2f}%")
Predicted solubility of random is 28.88%
Predicted solubility of DPO1_THEAQ is 46.57%

위 코드를 통해 각각의 서열에 대한 용해도값을 예측 할 수 있었지만 실제값과 동일한지는 실제 실험을 해봐야할 것입니다.