2021년 12월 2일 목요일

감염병 진단용 PCR 프라이머 설계 그리고 결과물의 평가

PCR을 이용한 감염병 진단 제품을 개발하는 기업에게는 대단히 죄송한 이야기가 되겠지만 나는 지금까지 프라이머 설계를 학술적으로 매우 중요하고 매력적인 일이라고 생각하지 않았었다. 지난 여름이 되어서야 주변의 요청에 의하여 COVID-19의 유전체 서열을 바탕으로 PCR용 프라이머를 설계하는데 사용할 conserved sequence를 추출하는 작업을 하게 되었고, 그 과정에서 파생되는 일들을 정리하고 활용 범위를 확대하기 위한 가능성을 모색하는 중이다. 당연히 강한 흥미를 느끼면서 이 일에 조금씩 빠져들고 있다.

온라인에는 웹서비스로 돌아가는 다양한 PCR 프라이머 설계 및 평가용 툴이 존재한다. 어떤 프로그램은 소스 또는 실행파일을 가져다가 로컬 컴퓨터에서 직접 설치하여 사용하도록 배려해 주기도 한다. 예를 들자면 나는 EMBOSS의 primersearch를 이따금씩 사용한다. 이것은 프라이머 설계용 도구가 아니라 프라이머쌍의 서열(한 줄에 한 쌍)과 DNA 염기서열 정보를 입력하면 결합 위치와 amplimer(=amplicon)의 크기 정보를 Whitehead primer3_core 프로그램의 출력 파일 형태로 작성해 준다. 조정 가능한 파라미터 중에서 가장 중요한 것은 allowed percent mismatch이다. 그러나 서로 다른 쌍에 속하는 프라이머가 상호작용을 하여 다이머(dimer)를 만드는지의 여부를 점검해 주지는 않는다. 다시 말하자면 한 쌍으로 주어진 프라이머가 template DNA 서열에 mismatch를 고려하여 서로 마주보고 결합하는 것을 가정하여 생성되는 amplimer에 대한 결과만을 제시할 뿐이다. 프라이머 정보 파일을 구성하는 각 줄에 대해서 독립적인 분석을 수행하며, multiplexed PCR은 알 바가 아니다. 다음은 primersearch의 전형적인 실행 결과 파일이다.

Primer name primers-0
Amplimer 1
        Sequence: NC_045512
        Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome
        TGATGGTGGTGTCACTCGTG hits forward strand at 8701 with 0 mismatches
        GGCACGACAAAACCCACTTC hits reverse strand at [21039] with 0 mismatches
        Amplimer length: 165 bp

Primer name primers-1
Amplimer 1
        Sequence: NC_045512
        Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome
        GCCGCTGTTGATGCACTATG hits forward strand at 17170 with 0 mismatches
        TGTCGTCTCAGGCAATGCAT hits reverse strand at [12567] with 0 mismatches
        Amplimer length: 168 bp

프라이머끼리 서로 결합하여 어떤 다이머를 만드는지 점검하고 싶다면 ThermoFisher Scientific의 Multiple Primer Analyzer를 이용하면 좋다. 여기에서는 입력한 프라이머의 기본 특성 외에 self-dimer와 cross-dimer 정보를 제공한다.

이보다 더욱 진보한 프라이머 점검 도구에는 MFEprimer라는 것이 있다(논문 목록). 2009년에 최초 버전이 나왔으니 꽤 오랫동안 공들여서 기능을 추가하고 다듬어 온 도구라고 할 수 있다. 인간 등 주요 모델 생명체에 대한 'Background Database'를 선택하여 최대 50개까지의 프라이머에 대한 specificity check를 하게 배려해 준다. 명령행 버전을 설치하면 background database를 사용자가 공급할 수 있고, 입력하는 프라이머 수에 대한 제한도 없다. 내가 생각하는 MFEprimer의 가장 큰 특징은 당초 프라이머 설계 과정에서 누가 짝을 이루고 있는지를 전혀 고려하지 않고 amplicon을 예측해 준다. 물론 헤어핀과 다이머에 대한 점검은 기본이다. 최신 버전의 기능 소개 및 다른 버전과의 비교는 여기에 잘 설명되어 있다.

다음은 내 컴퓨터에 설치한 MFEprimer 명령행 버전에 100 쌍의 프라이머 서열을 넣고 분석한 첫 번째 amplimer의 구조이다.

Amp 1: primer-33-forward + primer-63-reverse ==> hCoV-19_USA_...(너무 길어서 생략):26844-26879

    Size = 35 bp, GC content = 54.29%
    F: Tm = 61.24 °C, Delta G = -23.36 kcal/mol
    R: Tm = 62.74 °C, Delta G = -23.58 kcal/mol
    Binding sites: 26845(20/20) ... 26879(20/20)

>>>primer-33-forward
   1                  20
