2021년 12월 7일 화요일

6LQ8 싱글 엔디드 앰플리파이어의 보수

자작 6LQ8 푸시풀(push-pull) 앰플리파이어에 이어서 같은 진공관을 사용한 싱글 엔디드(single-ended) 앰플리파이어의 보수를 마쳤다. 이번에는 딱히 문제점이 발생한 것을 고친 것은 아니었다. 외부에서 12V 어댑터를 연결하여 히터의 전원을 공급하던 것을 앰프 본체 내부에서 직접 공급하게 만든 것이다. 이를 위하여 여분의 전원 트랜스(0/9/12/15/18V 1.2A)와 정류회로기판을 사용하였다. 6LQ8 데이터시트에 따르면 히터 전원은 6.3V 775mA가 되어야 하므로, 좌우 채널을 이루는 두 진공관의 히터를 직렬로 연결한 상태에서는 추가로 장착한 전원 트랜스의 용량(1.2A)이 부족하지 않다. 0-15V 탭에 정류회로를 연결하여 히터를 점화한 경우 최종 출력 전압은 12.5V 정도가 나와서 매우 적당하였다. 시판되는 직류 12V 어댑터를 사용하면 정격 수치인 12.6V의 95.2% 정도라서 문제는 되지 않는다. 일반적인 경우 +/- 5% 정도의 오차는 허용하기 때문이다. 물론 5%를 조금 넘어간다 해서 당장 진공관이 망가지거나 하지는 않는다.

DC로 만들지 않은 상태에서 12V 탭에서 히터 전원을 교류 상태로 공급해도 되지만, 그라운드 처리가 까다로와진다. 가장 바람직한 것은 전원 트랜스 2차의 센터탭을 앰프 회로의 그라운드에 한 점에서 연결하는 것이다. 그러나 이 트랜스에는 센터탭으로 쓸 만한 것이 없다. 만약 0-6-12V..의 전원 트랜스였다면 6V 탭을 쓰면 된다. 탭이 없으면, 저항 두 개를 이용하여 가상 중점을 만들 수도 있다. Valve Wizard의 Heater/Filament Supplies 문서에서 "Artificial Centre Tap" 항목을 참조하라.

트랜스 쪽에 적당한 탭이 있다 하여도 이 앰프의 PCV 자체에 한계점에 의해서 쓰기 어려운 사정이 하나 있다. 보드 내에서 히터 공급선의 하나와 그라운드를 연결하여 쓸 것을 전제로 만들었기 때문이다. 머리를 조금만 쓰면 직류를 만들지 않은 상태에서 가상 중점을 그라운드에 연결할 방법이 없는 것은 아니지만, 조금이라도 험을 줄이려면 직류가 더 나을 것이라는 생각에 예전에 만든 정류회로기판을 사용하기로 하였다.

두 번째의 개선 포인트는 두꺼운 4P 커넥터로 앰프 베이스로부터 바깥쪽을 빙 둘러서 상판으로 전원을 공급하였는데, 사실 이는 매우 거추장스런 모습이라서 상판에 구멍을 뚫고 내부에서 선을 연결하게 만든 것이다.

두꺼운 선을 여럿 연결할 때 나는 접속자를 즐겨 사용한다. DIY 오디오에서는 잘 쓰이지 않는 기법(?)이다.

무신호 상태에서 전원을 넣고 스피커에 귀를 가까이 대면 약하게 험이 들린다. 특별히 예민하게 생각하지 않는다면 음악 감상에는 큰 문제가 없는 수준이다. 전원쪽에는 잔류 리플을 제거하기 위한 특별한 조치를 취하지는 않았다. 예전에 전원 트랜스의 코어를 빼서 갭을 주어 싱글 엔디드 앰프용 출력 트랜스로 개조한 것이 남아 있었는데, 이를 초크 코일 대신으로 삽입하였다. 없는 것보다는 나은데, 제대로 만든 초크 코일을 사서 써 본 일이 없으니 그 효과가 어느 수준인지를 알기가 어렵다. 여력이 생긴다면 MOSFET를 이용한 노이즈 제거 회로를 만들어서 쓸 수도 있을 것이다.

