시뮬레이티드 어닐링(Simulated Annealing) 알고리즘으로 대학원생 친구 구원해주기

약 2달 전쯤 친구들이 모여있는 단톡방에서 포스텍 박사과정을 진행중인 친구가 여느때와 같이 우는소리를 하기 시작했습니다.
얘기를 들어보니 자기가 하는 연구가 대충 단백질 만드는거라는데 만들어진 단백질의 분자량만 가지고 그 단백질을 구성하는 아미노산들의 조합을 뽑아보고 싶다는겁니다. 
저런 화학지식에 대해 아는게 전혀없는 저에게는 뭔소린가 싶었지만 설명을 듣다보니 생각보다 할만하다? 라는 느낌이 들었습니다.
(물론 처음에 말한 조건과 나중에 말한 조건이 달라서 낭패를 보긴했습니다.. )
 
 

조건은?

처음에 저한테말해줬던 조건은 
1. 단백질 분자량 약 500~ 10000 정도의 값을 입력할 수 있어야한다. (double 값)
2. 아미노산의 종류와 각 분자량은 이미 20개 정도로 정해져있다.
3. 필수로 들어가는 아미노산을 입력할 수 있으면 좋다.
4. 아미노산의 조합을 목표값과 가까운 순서대로 20개 정도 뽑아줬으면 좋겠다.
 
이정도로 생각보다 간단해 보이는 주문이었습니다.

final aminoMap = {
  'G': 75.03,
  'A': 89.05,
  'S': 105.04,
  'T': 119.06,
  'C': 121.02,
  'V': 117.08,
  'L': 131.09,
  'I': 131.09,
  'M': 149.05,
  'P': 115.06,
  'F': 165.08,
  'Y': 181.07,
  'W': 204.09,
  'D': 133.04,
  'E': 147.05,
  'N': 132.05,
  'Q': 146.07,
  'H': 155.07,
  'K': 146.11,
  'R': 174.11,
};

 
위 20개가 계산에 쓰이는 아미노산의 고유 Seq 와 분자량을 의미한다. 즉, 입력값을 750 입력하면 GGGGGGGGGG (분자량 : 750.3) 이런식으로 출력이 되어야 한다는 얘기로 이해했습니다.
 
 

어떻게 구현할까?

개발을 해주기로 한 이상 어떤식으로 로직을 짤까 생각을 했는데, 이왕 하는거 저는 제 공부도 할겸, 그 연구실 학생들도 편하게쓸겸 MacOs, Windows 데스크탑 앱 배포를 결심했고 비지니스 로직을 수행해줄 주체를 고르기만하면 됐다.
1. dart 코드보다는 성능이 좋다는 python 코드로 짜고 Django로 서버를 띄워 요청을 보내고 response로 받는다.
2. flutter local 에서 파이썬 코드를 실행하고 결과를 받아오는 방법을 찾아본다.
3. dart code로 알고리즘을 구현한다.
 
솔직히 학창시절 알고리즘 공부는 java 아니면 python 으로만 짜봤고 flutter 개발자가 된 이후로 dart 언어를 통해 별다른 알고리즘을 해볼 기회는 없어서 고민을 하던중 아래와 같이 소거법으로 진행했습니다.
 
1번 방법은 서버를 띄워줘야 하는데 ec2 프리티어 서버는 만료된지 오래고 자원봉사하는데 서버까지 띄워주는건 투머치라고 판단.
2번 방법은 파이썬 코드를 실행하는 플러그인은 있지만 연산의 결과를 유연하게 받아오는건 안되는것으로 보임.
 
즉 그냥 생짜 dart 언어로 알고리즘을 짜보기로 했고 속도가 안나오면 그때 서버를 띄워도 늦지 않을거라 생각해 dart로 진행했습니다.
 
 

알고리즘을 골라보자!