5' GGACACcATCTAGGACGCTG 3'
   ::::::::::::::::::::              26879
5' GGACACCATCTAGGACGCTGTGACATCAAGGACCT 3'
   26845          ::::::::::::::::::::
               3' CGCTGTGAcATcAAGGAcCT 5'
                 20                  1
                     primer-63-reverse<<<

어? 이건 cross-dimer가 아닌가? 원래 다른 위치를 증폭하도록 설계된 프라이머가 5'-...CGCTG(끝)-3'라는 염기서열 매치를 통해서 다이머를 형성한 것으로 볼 수도 있다. ThermoFisher Scientific의 Multiple Primer Analyzer에서는 두 프라이머 사이에 단지 5개 염기가 연속적으로 결합하여 만들어지는 다이머는 표시하지 않는다. 아마도 두 프라이머 사이에 불연속적인 basepairing을 포함하여 최소한 6개 염기는 붙어야 다이머로 판단하는 것 같다. MFEprimer가 이를 amplimer로 판정하는 이유는, amplimer의 염기서열(5'-GGACACCATCTAGGACGCTGTGACATCAAGGACCT-3')이 template DNA 내 26845..26875 위치에 분명히 존재하기 때문일 것이다.

그러면 MFEprimer는 몇 % 혹은 몇 base의 mismatch까지 고려하여 amplimer를 예측하는가? 이렇게 갑자기 물으면 대답하기가 어려워진다. 원래 MFEprimer는 프라이머 서열 검색에 k-mer 알고리즘을 쓰기 때문이다. NGS short read의 조립이나 유전체 비교 계산에 널리 쓰이는 바로 그 k-mer 말이다. 그렇기 때문에 이전 버전에서는 k-mer 내부에서 일어나는 mismatch는 허용하지 못하였었다. 그러나 v3.0부터는 그 기능이 추가되었다고 한다. 다음의 사례를 보자. Reverse primer binding site 내에 1-bp mismatch가 보인다. 

>>>primer-61-forward
   1                  20
5' TACGCGTTCCATGTGGTCAT 3'
   ::::::::::::::::::::                                                    26937
5' TACGCGTTCCATGTGGTCATtcaatccagaaacta...tgctggacaccatctAGGACGCTGTGATATCAAGG 3'
   26765                                                ::::::::::::.:::::::
                                                     3' AGGACGCTGTGAcATcAAGG 5'
                                                       20                  1
                                                           primer-62-reverse<<<

3'-end의 mismatch는 PCR 반응에 바람직하지 않으므로 허용하지 않지만, 웹 버전에서는 이를 해제할 수 있다.

명령행 버전에서는 훨씬 상세하게 mismatch 관련 파라미터를 건드릴 수 있다. 자세한 사항은 논문의 'The k-mer mismatch search algorithm' 섹션을 상세히 읽어보기 바란다. 나도 아직 완벽하게 이해하지는 못했다.
  • --misMatch int   max allowed mismatches between kmer and its binding sites
  • --misStart int   mis-match starts from the position of 3' end, 1 is for the very end of primer. (default 1)
  • --misEnd int     mis-match ends to the position of 3' end, 9 (kvalue) for the 9th base from the 3'end (default 9)
Mismatch 관련 파라미터의 설정은 EMBOSS primersearch가 훨씬 더 이해하기 쉽고 직관적이다. 그러나 실행 속도는? MFEprimer가 말할 수 없이 더 빠르다.

primer3_core의 실행

꽤 오래 전에 primer3의 WWW 버전을 설치하여 조금 써 본 일이 있으나 명령행에서 primer3_core 실행파일은 이번에 처음 돌려 보았다. 프라이머를 고를 염기서열 정보를 인수(FASTA file)로 제공하는 것으로 생각했는데, 그게 아니라 입력 파일('input.txt') 내에 ID와 서열을 각각 한 줄로 써 넣어야 함을 알았다. 어떤 문서에서는 입력 파일 내에 FASTA file의 path를 넣으면 된다고 설명하였는데 아무리 애를 써 봐도 작동하지 않았다.

코로나바이러스의 유전체(29.9 kb)에 대하여 총 100개의 프라이머 쌍을 만들도록 input.txt를 짠 뒤 primer3_core를 실행해 보았다. 다음은 input.txt의 전문이다. 맨 마지막 줄의 '=' 문자는 생략해서는 안 된다.

PRIMER_TASK=generic
PRIMER_PICK_LEFT_PRIMER=1
PRIMER_PICK_INTERNAL_OLIGO=0
PRIMER_PICK_RIGHT_PRIMER=1
PRIMER_OPT_SIZE=20
PRIMER_MIN_SIZE=18
PRIMER_MAX_SIZE=22
PRIMER_PRODUCT_SIZE_RANGE=150-200
PRIMER_EXPLAIN_FLAG=1
PRIMER_LOWERCASE_MASKING=1
PRIMER_NUM_RETURN=100
PRIMER_THERMODYNAMIC_PARAMETERS_PATH=/etc/primer3_config/
SEQUENCE_ID=NC_045512
SEQUENCE_TEMPLATE=ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCA...(생략)
=