출력 트랜스에는 폐 컴퓨터에서 적출한 냉각팬 관련 부품을 씌웠다. 여기에 포맥스나 아크릴 조각을 잘라서 덮으면 깔끔하게 마무리가 될 것 같다. 원래는 하얀 정육면체 모양의 덮개가 있었는데 방송 촬영 후 없어진 상태로 돌려받게 되었다. 히터 점화용 DC 어댑터도 없어진 것이 이번 보수 작업의 직접적인 계기가 된 것이다. 소품 대여와 관련하여 특별히 계약서를 쓴 것도 아니고, 약간의 대여료를 받았기 때문에 일부 부속품이 없어지거나 망가지더라도 어차피 감수할 생각이었다.

무척 작은 출력 트랜스를 사용한 소출력 싱글 앰프이지만 밤 늦은 시간에 침실에서 조용하게 FM 방송을 듣기에는 무리가 없다. 작업을 마친 뒤 사진 기록을 남겨 본다. 내가 생각해도 참 별난 모습의 진공관 앰프이다.




아직 나에게는 8개의 6LQ8 진공관이 남아 있다. 이것을 다 쓰려면 아마도 앰프를 하나 더 만들어야 될 것이다.


언제쯤 되어야 6V6이나 EL84와 같이 유구한 역사를 지닌 오디오 신호 증폭 전용관으로 앰프를 만들어 소리를 들어 보게 될까? 그렇게 될 날이 언젠가는 오게 될 것이다. 아마도 앞으로 2년 이내에?




2021년 12월 6일 월요일

VCF 파일을 BED로 전환할 때 좌표 정보는 어떻게 바꾸어야 하는가

요즘 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의 공식 설명에 웹문서에서 가져온 것이다. 대단히 부끄럽지만 오늘이라도 무식에서 벗어난 것을 다행으로 생각하자.



6LQ8 푸시풀 앰프의 잡음 줄이기

나의 6LQ8 푸시풀 앰프 '티라미수 II'는 방송 촬영 일정에 맞추어 급하게 만든 것이라서 음질적으로는 완벽한 튜닝을 하지 못한 상태였다. 약 반년이 넘는 대여 기간 후 돌아온 앰프를 점검해 보니 생각보다 잡음이 많았다. 특히 오른쪽 채널에서 무신호 시에 '웅-'하는 소리가 크게 들렸다. 처음에는 오른쪽 채널 보드 바로 밑에 위치한 히터 전원 공급용 SMPS에서 잡음이 유도된다는 생각을 갖고 있었는데 그건 아니었다. 왜 한쪽 채널에서만 잡음이 발생하는 것인가? 좌우 진공관을 바꾸어 끼워도 잡음은 여전히 오른쪽에서만 났다. 출력 트랜스에 문제가 생겼나? 그럴 가능성은 매우 희박하다. 부품을 바꾸어 연결하는 테스트를 해 보니 여전히 오른쪽 보드가 문제인 것 같았다.

부품이 조밀하게 꽂혀 있지만 비교적 단순한 PCB에서 무슨 문제가 발생한단 말인가? 저항이나 캐패시터와 같은 수동 소자에 무슨 문제가 생기겠는가? 제이앨범에서 새 PCB를 새로 구해서 처음부터 다시 만들어 볼까?

앰프 상판을 열어 놓은 상태로 이곳 저곳을 건드려 보았다. 간헐적으로 오른쪽 채널에서 천둥 소리 비슷한 굉음이 들린다. 입력단에서 접촉이 좋지 않아서 플로팅 상태가 될 때 흔히 들리는 소리다. 음량 조절용 포텐셔미터(흔히 '볼륨'이라 부르는 것)의 납땜 상태에 의심이 가기 시작하였다. 보드의 신호 입력용 커넥터를 스크류 드라이버로 단락시키면 아무런 소리가 나지 않는다. 이 관찰 결과는 앰프 보드 자체와 출력 트랜스까지는 문제가 없음을 시사한다.

포텐셔미터는 Alpha 제품이라서 품질 자체가 불량하다고 판단하기는 어려웠다. 포텐셔미터 전후를 연결하는 상태가 불량하다 판단하고 다음과 같은 몇 가지 개선 사항을 적용하였다.

  1. PCB의 착탈식 커넥터 중 신호 입력쪽은 납땜 직결한다.
  2. 입력 실드선을 교체한다. 낡은 인터케이블을 잘라서 사용했었는데 품질이 영 마음에 들지 않았다.
  3. NFB 배선에는 실드선을 사용하고 그라운드 처리도 빼먹지 않았다. 
  4. 포텐셔미터와 입력단 사이에 커플링 캐패시터와 470K 그리드 리크 저항을 삽입한다. 아래 사진에 보인 만능기판이 바로 이 기능을 한다. 이 기판은 43 오극관 앰프의 초기 실험에서 사용하던 것이다.
빨간색 캐패시터는 예전에 사용하던 TPA3116에서 적출한 것이다. 

입력단에 포텐셔미터가 위치할 경우 그리드 리크 저항은 종종 생략된다. 그러나 포텐셔미터의 접촉이 좋지 않을 때에 그리드에 바이어스 전압이 걸리지 못하게 되므로 문제가 생길 수 있다. 커플링 캐패시터는 입력단에 DC가 유입되는 것을 막는 역할을 한다. 비정상적인 신호가 들어오지 않는 이상 커플링 캐패시터는 없는 것이 가장 좋을 수도 있을 것이다.

결과는? 언제 그랬냐는 듯이 잡음이 사라졌다. 결국은 좋지 못한 선재와 납땜 불량이 원인이었던 것 같다.

선재의 재활용에는 특히 주의해야 할 것 같다. 한국진공관앰프동호회의 고 안병원 회장께서 작성하신 '왕초보 따라하기'를 보면 폐 컴퓨터에서 쓸만한 전선을 많이 구할 수 있다고 하셔서 나도 이를 따라 했었다. 그런데 이렇게 재활용한 전선 중에는 유난히 납땜이 잘 되지 않는 것이 있었다. 주석도금선처럼 은색을 띤 연선에 왜 이렇게 납이 붙지 않을까? 어쩌면 이런 전선 중 일부는 알루미늄을 사용하며 제조한 것이 아닌지 모르겠다. 알루미늄은 구리보다 전도도가 떨어지지만 가볍기 때문에 송전선 등에서 흔히 쓰이는 것으로 알고 있다. 요즘 전기차에서도 선재 총 중량을 줄이기 위해 알루미늄을 많이 쓴다고 한다. 그런데 알루미늄과 구리를 같이 꼬아서 연결하면 알루미늄이 부식하는 현상(갈바닉 부식, Galvanic corrosion)이 일어난다고 한다. 이는 전기 배선뿐만이 아니라 서로 다른 종류의 금속을 직접적으로 접촉해야만 하는 모든 현장에서 만나게 되는 문제이다.

유튜브: How to identify copper vs. alunimum wire | DIY material knowledge

우리나라의 DIYer가 알루미늄이 든 선재를 접할 일이 얼마나 있을지는 잘 모르겠다. 단지 나의 납땜 실력 부족을 선재 탓으로 돌리는 것은 아니었을까? 어쨌든 선재는 가급적 새 것을 사용하는 것이 바람직할 것이다.

납땜 개선 작업이 끝난 티라미수-II 앰프.
진공관 앰프의 매력은 따뜻한 불빛을 바라보며 아날로그 감성을 느낄 수 있다는 것.

다음 숙제는 6LQ8 '싱글' 앰프를 개선하는 것. 하부 케이스와 상판 사이의 배선을 바깥쪽에서 매우 굵은 4P 커넥터로 연결하던 종전 방식을 바꾸려 한다. 전원부 하나로 두 앰프를 번갈아 쓰던 예전 방식을 더 이상 고수할 필요가 없기 때문이다. 나무로 만들어진 상판에 구멍을 뚫고 케이블을 그리로 지나가게 하면 된다. 두 번째로는 히터 전원 공급 방식을 바꾸는 것이다. 기존에는 12V 직류 어댑터를 외부에서 접속하는 방식을 썼었다. 위에서 보인 푸시풀 앰프에서는 어댑터를 앰프 케이스 안에 넣어 버렸었다. 그러려면 어댑터 하나를 완전히 희생해야 한다. 마침 남는 전원트랜스가(9/12/15/18V 1.2A) 하나 있으니 이를 사용하려고 한다. 필요하다면 DC로 정류하여 사용할 수도 있다. 







