요즘 SARS-CoV-2의 알려진 변이체 정보를 분석하여 변이 위치를 수집한 뒤 conserved region에 해당하는 영역에서 바이러스 검출용 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/
여기에서 VCF를 BED로 전환하는 awk 명령어가 마음에 들지 않았다.
- [VCF] NC_045512 204 . G T #나머지 컬럼 생략
- [BED] NC_045512 203 204
204번째 위치의 G가 T로 바뀐 단일염기치환을 BED 파일에서는 왜 (203, 204)라는 위치로 표시하는 것일까? 'awk 명령어에 오류가 있는 것 같다!'는 것이 내 생각이었다. BEDOPS의 vcf2bed라는 유틸리티를 사용하면 정확한 feature 구간 위치 정보의 전환이 이루어질 것 같았다. 그런데 결과는 같았다.
왜 이렇지?
그건 내가 무식했기 때문이었다. 일단 유전체 관련 파일 포맷에서 위치 정보를 기술하는데 0-based 혹은 1-based system 중 어느 것을 쓰는지 잘 알아 두어야 한다. GFF/BAM/VCF는 1-based이지만 BAM/BED는 전산인이 좋아하는 0-based system을 사용한다. 이를 설명하는 글을 몇 개 찾아보았다. 전부 10년 가까이 된 글이다. 이를 아직까지 모르고 있었다니!
0-based system(+ half open)은 프로그램 입장에서 길이를 계산하기가 더욱 쉽다. 겹침(overlap)의 확인도 더욱 쉽다고 한다.
따라서 '인간적인' 1-based 좌표 체계에서 204번째 염기에 단일염기치환이 일어났음을 BED로 표시하려면 203, 204가 된다.
- 1-based: [start, end]
- 0-based: [start-1, end)
다음은 vcf2bed의 공식 설명에 웹문서에서 가져온 것이다. 대단히 부끄럽지만 오늘이라도 무식에서 벗어난 것을 다행으로 생각하자.
댓글 없음:
댓글 쓰기