로직의 주체는 골랐으니 이제 어떤 로직을 채택할지가 문제였는데 당연히 모든 경우의수를 계산한다는 단순한 해결책은 통할리가 없었습니다. 입력값이 10000쯤 되면 아미노산 시퀀스의 평균을 150정도로 잡아도  아미노산 조합의 길이가 60~70개는 될거고 그것은 2^60 ~ 2^70개 만큼의 경우의수가 있다는 얘긴데 심지어 가장 분자량이 작은 G의 경우에는 130개정도 되는 길이가 나오게 되고 이를 메모리 공간에 다 저장해두고 있다가 마지막에 한바퀴 쭉 돌면서 가장 근사치 20개 찾아줘! 할수는 없는노릇이었습니다.
 
효율적으로 답을 찾아올 방법을 고민했는데 듣도 보도 못한 알고리즘들을 비교하다가 결국 고른게 시뮬레이티드 어닐링(Simulated Annealing) 알고리즘 이었습니다. 이 개념이 꽤나 단순하면서 신기했습니다.
 
 

시뮬레이티드 어닐링..? 그게뭐야?

우선 어닐링의 의미부터 알아야 하는데 흔히 드라마나 영화 또는 만화에서 대장장이들이 칼을 만들때
1. 뜨거운 용광로에 철을 달구고
2. 망치로 두들기고
3. 찬물로 식히는
위 step을 계속해서 반복하면 원하는 칼의 형태가 만들어지는데 이 과정을 어닐링이라고 합니다.
 
그렇다면 시뮬레이티드 어닐링이란? 위 과정을 코드로 구현해 최적의 해(올바른 칼의 형태) 를 구해내는 과정인 것이지요
이 방법을 간략하게 설명하면 우선 초고온상태에 있다고 가정을 하고 얼마나 낮은 온도가 되어야 어닐링을 멈출지 미리 지정을 합니다.
 

const double initialTemperature = 10000.0; // 초기 온도
const double absoluteTemperature = 0.00001; // 최소 온도
const double coolingRate = 0.99; // 냉각률

 
우 세 값을 통해 초기 온도 부터 시작해 냉각률을 곱해가면서 최소온도가 될때까지가 한번의 싸이클입니다.

const int saIterations = 100; // 시뮬레이티드 어닐링 반복 횟수

 
하지만 여기서 우리는 최적의 해를 1개가 아닌 20개를 뽑아와야하기 때문에 중복을 고려해 넉넉하게 100번 시행하도록 합니다.
 

List<String> randomSolution(double targetMass) {
  List<String> solution = [];
  double mass = 0.0;
  while (mass < targetMass) {
    String aminoAcid = dataMap.keys.elementAt(random.nextInt(dataMap.length));
    double aminoAcidMass = dataMap[aminoAcid]!;
    if (mass + aminoAcidMass > targetMass) break;
    solution.add(aminoAcid);
    mass += aminoAcidMass;
  }
  return solution;
}

 
가장 처음에는 위 코드와 같이 목표가 되는 분자량과 멀어도 좋습니다. 대장장이가 창고에서 대충 철괴를 꺼내오듯 아무값이나 도출을 합니다.

  List<String> currentSolution = randomSolution(targetMass);

  List<String> bestSolution = List.from(currentSolution);

 
그리고는 대충만든값을 최적을 해 라고 판단하고 bestSolution 에 넣어놓습니다.

  while (temperature > absoluteTemperature) {
    List<String> newSolution = neighborSolution(currentSolution, targetMass);
    double newEnergy = evaluate(newSolution, targetMass);

    if (acceptanceProbability(currentEnergy, newEnergy, temperature) >
        random.nextDouble()) {
      currentSolution = newSolution;
      currentEnergy = newEnergy;
    }

    if (currentEnergy < bestEnergy) {
      bestSolution = List.from(currentSolution);
      bestEnergy = currentEnergy;
    }

    temperature *= coolingRate;
  }

 