2021년 12월 3일 금요일

설계된 프라이머의 결합 위치를 그림으로 표시해 보자

BioPerl의 그래픽스 기능(Bio::Graphics)을 사용하여 서열 'track' 위에 여러 feature를 표시할 수 있다. 프라이머 정보 역시 feature로 간주하여 방향까지 표시되는 화살표로 표현하면 된다. HOWTO 문서를 참조하여 직접 손으로 Perl 코드를 짜도 되지만, 남들이 먼저 만들어 놓은 것을 살짝 고쳐서 쓰면 더욱 낫지 않겠는가?

프라이머를 설계하고 시각화하는 공개 프로그램을 뒤적거리다가 이런 논문을 발견하였다.

PrimerView: high-throughput primer design and visulaization (2015).  GitHub 링크

미국 조지 워싱턴 대학의 Damien M. O'Halloran이라는 사람이 단독으로 개발하여 발표한 것이다. 이 도구는 리눅스에 설치하여 쓰는 Perl module이다. 여러 유전자를 multi-FASTA file로 제공하면 primer3가 각 유전자를 증폭하는 프라이머쌍을 설계하고, template DNA에 결합하는 위치를 이미지 파일로 보여주는 소박한(?) 도구이다. O'Halloran은 2016년에 GUI tool인 PrimerMapper(GitHub 링크)라는 것을 발표하기도 하였다. 기능은 훨씬 많아진 것 같은데 아직 제대로 테스트를 해 보지 못하였다.

PrimerView는 내부적으로 primer3_core를 사용하도록 되어 있다. Template DNA의 염기서열은 사용자가 제공하면 되지만 각 염기서열에서 설계하는 최대 프라이머쌍의 수 등 사용자가 입맛에 맞게 파라미터를 맞추어 넣을 수가 없다. 핵심 구성요소인 PRIMERVIEW.pm을 아무리 들여다 보아도 primer3 파라미터를 넣는 곳은 없다. 

그러면 소스를 고쳐 보면 되겠지... 모듈에서 대부분의 기능을 구현하고 메인 소스 코드는 간결하게 만드는 스크립트 작성 스타일에는 별로 익숙하지 않아서 시행착오를 거듭하여 수정을 해 보았다. 그런데 뭔가 이상하다. 아무리 모듈을 고쳐도 결과는 변화가 없다. 그거 참 이상하다. 이유는 간단한 곳에 있었다. 스크립트가 참고하는 모듈의 위치는 /usr/local/share/perl/5.30.0인데 프로그램을 설치한 곳의 lib/ 이하 위치를 참조하는 것으로 생각하여 엉뚱한 것을 고쳤던 것이다. @INC 배열로 지정된 곳 중 어디에 PRIMERVIEW.pm이 있는지를 찾고 나서야 비로소 답을 알게 되었다. Perl의 모듈 관리가 원래 그렇다. FindBin과 use lib 디렉티브를 잘 사용하여 보다 깔끔하게 고쳐야 되겠다.

원본 모듈에서는 사용자가 입력한 염기서열만을 primer3로 보내게 된다. 서열 정보를 담는 input file의 앞부분에 내가 별도의 파일로 제공한 primer3 파라미터 파일을 삽입하게 만들었다. 다음은 그 결과로 얻은 그림이다. 입력 template는 신종코로나바이러스의 게놈 서열을 주었고 총 20쌍의 프라이머를  만들게 하였다. Product size는 150-200 bp range로 하였다.

Primer Map이 그려주는 그림보다는 나은가? 해상도가 다르다는 것이 가장 큰 차이점이라고 생각한다.
프라이머 위에 적힌 숫자는 프라이머 ID를 대신하는 것으로서, 5'-end의 위치에 해당한다. 여기서 주의할 것이 하나 있다. Primer3의 위치 체계는 0-based라는 점이다. 그리고 한쪽 프라이머가 (거의) 동일한 상태에서 반대편 프라이머는 별도의 것을 취하는 복수의 amplimer를 예측하기도 한다는 점이다.

이 그림에서는 어떤 프라이머가 서로 짝을 이루고 있는지 한눈에 보이지는 않는다. 그러나 잘 관찰하면 어렵지 않게 파악할 수 있다.

PrimerView의 독특한 점은 설계한 프라이머가 결합하는 위치를 muscle alignment를 돌려서 찾는다는 것이다. 아니, primer3 결과 파일에 결합 위치가 다 있는데 왜 이런 수고를? 염기서열이 그대로 보이는 출력물을 이미지로 만들기 위해서 그런지도 모르겠다.

만약 다른 프로그램으로 프라이머를 설계한 뒤 결합 위치를 표시하고자 한다면 이러한 정보를 수록한 텍스트 파일을 먼저 만든 다음, PRIMERVIEW.pm의 _graphics_all_primers 함수를 참조하여 코드를 직접 만들어야 한다. 이것은 별 문제가 되지 않지만, 정보 파일을 만드는 일은 각자 알아서 해야 한다. EMBOSS의 primersearch가 되었든, 최근에 소개한 MFEprimer가 되었든 다른 프라이머 평가 도구를 써서 만든 결과 파일을 다듬어야 한다.


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'라는 글귀가 선명하다.

2021년 11월 25일 목요일

PubMLST 자료를 쉽게 가져오는 방법 - mlst_check(Sanger)

오늘도 나의 2021년 11월 하순을 뜨겁게 달구고 있는 fIDBAC과 관련한 이야기를 써 내려가려고 한다. fIDBAC 파이프라인에서는 16S rRNA 서열 분석과 KmerFinder를 통해서 top 20 closest genome을 고르고, 이를 대상으로 fastANI를 실시하여 일정 컷오프를 통과하는 결과로부터 query genome이 어떠한 species에 해당하는지를 판별하여 그 정보를 OUTDIR/fastANISpecies 파일에 출력한다. example/example.fa를 사용한 경우는 이러하다. 

$ cat example_2021_11_25_08_20_15/fastANISpecies 
Salmonella enterica

fastANISpecies 파일에 여러 라인이 기록될 수도 있을까? fIDBAC.pl(v20181126)의 133번째 줄에서 다음과 같이 fastANISpecies 파일의 첫 줄만을 뽑아내도록 한 것을 보면 가끔 그런 일이 생기기도 하는 모양이다.

my $template=`head -n 1  $result/fastANISpecies`;chomp $template;

fastANISpecies 파일을 만드는 과정에 대해서도 살펴볼 것이 많지만 이에 대해서는 나중에 알아보기로 하고, 오늘은 fIDBAC.pl이 이 파일을 이용하는 다음 단계에 관해서 탐구해 보자. 이 단계에서 fIDABAC.pl 스크립트는 표준 출력으로 다음을 내보낸다.

Salmonella enterica...Salmonella enterica...

마치 강도에 피습당한 피해자가 '나를 칼로 찌른 놈은 아무개였어...'하고 간신히 내뱉고는 숨을 거두는 영화의 한 장면같다. 똑같은 종명(Salmonella enterica)를 두 번에 걸쳐서 출력한 것으로 보이지만 앞의 것은 $name 변수이고 뒤의 것은 $template 변수이다. 어느 종에 대한 MLST 분석을 해야 하는지 결정하는 열쇠가 되는 것이 바로 $name이다.

config_db.txt에서 설정한 LPSN(Species.tax 파일)의 첫 번째 컬럼은 'Salmonella_enterica'와 같은 형식으로 되어 있다. 공백을 underscore로 치환한 것이다. 왜냐하면 MLST DB를 구성하는 각 종명이 디렉토리 형식으로 그대로 쓰여야 하기에, 공백을 그대로 두면 곤란하기 때문이다. 이 문자열은 최종적으로 $species 변수에 저장된다. 일단 Species.txt 파일에서 첫 컬럼을 하나씩 읽어서 underscore를 다시 공백으로 바꾸어 $name과 비교하고, 일치하는 것이 있으면 이를 OUTDIR/Taxonomy.txt에 기록한다. 이는 MLST 함수가 하는 일 중의 하나이다.

