====== SARS-CoV-2 검출을 위한 PCR 프라이머 설계 ======
  * Unordered List Item<color #22b14c>아직 위키문서 작성이 끝나지 않았습니다!  --- //[[hyjeong@kribb.re.kr|Haeyoung Jeong]] 2021/12/02 08:50//</color> 
  * 제 블로그에 별도의 글인 [[https://blog.genoglobe.com/2021/12/pcr.html|감염병 진단용 PCR 프라이머 설계 그리고 결과물의 평가]]를 게시하였습니다.
  * Genes & Genomics에 제 논문이 게재되었습니다. [[https://link.springer.com/article/10.1007/s13258-022-01264-7|Identification of conserved regions from 230,163 SARS-CoV-2 genomes and their use in diagnostic PCR primer design]] 본 위키페이지에서 다루는 내용과 100% 일치하지는 않습니다.  --- //[[hyjeong@kribb.re.kr|Haeyoung Jeong]] 2022/06/16 08:09//
===== 들어가는 글 =====
COVID-19 진단을 위해 쓰이는 PCR 기법은 빠르게 발생하는 변이체에 의해 민감도가 떨어질 우려가 크다. 새롭게 보고된 바이러스 변이체의 게놈 정보를 이용하여 돌연변이가 발생한 위치를 피하여 PCR 프라이머를 설계하는 방법이 다음과 같이 소개되어서 이를 재현해 보는 것이 이 문서의 작성 취지이다.
  * PCR Primer Design for the Rapidly Evolving SARS-CoV-2 Genome. Methods Mol Biol. 2022;2392:185-197. doi: 10.1007/978-1-0716-1799-1_14. https://pubmed.ncbi.nlm.nih.gov/34773624/
이 논문 말고도 엄청난 수의 SARS-CoV-2 검출용 PCR 기법이 넘쳐나고 있을 것이다. 이 논문을 택한 것은 이해하기 쉬운 기본 원리에 바탕을 두고 있고, 리눅스 컴퓨터에서 실행할 수 있는 튜토리얼 형태로 작성되어 있어서 따라하기에 좋기 때문이다. 

여러 DNA target을 공통적으로 증폭할 수 있는 프라이머를 설계하기 위하여 conserved region을 찾는 가장 보편적인 방법은 multiple sequence alignment(MSA)를 사용하는 것이다. 예를 들어서 EasyPrimer([[https://pubmed.ncbi.nlm.nih.gov/31992749/|논문]] [[https://skynet.unimi.it/index.php/tools/easyprimer/|웹사이트]])는 hypervariable gene을 공통적으로 증폭할 프라이머를 설계하기에 적당한 위치를 찾아주는 역할을 하며, PrimerDesign-M([[https://pubmed.ncbi.nlm.nih.gov/25524896/|논문]] [[https://www.hiv.lanl.gov/content/sequence/PRIMER_DESIGN/primer_design.html|웹사이트]])은 MSA를 입력하여 프라이머를 만들 수 있다. 

그러나 이 글에서 소개하는 Methods Mol Biol.의 방법은 시간이 많이 걸리는 MSA에 의존하지 않고 각 변이체 바이러스 게놈을 레퍼런스에 비교하여 변이의 위치를 VCF로 추출한 뒤 이를 병합한 다음, BED 파일로 전환하여 프라이머 설계용 reference 서열 위에 대/소문자로 변이 발생 위치를 표시하게 만든다. 이어서 Primer3의 파라미터를 조정하여 소문자(변이)가 3'-end에 오는 프라이머가 아닌 것을 고르게 한 뒤 MFEprimer로 최종 점검을 하는 방식을 취한다.
===== 분석 과정 요약 =====
  - [[https://pubmed.ncbi.nlm.nih.gov/32670259/|MicroGMT]]를 이용하여 reference genome에 대한 point mutation과 indel을 검출' Raw sequencing read와 assembled genome을 전부 사용할 수 있다. SARS-CoV-2에 최적화되어 있으나 다른 미생물 유전체에도 적용 가능하다고 한다. 따라서 snippy와 비교해 보는 것도 흥미로울 것이다.
  - '유명한' Primer3를 사용하여 프라이머를 설계
  - [[https://mfeprimer3.igenetech.com/spec|MFEprimer]]([[https://pubmed.ncbi.nlm.nih.gov/31066442/|논문]]) 명령행 버전을 사용하여 프라이머 평가

===== 실제 방법 =====
논문에서 소개한 명령어를 그대로 복사하여 실행하면 오류가 발생하기도 한다. 이를 바로잡아 나가도록 하자. 나는 <color #22b14c>sars_cov_2라는 conda environment(python 3.10.0)</color>를 마련하였는데, 이 분석 과정에서 필요로 하는 대부분의 프로그램이 시스템 수준에서 설치되어 있어서 큰 의미는 없었다. 전 과정은 /data/apps/MicroGMT(=${MICROGMT}) 디렉토리에서 실시하였다. 작업의 병렬화를 위하여 GNU parallel을 종종 사용하여 왔는데, 여기에서는 [[https://github.com/shenwei356/rush|rush]]라는 것을 사용한다. 
==== 1. Identify the mutations compared to the reference gnome ====
SARS-CoV-2 유전체 염기서열(FASTA file)을 다운로드하여 ${MICROGMT}/data에 보관한다. 나는 GISAID에서 다운로드한 1000개의 유전체 서열을 사용하였다(파일명: 'test_1000.fa'). 서열 ID에 슬래쉬와 파이프 기호('|')가 있어서 이를 각각 '_'와 '%%__%%'로 치환해 놓았다. 레퍼런스(NC_045512.2), hg19 인체 유전체 서열 및 인플루엔자 데이터베이스 설치는 본문을 참고하여 진행하여라.

  >hCoV-19/USA/ME-HETL-J0115/2020|EPI_ISL_755317|2020-09-02   #원본
  >hCoV-19_USA_ME-HETL-J0115_2020__EPI_ISL_755317__2020-09-02 # 수정 후
  
첫 명령어는 다음과 같다. variants 디렉토리를 미리 만들어 두어야 한다. 여러 쓰레드를 쓰고 싶다면 -t THREAD 옵션을 사용하라.
  $ python sequence_to_vcf.py -r data/NC_045512.fa -i assembly -fs data/test_1000.fasta -o variants
variants 디렉토리에 생성된 모든 VCF 파일([[https://samtools.github.io/hts-specs/VCFv4.3.pdf|VCF v4.3 specification 문서]])을 한 곳에 모은다. 원문에서는 '{}<color #ed1c24>s</color>'라고 표현되어 있는데, 's'가 있을 이유가 없다. rush가 여기에서 비로소 처음 쓰인다.
  $ find variants -name "*.vcf" | rush 'grep -v "#" {} >> all.variants.vcf'
다음과 같은 에러 메시지는 해당 게놈에 변이가 없음을 의미하므로 무시하면 된다.
  10:46:08.343 [ERRO] wait cmd #905: grep -v "#" variants/hCoV-19_Indonesia_YO-EIJK08_2020__EPI_ISL_766032__2020-03-16.vcf >> all.variants.vcf: exit status 1
다음 단계에서는 VCF 파일(1-based)을 BED 포맷(0-based)으로 전환한다. 본문에서 소개한 명령어의 첫 줄("awk '!/\#/' all.variants.vcf")는 정확하지도 않고 필요하지도 않다. 아마도 #로 시작하는 라인을 배제하기 위함으로 보이는데, 이미 이전 단계 명령어에서 이러한 라인은 전부 제거되었다. 또한 그러한 목적이라면 "awk '!/^#/' all.variants.vcf"라고 입력하는 것이 가장 적합할 것이다. 읽기 쉽도록 적당히 줄바꿈을 실시하였다.
  $ awk '{\
         if(length($4) > length($5)) \
                print $1"\t"($2-1)"\t"($2+length($4)-1); \
         else \
    	    print $1"\t"($2-1)"\t"($2+length($5)-1)\
         }' all.variants.vcf > all.variants.bed
과연 이 awk 명령어가 제대로 작동하고 있는 것인가? 다음은 204번째 위치에서 G가 T로 치환된 것을 표시한 것이다. 샘플로 사용한 100개 게놈 서열에서는 모두 158회 출현한다. 0-based BED에서는 2,3번째 컬럼에서 각각 203(chromStart), 203(chromEnd)으로 표현해야 할 것 같다. 그런데 왜 chromEnd는 204인가? 반드시 chromStart < chromEnd여야 할 필요는 없을 것이다. 왜냐하면 길이가 1 bp에 불과한 feature도 얼마든지 있을 것이기 때문이다. 이것이 대단히 석연치 않다. 차라리 [[https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/vcf2bed.html|vcf2bed]]와 같은 프로그램을 설치하여 테스트하는 것이 더 나을 수도 있다. <color #ed1c24>이것은 BED의 형식을 제대로 알지 못한 것에서 기인한 오해였다. 상세한 설명은 블로그 기록 [[https://blog.genoglobe.com/2021/12/vcf-bed.html|VFC 파일을 BED로 전환할 때 좌표 정보는 어떻게 바꾸어야 하는가]]를 참조하기 바란다.</color>
  [VCF 파일] NC_045512	204	.	G	T
  [BED 파일] NC_045512	203	204
  
다음 단계는 BED 파일에서 정의한 구간을 정렬하여 서로 겹치는 것을 병합한다.
  $ cat all.variants.bed | bedtools sort | bedtools merge > all.variants.merged.bed
==== 2. Design primers ====
먼저 ${MicroGMT}/design 디렉토리부터 만든다. BED 파일에 정의된 변이 발생 위치에 해당하는 염기를 소문자로 전환하는 "soft" masking을 실시한다. 보통의 마스킹에서는 N으로 치환해 버리므로 염기 정보를 잃어버리게 된다.
  $ bedtools maskfasta -soft -fi data/NC_045512.fa -bed all.variants.merged.bed -fo data/NC_045512.softmask.fa
  $ head -n 5 data/NC_045512.softmask.fa 
  >NC_045512
  ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTcTAAA
  CGAACTTTAAAATCTGTGTGGCTGTCACTcGGCTGcATGCTTAGTGCAcTCACGCAGTATAATTAATAAC
  TAATTACTGTCGTTGACAGGAcACGAGtAACTCgTCTATcTTCTGcAGgCTGCTTacggtttcgtccGtg
  TTGCAGCCGATcATCAGCAcATcTAGGTTTcGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGcCTTGTC

이 과정에서는 primer3_core 명령을 실행하는 것이 핵심이다. primer3가 유명한 프로그램이지만 명령행 환경에서 직접 실행할 일은 별로 없다. 따라서 사용법을 이번 기회에 알아 두는 것도 유익할 것이다. primer3는 프라이머 설계를 위한 파라미터와 염기서열을 하나의 텍스트 파일(input.txt)에 수록하여 STDIN으로 제공해야 한다. 즉, 염기서열 FASTA file을 primer3의 인수로 주는 것이 아니라는 뜻이다. primer3용 input.txt 파일에는 염기서열을 'SEQUENCE_TEMPLATE=' 뒤에 한 줄로 삽입하거나 혹은 한 줄로 만들어진 파일의 전체 경로를 넣어 주어야 한다. 따라서 변이 발생 위치를 soft masking한 FASTA file을 다음 명령을 이용하여 한 줄로 전환한다.
  $ awk '/^>/ {printf("% s\n",$0);next;} {printf("%s",$0);} END {printf("\n");}' \
            < data/NC_045512.softmask.fa \
	    > data/NC_045512.softmask.one-line.fa.txt

primer3_core 실행에 사용할 input.txt 파일의 샘플은 /usr/share/doc/primer3/examples/example이다(apt로 primer3를 설치한 경우). 이를 가져다가 수정과 추가를 거쳐서 ${MicroGMT}/design/input.txt로 저장한다. 그런데 정말 납득이 가지 않는 점이 있었다. 본문에서는 'SEQUENCE_TEMPLATE=../data/NC_045512.softmask.one-line.fa.txt' 형식으로 염기서열 파일을 지정하라고 하였지만 이렇게 해서는 primer3_core가 제대로 돌지 않는다. 이 뒤에 서열을 한 줄로 직접 붙여넣지 않는이상 primer3_core는 절대로 작동하지 않았다. 그리고 웹을 아무리 뒤져도 input.txt에 염기서열파일의 path를 넣는다는 설명도 일절 찾을 수가 없었다. 따라서 input.txt를 만든 뒤 data/NC_045512.softmask.one-line.fa.txt을 붙여 넣은 다음, SEQUENCE_ID와 SEQUNCE_TEMPLATE를 수정해야만 했다. input.txt 파일에 one-line FASTA file을 추가하면서 수정을 하는 간단한 스크립트를 작성하여 사용하면 좋을 것이다. 다음은 실제로 사용한 input.txt이다. /usr/share/doc/primer3/examples/example 파일과 비교하였을 때 수정 또는 추가한 라인은 *로 표시하였다. 맨 끝의 '='을 빠뜨리면 안된다. 그리고 PRIMER_THERMODYNAMIC_PARAMETERS_PATH는 굳이 지정할 필요는 없을 것 같다. **PRIMER_LOWERCASE_MASKING**  옵션은 3'-말단에 소문자가 있는 프라이머를 배제하게 한다. 우분투에서 primer3 패키지를 설치했다면 PRIMER_THERMODYNAMIC_PARAMETERS_PATH는 따로 설정할 필요가 없을 것이다.
  
  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_EXPLAIN_FLAG=1
  *PRIMER_LOWERCASE_MASKING=1
  *PRIMER_THERMODYNAMIC_PARAMETERS_PATH=/etc/primer3_config/
  *PRIMER_NUM_RETURN=100
  *PRIMER_PRODUCT_SIZE_RANGE=150-200
  *SEQUENCE_ID=NC_045512
  *SEQUENCE_TEMPLATE=ATTAAAGGTTTATAC...
  =
primer3는 다음과 같이 실행한다. 결과 파일(output.txt)에는 이해하기 매우 쉬운 구조로서 한 줄에 하나씩의 정보가 'TAG=VALUE'의 형식으로 기록된다. 
  $ primer3_core < input.txt > output.txt
  $ cat output.txt
  ...
  PRIMER_PAIR_NUM_RETURNED=100
  ...
  PRIMER_LEFT_0_SEQUENCE=TGATGGTGGTGTCACTCGTG
  PRIMER_RIGHT_0_SEQUENCE=GGCACGACAAAACCCACTTC
  PRIMER_LEFT_0=8700,20
  PRIMER_RIGHT_0=8864,20
  ...
  PRIMER_PAIR_0_PRODUCT_SIZE=165
  ...
output.txt 파일을 열어서 소문자가 없는 프라이머쌍을 선택하여 FASTA 파일로 저장한다. 마우스로 복사하여 편집하는 것은 매우 번거로우므로 다음과 같이 awk로 처리하면 된다. 프라이머 ID인 'PRIMER_LEFT_0'은 'primer-0-forward'의 형식으로 바뀐다. 프라이머쌍을 구별하는 숫자 정보는 프라이머 이름에서 중간을 차지하는 것이 바람직하다고 생각한다.
  $ grep '^PRIMER.*SEQUENCE' output.txt |\
    awk -F= '{\
             if(/LEFT/) dir="-forward"; else dir="-reverse";\
             split($1, a, "_"); id = ">primer-" a[3] dir;\
             printf("%s\n%s\n", id, $2)
             }' > p.fa
  $ head -n 4 p.fa
  >primer-0-forward
  TGATGGTGGTGTCACTCGTG
  >primer-0-reverse
  GGCACGACAAAACCCACTTC
==== 3. Evaluate primers ====
마지막 단계에서는 선별한 프라이머가 인간 유전체 등 nontarget sequence에 대하여 증폭산물을 생성하는지 점검하여 특이성이 우수한 것을 최종적으로 고르도록 한다. 이 과정에서는 mfeprimer 명령어가 쓰일 것이다. SARS-CoV-2, 인체 및 인플루엔자 유전체 데이터에 대한 인덱스를 먼저 생성해 놓는다.
  $ mfeprimer index -i data/test_1000.fasta
  $ mfeprimer index -i data/hg19.fa
  $ mfeprimer index -i data/influenza.fna
mfeprimer를 이용하여 SARS-CoV-2 database에 대한 프라이머(p.fa)의 coverage를 산출한다. EMBOSS primersearch를 사용했던 것에 비하면 훨씬 나은 것 같다. 그런데 primersearch에서 적용한 mismatch(0, 5, 10%)에 상응하는 파라미터는 무엇일까? 아래의 사례에서 -S는 maxSize(default 2000)을 의미하고 -t는 amplicon의 Tm cutoff(default 30)을 의미한다. 또한 '%%--%%virus' 옵션(
to check these primers amplify the percent of sequences in db, like SARS-Cov-2')의 의미도 알기가 어렵다. viral DB를 이 명령행에서 제공하지 않기 때문이다. 기본적으로 선택하는 DB가 있단 말인가?
  $ mfeprimer -d data/test_1000.fasta -i p.fa -S 300 -t 55 --virus -o p.mfe.txt
결과물인 p.mfe.txt에서 Hairpin List 및 Dimer List에 어떤 프라이머가 있는지 한번 확인해 보라. 가장 중요한 수치는 Primer Coverage에 관한 것이다. 이 샘플의 경우 검출하지 못하는 virus genome은 없었다.
  Primer Coverate Rate (virus mode, e.g., to check whether these primers can amplify all the variants of SARS-CoV-2 genomes)
  
  Coverage 100.00%, total records in database: 1000, amplicon number: 162663, missed number: 0.
검출에 실패하는 바이러스가 있다면 그 ID가 p.mfe.txt.virus.Failed.txt 파일에 기록된다. 

사람과 인플루엔자 DNA에 대해서도 마찬가지 방법으로 점검을 한다.
  $ mfeprimer -d data/hg19.fa -d data/influenza.fna -i p.fa -S 300 -t 55 -o p-against-human-influenza.mfe.txt
===== COVID-19 정보 관련 웹사이트 =====
  * [[https://www.gisaid.org/|GISAID]]
  * [[https://www.ncbi.nlm.nih.gov/sars-cov-2/|NCBI SARS-CoV-2 Resources]]
  * [[https://covid19datahub.io/|COVID-19 Data Hub]]
  * [[https://www.who.int/en/activities/tracking-SARS-CoV-2-variants/|WHO Tracking SARS-CoV-2 variants]]: VOC, VOI 공식 정보