그리고 나서는 10000.0(최고온도) 부터 시작된 두들김이 0.99를 반복해 곱해서 0.00001(최저온도) 가 될때까지 로직을 반복하며 망치를 두들깁니다.
현재 코드에서는 acceptanceProbability 함수가 1번의 망치 두들김이라고 보시면 되고 neighborSolution 함수는 기존 bestSolution 에서 하나의 아미노산 조합에서 랜덤으로 하나를 교체해 비슷하지만 다른 아미노산 조합을 만드는 함수입니다.

double acceptanceProbability(
    double currentEnergy, double newEnergy, double temperature) {
  return newEnergy < currentEnergy;
}

 
새로 만들어진 값과 목표 분자량과의 차이가 기존 bestSolution에 비교해서 적다면 새로운 값을 bestSolution에 넣어주는 것입니다.
여기서 이 방법의 장점이 나오는데 그것은 메모리에 하나의 값만 저장하고 있으면 된다는것입니다.
계산을 수도없이 하지만 메모리에 1개의 값만 가지고 비교해서 더 나은값을 찾으면 기존값은 몰라도 되고 새로운값만 알면 되기때문에 이점이 있다고 볼 수 있습니다.
 
이런식으로 어닐링을 100회정도 하고나면 그때그때 다른 100개의 최적의 해가 나오는데 여기서 중복을 제거하고 

  // 에너지(오차)에 따라 해들을 정렬하고 상위 20개를 선택
  bestSolutions.sort((a, b) => a.values.single.compareTo(b.values.single));
  bestSolutions = bestSolutions.take(topSolutionsCount).toList();

 
위 코드를 통해 20개를 잘라서 최종 결과를 도출해냅니다.
.
.
.
근데 여기서 잠깐.. 항상 그렇듯 기획자는 개발자에게 숨기는게 있기 마련이죠
대충 다 만들고 나니 한다는 말이 '어 근데 거기서 물 증발량 빼야해 그건 금방하지?' 라더군요 얘기를 들어보니 아미노산의 갯수 * 18.01 만큼은 최초 입력값에 포함 되어있지 않다는겁니다.
 
즉, 'G': 75.03 로 예시를 들면
입력값을 750.3 을 입력했으면 75.03개 10개인 GGGGGGGGGG 가 나오는게 아닌겁니다.
GGGGGGGGGG 의 분자량은 750.3 - (18.01 * 10) = 570.2 가 되어버리는겁니다.
여태 작업했던 코드 기준으로는 생뚱맞은 결과값밖에 얻을수가 없는거였습니다. 여기서 더 문제인건 아미노산조합의 길이는 아미노산마다 분자량이 다르기때문에 결과가 도출될때까지 알 수 가 없습니다. 
 
분자량 1000으로 따지면 G 기준으로는 10개가 넘지만 W 기준으로는 5개밖에 되지 않기때문에 들쭉날쭉하여 가장 분자량이 적은 아미노산, 가장 큰 아미노산으로 가능한 최소 길이와 최대 길이를 지정해서
만약 길이가 5~ 10이라면 물 증발량은 각각 90.05, 108.06, 126.07, 144.08, 162.09, 180.1 이 되고 
최초 입력값에서 위 6개 값을 순서대로 뺀 값을 기준으로 알고리즘을 실행해 결과를 구해야 했습니다.
 
사실 위 문제 말고도 이후에 다양한 추가 요구사항이 나왔지만 본문엔 안적는거로 하겠습니다...
 

결과는? 

아무튼 이런 우여곡절 끝에 포항공대 바이오 연구실에서 사용될 아미노산 시퀀스 유추 툴이 완성되었습니다.
현재도 업그레이드 중이고  친구 연구실 교수님이 제가 만든 프로그램을 보시더니 친구한테 저랑 같이 논문하나 내보는거 어떠냐고 요구사항을 몇가지 준비하시고 있계시다던데 벌써 두렵습니다.
 
만들어진 계산기는 현재 아래 링크에서 업데이트 되고있습니다. 
 
https://github.com/Gyeony95/Mass-Finder/releases

Releases · Gyeony95/Mass-Finder