MLST 분석은 fIDBAC.pl의 199번 라인부터 시작한다. 

my ($species, $flag) = split (/\t/, MLST ($template));
 if ($flag == 0){
      my $mlstlist = glob "$pubMLSTdb/$species/*.txt" ;
      system("perl $Bin/run_MLST.pipeline.pl -s $species -t $mlstlist -q $input -o $result \n");
  }else{
      print STDERR "This species is not exists in the pubMLST db\n";
} 

여기까지 살펴보면 $template와 $name 변수가 제대로 쓰였음을 알 수 있다. pubMLSTdb가 설치된 디렉토리 하위에는 $species에 해당하는 서브디렉토리가 존재하고, 또 그 밑에 세부적인 파일들(allelic progiles & sequenes)이 있어야 한다. 정확히 말하자면 여기까지가 서론이었다.

config_db.txt에서는 pubMLSTdb의 설치 위치를 지정해 두어야 한다. 그리고 fIDBAC README.md 파일에 의하면 'Download PubMLST database here.'라고 하였다. 그러나 'here'를 클릭하면 XML 파일이 나올 뿐이다. PubMLST 공식 웹사이트의 다운로드 페이지를 가 보아도 전체 자료를 하나로 묶어서 클릭만 하면 쉽게 내려받을 수 있게 만들지 않았다. 흠... 공식으로 제공하는 API를 이용해서 활용할 재주는 나에게 없다. 어떻게 해야 오늘 날짜 기준으로 153개나 되는 데이터베이스를 한꺼번에 가져올 수 있을까? 'bactopia datasets --available_datasets'으로 확인하면 무려 683개나 되는 자료가 있다. 아마 후자의 숫자는 PubMLST가 호스팅하지 않는 다른 것까지 포함해서 그런지도 모른다.

Torsten Seemann의 mlst를 사용하여 다운로드한 자료를 쓸 수 있을까? 이는 직접 실행할 일은 별로 없고, TORMES를 통해서 이용하게 된다. 나의 경우에는 다음의 위치에 파일들이 전부 저장되어 있다.

/opt/miniconda3/envs/tormes-1.3.0/db/pubmlst

그러면 이 디렉토리를 config_db.txt의 pubMLSTdb로 설정하면 될까? 전혀 그렇지 않다. 이 디렉토리의 목록은 다음과 같다. fIDBAC.pl에서는 Salmonela_enterica라는 하위 디렉토리를 이로부터 찾으려고 할 것이지만, 실제 존재하는 디렉토리 명칭은 약칭인 senteriaca이다. 이는 pubMLST 체계에서는 공식적으로 쓰이는 문자열이지만 fIDBAC.pl은 완전한 종명을 원한다.


어쩌겠는가? 완전한 종명 체계로 pubMLST를 다운로드해 주는 유틸리티를 찾아 보아야 한다. 검색 끝에 Sanger Institute에서 개발한 mlst_check라는 것을 발견하였다. 설치는 쉽게 하였으나 역시 SSL 인증서 문제로 직장 내 전산망에서는 다운로드 명령('download_mlst_databases')이 먹히질 않았다. 이를 우회하는 가장 간단한 방법은 노트북 컴퓨터로 집에서 설치한 뒤 직장으로 가지고 나와서 복사하는 것이다. 

이제 fIDBAC이 잘 실행되겠지... 그러나 이는 잘못된 기대였다. T. Seemann의 mlst 명령을 이용한 것과 Sanger의 mlst_check를 이용한 것은 종명 서브디렉토리뿐만 아니라 그 하위의 파일 구조도 다르다는 것을 발견하였다. Salmonalle enterica의 예를 들어 보자.

  • T. Seemann의 스타일: senterica 디렉토리 하위에 aroC.tfa  dnaN.tfa  hemD.tfa...의 파일이 있음
  • Sanger의 스타일: Salmonella_enteria 디렉토리 하위에 profiles/profiles_csv와 alleles/alleles_fasta 두 개의 파일이 있음. 이는 fIDBAC에 의해 쓰이기가 곤란하다. 왜냐하면 '$mlstlist = glob "$pubMLSTdb/$species/*.txt'라는 코드를 상기해 보라.