Template 서열 내에서 각 프라이머가 어떻게 위치하는지를 시각적으로 보여주는 방법이 있다면 참 좋겠다는 생각이 든다. 물론 염기서열 정보를 다루는 일반적인 상용 GUI 소프트웨어에는 프라이머를 설계하여 그 위치를 표시하는 기능은 다 들어 있을 것이다. 그러나 프라이머 정보(염기서열 및 위치)를 별도로 갖고 있을 때 이를 일괄적으로 임포트하는 기능은 아마 없을 것이다. 나의 프로그래밍 기법이 충분하다면 아래와 같은 수준의 이미지를 만들어 봄직도 한데 말이다. Perl의 Bio::Graphics를 이용하여 이런 출력물을 만드는 방법을 공부해 볼까? 아이고, 인공지능과 기계학습, 신경망, 딥 러닝과 같은 buzzword가 난무하는 2021년 마지막 달에 과연 합당한 결심일지 모르겠다. '메타버스'라는 말도 있구나! 정말로 관심이 있다면 Perl로 짜여진 PrimerView(2015)부터 공부를 하는 것이 좋을 것이다.

출처: PrimerDesign-M

Primer3_core가 설계한 프라이머가 만들어낼 100개의 amplimer는 서로 조금씩 위치가 겹칠까? 한쪽 프라이머를 고정해 놓고, 반대편 프라이머를 여러 개 선발하여 이를 별도의 amplimer로 취급할까? 이런 세부적인 사항을 건드리는 파라미터는 위 input.txt 파일에서 찾아볼 수 없다. primer3의 논문(NAR 2012)에서 새로 도입된 파라미터인 SEQUENCE_OVERLAP_JUNCTION_LIST는 내가 궁금해하는 사항과는 거리가 멀다. 단, left primer를 하나 지정해 두고 여러 right primer를 골라서 별도의 amplimer로 취급하는 경우는 확인을 하였다. left primer의 이름은 달리 부르고 있었으나 염기서열과 위치를 보니 전부 동일한 것이었다. 예를 들어 PRIMER_LEFT_61는 이름은 다르게 불리지만 총 7개의 amplimer에서 공통적으로 쓰이고 있었고, amplimer size는 159~180였다. 프라이머의 길이가 20 bp이므로 7개의 right primer 중 많은 것이 서로 겹치고 있음이 명백하다. 

Sequence Manipulation Suite: Primer Map에서 그린 소박한 그림. Primer3_core가 어떤 생각을 갖고 프라이머를 설계했는지 한눈에 알아볼 수 있다. 이렇게 만든 프라이머를 전부 때려넣고 MFEprimer를 돌리면 어찌 헷갈리지 않겠는가? 

만약 이 프라이머 서열을 MFEprimer에 넣어서 점검하려면, 동일한 프라이머는 제거하고 하나만 넣는 것이 옳을 수도 있다. 만약 모든 것이 성가시고 multiplexed PCR을 하지 않을 것이 확실하다면, 프라이머 쌍 하나에 대하여 개별적으로 MFEprimer를 실시하는 것도 바람직하다.

MFEprimer에 대해서 상세히 조사하게 된 것은 PCR Primer Design for the Rapidly Evolving SARS-CoV-2 Genome(Methods Mol Biol 2022? 내년에 정식 출판될 예정인가?)이라는 논문을 따라 하면서 많은 흥미를 느꼈기 때문이다. fIDBAC을 설치하여 사용하려다가 많은 실망을 했었는데, 이번에는 꽤 좋은 도구를 만난 것 같아서 반갑게 느껴진다. 튜토리얼 형식으로 되어 있어서 매우 쉽게 실행 가능하였는데, 변이 위치를 BED의 구간으로 전환하는 awk 공식은 좀 잘못되어 있는 듯하다. 이에 대해서는 별도의 위키문서 SARS-CoV-2 검출을 위한 PCR 프라이머 설계에 기록해 놓았다.

Multiple sequence alignment를 이용한 PCR용 프라이머 설계 방법에 대해서는 다음 기회에...

확인할 사항

Primer3_core의 coordinate system은 BED 파일처럼 0-based 같다는 생각이 든다. 확인해 볼 필요가 있다. 'primer3_core -format_output' 옵션을 이용하여 출력 파일을 읽기 쉬운 형태로 바꾸어서 실행해 보니 내 예상이 맞았다. 출력 파일 넷째 줄에 'Using 0-based sequence positions'라는 글귀가 선명하다.

댓글 없음: