Perl을 이용하여 multiple alignment에서 컬럼 단위의 conservation 추출하기

다중서열정렬(multiple-sequence alignment)에서 conserved region을 뚝딱 뽑아주는 프로그램이 널려 있을 것 같지만 현실은 그렇지 않다. trimAl이나 ClipKIT에서 그런 기능을 제공하지는 않는다. 원래 이들 프로그램의 용도는 다른 곳에 있기 때문이다. ClipKIT의 log 파일에서는 각 컬럼의 위치에 대하여 constant site 여부를 보여주기는 하지만(옵션에서 constant site를 추출하도록 지정하는 경우에 한함), 이를 수치로 출력한 뒤 사용자가 고르게 해 주지는 않는다. 그도 그럴 것이, 모든 서열에 대하여 동일한 잔기가 존재하는 경우 계통발생학적 분석에서 많은 정보를 주지 못하기 때문이다.

Biostars를 검색하니 이런 목적으로 Perl 코드를 이용하는 방법이 소개되어 있었다. 무려 10.5년 전의 질문과 답이었다!

Per-Column Conservation of Multiple Alignment in Perl

여기에 대하여 답을 제공한 Neilfws라는 사용자의 프로필을 따라가 보니 그는 내가 관심을 갖고 체크해 두었던 논문 "Pathogenic archaea: do they exist?" BioEssays 25:1119의 공저자인 Neil Saunders라는 사람이었다(구글 스칼라). 시드니에 거주하며 Curtin University에 근무하는 것으로 나온다.

파이썬의 세계에서는 이러한 요구사항을 충족하는 매우 좋은 방법이 있는지는 모르겠다. Neil의 답글에 따르면, Perl에서는 Bio::SimpleAlign 오브젝트를 이용하여 차근차근 풀어 나가는 방법이 전부이다. 

출처: Biostar

이 해법에서는 alignment를 컬럼 단위로 slice하여 처리하는 것이 핵심으로 보인다. BioPerl로 늘 FASTA 파일만 다루다가 이 스크립트를 따라하면서 AlignIO 및 SimpleAlign을 처음으로 다루어 본다. 아, 부끄러워라.

3만개가 넘는 SARS-CoV-2 genome의 MSA를 처리하는 것이 목표라서 일단은 100개 서열만 취하여 테스트를 해 보았다. 현재 사용하는 리눅스 서버에서 약 13분이 걸렸다. 이를 32,483으로 확장했을 때 단순히 서열의 수에 비례하여 시간이 증가한다고 가정하면... 어휴, 71시간이 걸린다. CPU를 하나밖에 사용하지 못하는 스크립트라는 것이 문제이다. 이를 병렬화하려면 어떻게 해야 하나? MSA를 위치 기반으로 썰어서 20개 정도로 나누어서 각각을 별도의 job으로 돌리고, 다시 이를 합치면서 coordinate를 재계산하고... 아, 귀찮아.

바로 곁에 있는 '최신' 컴퓨터엔 AMD Ryzen 9 5950x에서 시험 삼아 100개의 서열을 처리해 보았다. 5분 정도가 걸렸으므로 구형 Xeon 서버보다는 2.7배 정도 빠르다. 그래, 여기에서 돌리자. 어제 늦게 퇴근하면서 Ryzen 서버에서 작업을 시작하였는데 오전 10시 23분 현재 17500 site에 대한 계산이 끝났다. 오늘 퇴근 무렵이면 전체 29 kb에 대한 conservation 수치를 얻을 수 있을 것이다. 이를 사용하면 내 입맛에 맞게 기준을 정하여 conserved sequence를 뽑는 것이 가능하다. 100% conserved sequence를 추려내면 길이가 짧은 서열만이 나와서 후속 분석 작업에 쓰기 곤란하다는 것은 이미 ClipKIT 결과를 통해서 알고 있다.

trimAl은 정말 느리다. 그러나 alignment의 포맷 전환이라든가 위치 기반 추출 등 유용한 부속 기능이 있어서 가끔씩 쓰기에는 좋다.

독후감 쓸 거리가 몇 개 있는데 주말에 차분히 앉아서 블로그 작업을 할 수 있을지 모르겠다. 요즘은 뇌세포가 퇴화했는지 생각을 글로 풀어나가는 것이 예전보다는 어렵다.