지금까지 테스트한 바에 따르면, T. Seemann 스타일로 다운로드하되 종명 디렉토리를 Salmonella_enterica로 바꾸어 놓은 것만 제대로 결과를 내었다. 어느 것이 더 바람직한가? (1) T. Seemann 스타일로 다운로드한 뒤 senterica => Salmonella_enterica로 디렉토리명을 바꾸는 것과, (2) Sanger 스타일로 다운로드한 뒤 alleles_fasta 파일을 처리하여 각 유전자별로 aroC.tfa  dnaN.tfa... 파일을 만드는 것. 종명과 그 약자를 연결한 파일이 있다면 (1)의 작업은 아주 쉬울 것이다. 그러나 이는 간편한 텍스트 파일이 아니라 XML로 주어진다. 

이 문제를 해결한 뒤에도 복수의 scheme이 존재하는 종이라든가(예: Acinetobacter_baumannii_1 & Acinetobacter_baumannii_2), spp로 끝나는 scheme을 fIDBAC이 제대로 처리할지는 의문이다.

여기까지 적고 나니 회의감이 밀려온다. 내가 지금 월 하고 있는거지?

회의감을 정리하는 도중 또 경악을 금치 못할 일을 발견하였다. run_MLST.pipeline.pl 스크립트가 pubMLST가 지정된 위치를 별도로 지정하여 사용하는 것이다! config_db.txt는 어디다 갖다 버리고!

2021년 11월 24일 수요일

fIDBAC 실행용 펄 스크립트의 use compiler directive 선언 오류, 그리고 플러스 알파(+베타, 감마...)

fIDBAC(어제 작성한 소개의 글 링크)에서 오류를 찾는 것이 취미가 될 수준에 이르렀다. 나는 파이썬에는 까막눈이지만 Perl은 조금 아는 편이라서 fIDBAC의 주요 스크립트를 뜯어보는 수준은 된다. 메인 스크립트는 fIDBAC.pl이다.

fDIBAC.pl을 처음 실행하게 되면 GACP.pm 모듈을 인식하지 못한다는 에러가 발생한다. 

Can't locate GACP.pm in @INC (you may need to install the GACP module) (@INC contains: ...)

스크립트의 위치로 이동하여 실행을 하거나, 혹은 PERL5LIB 환경변수를 스크립트가 위치한 곳으로 선언하면 일단은 에러가 없어진다. GACP.pm은 config_db.txt 설정 파일을 읽어들여서 중요한 프로그램과 파일을 변수로 저장하는 역할을 한다.

왜 이런 사소한 오류가 뜨는 것일까? fIDBAC.pl 스크립트의 시작 부분을 확인해 보았다.

#!/usr/bin/perl -w
#author:liangq at 20181126 ,modified by liangqian ,at 20181206
use strict ;
use Getopt::Long ;
use Cwd 'abs_path';
use File::Basename ;
use File::Path 'mkpath';
use FindBin '$Bin';
use Data::Dumper;
use List::Util qw/max min/;
use lib $Bin;
use threads;
use GACP qw(parse_config);

특별한 문제는 없다. FindBin을 이용하여 원본 펄 스크립트의 위치를 찾아서 $Bin 변수에 저장하고(작은 따옴표로 둘러싼 것은 별로 마음에 들지 않음), 이어서 use lib $Bin 디렉티브를 선언했으니 펄 스크립트가 있는 디렉토리를 @INC의 맨 앞에 추가한 셈이 된다. 그리고 가장 마지막에서 use GACP를 선언하여 GACP.pm 모듈을 선언하였다. 훌륭하다!

그런데도 여전히 같은 에러가 발생하였다. GACP.pm을 로드하는 스크립트가 또 있나? fIDBAC.pl에 의해 내부적으로 구동되는 다른 스크립트 중에서 select_kmerfinder_16SAndANI.pl와 AR_VF.run.pl도 GACP.pm을 필요로 한다. 스크립트에 오타 같은 것은 없는데 왜 에러가 사라지지를 않는 것일까...

이런! select_kmerfinder_16SAndANI.pl 파일을 열어보니 'use lib $Bin' 라인이 없었다. 무슨 이런 실수를... 개발자는 항상 스크립트 설치 디렉토리에서만 테스트를 했단 말인가? 수준 이하의 오류를 찾아내어 수정하느라 인생을 낭비하는 것만 같아서 이제는 서글프기까지 하다. 빠진 라인을 삽입하였더니 이 오류는 사라졌다.

그러나 아직 끝이 아니다. 

  1. run_rgi.sh가 호출하는 format_RGIresult.pl는 도대체 GitHub 사이트 어디에 숨어 있는가?
  2. config_DB.txt에서 'orthoANI'라는 이름으로 부르는 프로그램은 도대체 무엇을 의미하는가? 내가 발견한 fIDBAC 파이프라인의 문제점 중 가장 심각한 것은 바로 여기에 있다. 
orthoANI의 문제를 좀더 상세하게 알아보자. config_DB.txt 파일을 열어보면 orthoANI는 ANI.pl 펄 스크립트를 지정하는 것으로 보인다. 그러나 OrthoANI라 하면, 보통은 천랩에서 개발한 ANI 계산용 알고리즘을 뜻한다. 이 알고리즘을 프로그램으로 구현한 것은 자바 애플리케이션인 OAT(Orthologous Average Nucleotide Identity Tool)이다. 'ortho-'라는 접두사를 붙여서 불필요한 혼동을 불러 일으키고 있다. 이것이 (2)에 따르는 첫 번째 문제이다.

ANI.pl은 원래 Jaipeng Chen이라는 사람이 JSpecies(GUI Java application)를 참조하여 legacy blast + Perl로 만든 ANI 계산용 스크립트이다(GitHub). 9년 전에 업로드된 상태 그대로 수정되지 않았다. fIDBAC의 script/Average_Nucleotide_Identity/readme.txt 파일에서도 이 GitHub를 인용하고 있으니 내 예상이 틀리지는 않을 것이다. ANI.pl이 호출되는 순서는 다음과 같다.

fDIBAC.pl(main script) -> OrthoANI.all_tre.new.py -> orthoAni.sh라는 스크립트를 실행 단계에 작성하여 활용함

OrthoANI.all_tre.new.py에서 orthoAni.sh 파일을 작성하기 위하여 다음과 같은 문자열을 만들어 ANI.pl을 실행하도록 만드는데, 이게 또 이상하다. 


ANI.pl의 필수 옵션인 '-fd formatdb -bl blastall'이 빠진 상태이다.  이 옵션은 formatdb와 blastall 실행 파일을 지정하기 위한 것이다. $PATH에 위치한다고 하여 생략해서는 안 된다. 따라서 그림에서 보인 cmd로는 ANI.pl이 제대로 실행되지 않는다. 이것이 두 번째 문제이다. 그러면 cmd를 조합하는 명령어 라인에 '-fd formatdb -bl blastall'을 삽입하면 되지 않을까 생각할 수 있다. 그러나 전혀 그렇지 않다. ANI.pl은 출력 파일을 쓰지도 않고, 오로지 표준 출력에 두 genome으로부터 계산한 ANI 수치를 표시할 뿐이다 OrthoANI.all_tre.new.py 스크립트의 후반부를 보면 ANI 수치 쌍을 전부 조합하여 하나의 OrthoANI.txt 파일을 만드는 것으로 되어 있는데, 내가 알고 있는 ANI.pl의 출력 형식으로는 이를 어떻게 만드는지 이해하기 어렵다. 이것이 세 번째의 문제이다.

full_path_to_file_1 VS full_path_to_file_2
 ANI: 93.0621988037596

세 번째 문제의 해결을 위해서 지금까지 미루어 두었던 파이썬 공부를 시작할 용의가 있다. 지금이 아니라면, 앞으로 영원히 파이썬 문맹으로 남게 될지도 모른다.