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

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

2021년 진공관 앰프 제작을 위한 정보 수집

세상은 넓고 고수는 많다. 인터넷을 뒤져서 단편적으로 얻은 지식에 내 생각을 약간 달아서 블로그에 적는다고 하여 세상의 지식이 늘어나는데 보탬이 될 일은 없을 것이다. 그저 내 기억력을 보조하기 위한 수단으로 기록을 할 뿐이다.

지금까지 예닐곱대의 진공관 앰프를 다루어 보았다. 여기에는 외부에 제작을 의뢰한 것, 알리익스프레스에서 완제품으로 구입한 프리/헤드폰 앰프가 전부 포함된다. 완성품 보드를 구입하여 주변 부품만 붙인 것도 있고, 빈 PCB만 구해서 부품을 손수 납땜하여 만든 것도 있으며(6LQ8 PP & SE), CAD를 이용한 상판 설계와 point-to-point wiring으로 전 과정을 만든 것도 있다. 매우 낮은 출력이지만 푸시풀 앰프도 있었다.

내년에는 '제대로 된' 푸시풀 앰프를 point-to-point 결선 방법으로 만들고 싶다. 기왕이면 고정 바이어스 회로를 채택하여 서로 매칭이 되지 않은 출력관이라 하더라도 무난하게 쓰도록 만들자는 목표를 삼아 보았다. 아직 출력관은 정하지 않았다. 갖고 있는 PCL86을 그대로 쓸 수도 있고, EL84나 6V6도 고려 대상에 들어간다. EL84 또는 6V6은 완제품 혹은 반제품 형태의 PCB가 알리익스프레스 등에 흔하게 팔린다. 이걸 그대로 사서 쓴다면 만드는 재미는 반감될 것이고, 바이어스 조정도 힘들어진다.

푸시풀 앰프를 만드는데 필수적인 요소인 위상분리기(phase splitter) 또는 위상반전기(phase inverter)의 작동 원리를 이해하는 것도 어렵지만, 요즘은 고정 바이어스 공급 회로의 오묘함에 더욱 매력을 느끼는 중이다.

필요한 정보를 찾아서 목록으로 만들어 보았다. 특히 세상에 없던 6LQ8 사용 오디오 앰플리파이어가 자작인들을 통해서 2010년대 중반에 등장하게 되었으므로 그 개발 뒷이야기를 알아두는 것이 유익할 것이다.

[제이앨범] 작업실 > 앰프 자작 > 11LQ8/6LQ8 카테고리의 모든 글

[6LQ8 PP 관련 회로도 자료]  

[다음 블로그: 끝없는 여정] 6V6 푸시풀 파워앰프 - 회로도

[소리전자] PCL86 PP ('시나브로') PCB형 인티앰프 일반형 부품키트 - 회로도 참고용. 어떤 타입의 위상분리회로를 쓴 것인지 상당히 헷갈린다.

[Franz Wichlas] 4 x PCL86-PP - Amp - Jogis Röhrenbude의 웹사이트에 게시됨

[DIY Audio] New build of P-P PCL86 Amp - 로그인 필요. 회로도는 맨 아래 포스팅 참조

[DIY Audio] PC86 - Worth a try?

[The Tube CAD Journal] Phase splitters

[The Bona's HomePage] Phase Splitter - 매우 간결하게 설명을 잘 하였다. 이에 의하면 제이앨범이 디자인한 6LQ8 PP 앰프는 cathodyne phase splitter를 사용한 것이다.

[Audio Asylum] ECC802S SRPP / EL84(6BQ5) Push-Pull Amplifier - 25R 가변저항을 이용하여 두 출력관 사이의 바이어스 전류 균형을 맞춘다. 이런 용도로 사용할 가변저항의 와트 수는 어느 정도가 적당할까? 예를 들어 이런 제품의 규격은 0.5W이다. 2W를 견디는 다회전 볼륨은 훨씬 비싸다.

[Atra-Audio] PP Fixed Bias circuit design and calculator

[The Valve Wizard] Bias Supplies

[VTADIY: Vacuum Tube Amplifiers DIY] 5.3 Power supply for the fixed grid bias of a vacuum tube amplifier

2021년 11월 23일 화요일

fIDBAC: A Platform for Fast Bacterial Genome Identification and Typing 공부하기

유전체 서열 정보를 이용한 미생물 균주의 정확한 동정은 기초 생명과학과 감염병 관련 연구 분야에서 매우 중요한 일이다. 종(species) 또는 더 세분화된 수준으로 균주를 동정하거나 - 후자를 보통 타이핑(typing)이라고들 한다 - 유전체가 보유한 항생제 내성(AMR, antimicrobial resistance) 및 병원성 인자(VF, virulence factor)를 예측하는 것도 중요하다. 만약 유전체 해독을 실시한 미생물 균주가 쓸모있는 이차대사물을 만들어내는 것 같다면, 이를 책임지는 유전자 클러스터를 예측하는 것도 필요하다.

이러한 분석을 한꺼번에 실시하는 '입맛에 딱 맞는' 파이프라인은 아직 눈에 잘 뜨이지 않는다. 다만 병원성 미생물의 동정과 관련된 도구들은 가끔씩 눈에 보여서 이를 설치하고 테스트하면서 나름대로 평가를 하고 있다.

Frontiers in Microniology에 최근 공개된 fIDBAC라는 분석 플랫폼을 요즘 열심히 테스트하고 있다. fIDBAC는 면밀한 검토를 거쳐서 작성된 type strain 유전체 DB(N=12,784; 전체 목록은 12784.genome.txt를 참조)을 바탕으로 하여 미지의 유전체가 주어졌을 때 이것이 어떤 종에 속하는지를 정확히 판별해 주는 것을 목표로 한다.

다음 그림은 논문에서 제공한 graphical abstract이다.
fIDBAC의 프레임워크. 출처 링크
fIDBAC 개발에서 가장 심혈을 기울인 것은 type strain의 유전체 DB curation이었을 것이다. 논문을 보면 NCBI에서 균주의 source에 대한 정보를 참고하여 유전체 정보를 수집한 뒤, LPSN과 Bergey's Manual 및 IJSEM 논문을 참고하여 입증된 종명 및 type strain 정보를 완성했다고 한다. 이명(synonym) 정보를 정리하는 수고스러운 작업도 포함하였다. 

여기까지는 이름이 제대로 되어 있는지를 점검하는 것에 지나지 않는다. 다음 단계에서는 해당 균주의 유전체가 얼마나 온전한지를 확인하는 것이 필요하다. CheckM으로 점검하여 contamination이 5%를 넘거나 completeness가 90% 미만인 것은 버리고, 16S rRNA sequence를 추출하여 LTP database와 비교하여 genus level에서 불일치하는 것은 제외하였다.  

SILVA, RDP, GreenGenes DB는 많이 들어 보았지만 LTP('The All-Species Living Tree' Project)는 생소하여 검색을 해 보았다. 이곳에서 LTP의 상세한 연혁 소개와 함께 다운로드 기능도 제공한다. 앞서 소개한 잘 알려진 ribosomal RNA database가 있음에도 불구하고 LTP가 별도의 프로젝트로서 존재해야만 했던 이유가 있었을 것이다.

큐레이션의 마지막 단계에서는 모든 유전체에 대한 pairwise ANI 계산을 실시하여 동일한 종명을 갖는 것끼리를 클러스터링하였다. 여기에서는 SciPy 파이썬 패키지를 사용하였다고 한다. 이렇게 만들어진 type strain의 유전체 염기서열은 KmerFinder를 통해서 k-mer (다운로드 링크)DB로 변환하였다. KmerFinder에 대한 상세한 정보는 관련 논문을 참조하라.

Reference DB 구축이 완료되었으니, 이를 활용하여 query genome을 동정하는 방법을 수립하는 것이 fIDBAC platform의 나머지 절반 과정이 되겠다. 메인 스크립트(fIDBAC.pl)을 통해서 작동 순서를 알아보면 다음과 같다. 이는 논문에서 소개한 순서와는 조금 다르다.

  1. Reference genome의 k-mer DB에 대하여 검색 실시(findTemplate without -w). 결과 파일은 $sample.out.txt이다. 
  2. Genome annotation 실시(prokka) 후 16S rRNA와 유전자 서열 추출. 아미노산으로 번역한 유전자는 항생제 내성 및 병원성 인자를 탐색하는데 쓰임(AR_VF.run.pl)
  3. 16S rRNA에 대한 검색 및 fastANI를 이용하여 query genome의 종 동정 정보를 fastANISpecies라는 파일에 한 줄로 출력. 이것은 select_kmerfinder_16SAndANI.pl 스크립트를 통해서 이루어지는데 내부 구조가 좀 난해하다. [1]번 과정의 결과 파일($sample.out.txt)로부터 top1/3/10/20의 목록을 취한다. 최고 스코어를 보이는 genome이 3개라면 top1 목록은 세 줄이 된다. FASTA file 이름으로부터 목록 파일을 참조하여 taxonomic name을 얻는 과정이 꽤 된다. 16S rRNA 서열은 16S.V3.fa(관련 파일은 전부 여기에 있음; 이것이 LTP3132_SSU에서 유래한 것으로 추정됨)에 대하여 blastn을 실행하여 top 20을 얻은 뒤, local FASTA file에서 hit에 해당하는 것을 찾아 놓는다. fastANI에 계산에 투입되는 top 20은 16S 분석과 KmerFinder에서 얻어진 top genome을 추려서 얻어지는 것이 아닐까 생각하는데, Perl 스크립트를 뜯어보면서 정말 그러한지를 확인하지는 못하였다. 논문에서는 "The top 20 closest species were extracted from the results of the above-mentioned methods(즉 16S rRNA 서열 분석과 k-mer 분석)."이라고 하였으니 두 가지 결과를 종합하는 것이 맞을 것이다. 단, 두 가지 분석을 통해 얻어진 top 20 closest genome 안에서 종명이 consistent하지 않은 것이 섞여있을 경우 어떻게 처리하는지는 아직 잘 모르겠다.
  4. OrthoANI.all_tre.new.py 스크립트를 사용하여 ANI.pl 실행. 이 Perl wrapper script는 엄밀히 말해서 천랩이 개발한 OrthoANI와는 다른 것인데 왜 이런 이름으로 부르는지를 모르겠다. 이 wrapper script는 아직 잘 작동하지 않는 것 같다. 왜냐하면 주요 결과 파일인 OrthoANI.txtd와 ANI.all.txt가 빈 상태이기 때문이다. 왜 그런가 추적을 해 본 결과 이유는 너무나 허무한 곳에 있었다. OrthoANI.all_tre.new.py 스크립트 내부에 ANI.pl의 경로가 지정되어 있었다. 아니, config_db.txt 파일은 뭐하러 만들었단 말인가? 부글부글...
  5. ANIcaculator를 사용하여 gANI 계산. 유전자는 prokka가 만든 것을 이용한다. Step [4]와 [5]의 상호 연관성은 아직 파악하지 못하였다. ANI > 95%인 것만을 선별하는 것은 fastANI 레벨인가, gANI 레벨인가?
  6. MLST analysis
  7. SNP and MST analysis
  8. ....어휴!
스크립트를 뜯어보면서 작동 순서를 정리하려다가 머리에 쥐가 날 지경이 되었다. 이 목록은 앞으로 시간이 나는대로 틈틈이 업데이트를 해 보겠다. fIDBAC을 로컬 머신에 설치하면서 약 이틀 이상을 삽질(?)한 것을 생각하면 저절로 주먹에 불끈 힘이 솟는다. 설명은 너무나 부실하고, 필수 스크립트가 GitHub 사이트에 빠져 있고, 일반적인 tab-delimited file을 .xls 파일이라 부르는 등 사용자로서는 다소 과격한 의견이 되겠지만 프로그램 배포의 기본을 갖추고 있지 못하기 때문이다. MSTgold(MST for 'minimum spanning trees') 프로그램도 설치해야 fIBAC이 돌아가는데, 전혀 설명이 없다. 

아주 작은 사례를 보자. 사용자가 직접 수정해야 하는 config_db.txt 파일을 제공하고 있기 때문에 이것만 고치면 fIDBAC.pl 스크립트가 돌아갈 것으로 착각하기 쉽다. 그런데 엉뚱하게도 개별적인 스크립트 안에는 필수 프로그램들이 다른 경로로 지정되어 있는 것이 아닌가. 그것도 개발자 환경의 full path로 말이다. 테스트 러닝을 하면서 실행이 안되는 프로그램을 찾아 스크립트 내부를 수정하고, 또 실행 후 수정을 반복하고... 몇 차례의 시행 착오를 거쳐서 최초의 결과물을 손에 쥘 수 있었다. 그리고 본 블로그에서 자세히 설명하지는 않았지만, fIDBAC 개발자가 제시한 12,784 유전체 염기서열을 GenBank에서 다운로드하는 것도 다소 번거로왔다. 개발자는 accession number만 제시했을 뿐, 분석 과정에 필요한 유전체 및 유전자 염기서열을 따로 준비하여 목록 파일을 만들어 두어야 한다. 목록 파일(샘플)을 만드는 설명도 부실해서 애를 먹었고, 이들이 제시한 유전체 중 이미 수십 개는 GenBank에서 탈락하거나 업데이트된 것이 있어서 이를 다시 찾아야 한다. Genome assembly의 업데이트는 GCA_#########. 1뒤의 숫자가 하나 증가하는 것(1 -> 2)으로 표현되는 것도 있으나. 어떤 것은 아예 다른 accession no.로 대체되는 것도 있어서 이를 파악하려면 수작업에 의존하지 않을 수가 없었다.

이렇게 불평을 하기에 앞서서 과연 나는 다른 연구자들에게 도움이 될 수준의 완성도를 갖춘 스크립트를 공개한 적이 있었나? 블로그를 통해서 팁 수준의 짤막한 코드 조각을 이따금씩 올렸을 뿐이다. 남을 비판하기에 앞서서 나부터 반성을 하자.

논문으로서는 좋은 설계 개념을 제시하였지만 프로그램을 가져다 쓰려는 사람에게는 불편한 것이 사실이다. 그런 사람을 위하여 웹사이트를 구축해 놓았으니 꼭 필요한 사람은 이를 사용하면 된다. 앞으로 신종 박테리아는 계속 나올 것이고, 이를 얼마나 충실하게 업데이트를 할지는 지켜보아야 할 일이다. 업데이트는커녕 2년 정도 유지되다가 슬쩍 사라지는 웹사이트도 많이 있기 때문이다.

치명적인 약점(?) 또는 개선할 점을 생각해 보았다. 우선 archaea는 레퍼런스 DB에 전혀 반영되어 있지 않다는 점을 들고 싶다. Archaea는 인체에 병을 일으키는 것이 없다고는 하지만, 일반적인 species 동정을 목적으로 fIDBAC을 이용하려는 사람에게는 치명적인 단점이다. 만약 고균을 포함한다면 fIDBAC은 fIDBAAC(...Bacterial and Archaeal...)이 되어야 할 것이다. k-mer DB는 보다 현대적인 Mash 기반 DB로 바꾸는 것이 바람직할 것이다. 12,000개 정도의 유전체 서열을 다운로드하여 k-mer DB(maketemplatedb 명령)을 어젯밤부터 실행했는데 아직 400개도 진행을 하지 못했다. 파이썬 2.7을 써야 하고(버전 3에서 작동하도록 고쳤다는데 잘 안됨), 다중 쓰레드도 쓰지 못하니 너무나 답답하다.

그래도 fIDBAC을 뜯어보면서 얻은 성과가 있다면 ReadSeq이나 phylip 등 고전적인 프로그램을 다룰 기회가 생겼다는 점이다. 








2021년 11월 17일 수요일

43 오극관 싱글 앰프의 개선 작업 마무리

전원 트랜스포머의 교체 및 B+ 전압 조정 작업이 얼추 끝났다고 생각하고 음악을 듣는데, 좌우 채널 전체에서 '버석 버석'하는 소리가 들리기 시작했다. 접촉 불량인가, 새로운 타입의 발진인가? 또다시 좌절감에 휩싸였다. 겨우 쌍삼극관(6N2P) 하나와 핀이 6개 달린 구형 power pentode 두 개를 가지고  점대점 결선(point-to-point wiring, 나는 하드 와이어링이라는 표현이 옳지 않다고 믿음)으로 만든 싱글 앰프 하나를 제대로 못 만들다니! 내가 X손이라니!

앰프를 뒤집어 놓고 이곳 저곳을 건드려 보았다. B+ 전원 공급을 위해 새로 만들어 넣은 기판을 움직일 때 이에 맞추어 잡음도 같이 발생하는 것을 눈치챘다. 만능기판에 납땜한 스크류 터미널의 접촉이 아무래도 문제인 것 같았다. 나사를 풀고 기판의 동박에 납땜으로 직결을 하였더니 잡음이 사라졌다.

차폐용 커버가 없는 일반(오디오용이 아닌) 전원 트랜스포머를 쓰게 되면서 전원부에서 유도되는 잡음이 이전보다는 약간 늘어난 것 같았다. 리플 제거 보드를 사용함에도 불구하고 험이 들린다는 것은 몹시 자존심이 상하는 일이다. 전원을 넣으면 전원 트랜스 자체에서도 미약하나마 울리는 소리가 난다. 아마도 유도 잡음을 유발하는 부품의 배치 문제일 것으로 여겨진다. 실용상으로는 별 문제가 없으므로 당장은 그냥 쓰기로 하였다.

잡음의 원인을 찾는답시고 출력관을 좀더 상태가 좋은 실바니아 것으로 바꾸어 끼웠다.
대단한 것을 이룬 것은 아니지만 진공관 앰프의 기본에 대하여 이해도를 높인 좋은 계기가 되었다고 생각한다. 아직도 앰프를 뒤집어 놓고 보면 손을 대고 싶은 곳이 한두 군데가 아니다. 하루 아침에 해결될 문제는 아니다. 빨리 끝내려고 할 일이 아니라 과정을 즐기는 것, 그것이 진정한 아마추어의 자세이자 특권이 아닌가 싶다. 

아마추어는 과정을 즐기고, 프로는 결과로 말한다.

2021년 11월 16일 화요일

Roary 'query_pan_genome'의 고약한 출력물

Roary: the Pan Genome Pipeline은 내가 밥 먹듯이 즐겨 사용하는 프로그램이다. 최신판은 2019년 11월 6일에 릴리즈되었으니 더 이상 개선의 여지도 없는 프로그램이 아닐까 한다. MUMmer와 같이 거의 수정도 거치지 않으면서 오랜 기간 사랑받는.

Roary의 부속 유틸리티인 query_pan_genome은 두 그룹 간의 교집합, 합집합, 또는 차집합에 해당하는 gene cluster를 추출하는 다재다능한 프로그램이다. 그룹 정보는 GFF 파일로 제공하면 된다. 예를 들어 그룹 one에는 있지만 그룹 two에는 없는 유전자를 찾고 싶다면, 다음과 같이 실행하면 된다.

$ query_pan_genome -a difference --input_set_one 1.gff,2.gff --input_set_two 3.gff,4.gff,5.gff

만약 각 그룹을 구성하는 GFF 파일이 많다면 명령행에 이를 일일이 타이핑하기가 곤란할 것이다. 그런 경우에는 fofn(file of file name)을 만든 뒤 paste 명령을 써서 쉼표를 사이에 두고 이어 붙이면 된다.

$ cat one.fofn
1.gff
2.gff
3.gff
$ one=$(paste -sd, one.fofn)

두 그룹에 대하여 이렇게 변수를 만든 다음, '--input_set_one $one'과 같은 형식으로 옵션을 주면 된다. GFF 파일이 현재 디렉토리에 없다면, fofn을 만들 때 find 명령어를 잘 이용하면 된다.

'query_pan_genome -a difference' 명령을 수행하면 'set_difference_'로 시작하는 11개의 파일이 생긴다. 그룹 one에만 있는 유전자에 대한 정보는 다음의 파일에 수록된다.

  • set_difference_unique_set_one
  • set_difference_unique_set_one_reannotated
  • set_difference_unique_set_one_statistics.csv
여기서 한 가지 주의할 것이 있다.  그룹 one에만 있는('unique') 유전자란, 그룹 one의 모든 멤버에 다 있는 것만을 포함하는 것이 아니다. 예를 들어 그룹 one이 10개의 균주(A, B, C ... J)로 이루어졌다고 하자. 그러면 set_difference_unique_set_one에는 A에만 있는 것, A와 B에 있는 것 등 별의별 녀석들이 다 섞여 있다. 만약 그룹 특이적인 유전자를 찾는 것이 목적이라면, A부터 J까지 모든 멤버에 다 있는 것을 파악해야 한다. 그러려면 위에서 나열한 세 개의 결과 파일 중 set_difference_unique_set_one_statistics.csv를 참조해야 한다.

여기에서 문제가 발생한다. 바로 gene ID가 바뀐다는 것이다. set_difference_unique_set_one의 첫 번째 컬럼은 Roary 실행 결과 파일에서 보편적으로 쓰이는 unique cluster ID를 유지하지만(cluster는 accessory gene도 포함, 즉 pan genome을 전부 아우름), set_difference_unique_set_one_reannotated에서는 GFF 파일을 참조하여 일부가 gene name으로 바뀌고, 이것이 set_difference_unique_set_one_statistics.csv로 그대로 전달된다.

실제 사례를 들어 보자. 오늘 아침에 query_pan_genome 스크립트로 Strepcococcus genus에 속하는 두 종(N=301, set one은 232개, set two는 69개)의 roary 결과를 투입하여 13683개의 set one-specific gene을 얻었다. 그런데 gene ID를 sort하여 uniq 명령어를 통과시키면 13553개만 남는다. reannotated gene ID에서 중복이 발생했다는 뜻이다.

set_difference_unique_set_one과 set_difference_unique_set_one_reannotated 파일의 나머지 컬럼에는 실제 이를 구성하는 유전자들의 ID가 있으니 그룹 자체를 건드리지는 않는다. 그러나 pan_genome_reference.fa 파일에서 set one에 특이적인 유전자(대표 서열로서)를 찾으려면 여간 번거로운 일이 아니다.

query_pan_genome 스크립트를 전혀 쓰지 않고 roary의 기본 결과 파일인 gene_presence_absence.csv를 R에서 조작하여 각 그룹 특이적인 유전자를 찾을 수는 있는데 여간 성가신 일이 아니다. 왜 query_pan_genome 스크립트는 다시 annotation을 하여 식별자를 유일하지 않은 것으로 만들어 버리는지 알 수가 없다. QueryRoary.pm를 내가 고칠 수 있는 수준의 일도 아니고...


2021년 11월 11일 목요일

43 싱글 앰프 튜닝 - 전원부 개선

전원 트랜스포머를 220V 50W 절연트랜스로 교체하고, 초단관을 12DT8에서 원래 구성이었던 6N2P로 바꾸었다. 그리고 각 진공관은 최적 B+ 전원 조건에서 구동하고자 실험을 시작하였다. 일반적인 진공관 앰프에서는 전원장치에서 출력트랜스포머를 우선적으로 연결하고, 저항으로 전압 강하를 하여 초단관의 플레이트로 보내게 된다. 그런데 구식 라디오에 쓰이던 43이라는 오극관은 요구하는 플레이트 전압이 매우 낮다. 최대 160V, 찌그러짐이 가장 적은 조건(2와트 출력)에서는 135V이다. 콘트롤 그리드에 걸리는 바이어스 전압 -15V를 감안허면 155V가 가장 적당한 수준의 전압이 된다. 그러나 드라이브단의 6N2P의 플레이트 전압은 ~250V 정도가 되어야 한다.

인터넷에 43 오극관의 데이터시트가 있다는 것이 얼마나 고마운 일인지 모르겠다. 작성일은 무려 1936년 6월 26일다.

지금까지는 43에 공급한 전압을 약간 낮추어서 6N2P에 공급하는 방식이었다. 이렇게 낮은 플레이트 전압에서 6N2P가 작동을 하지 않는 것은 아니지만, 아무리 들어 봐도 소리가 작은 것은 초단관의 동작 상태가 최적이 아닌 것에 그 원인이 있는 것 같아서 개선 작업에 착수하였다. 별 것 있겠는가? 전압을 거는 순서를 바꾸면 되지 않는가. 특히 이번 개선 작업에서는 리플 제거 보드에서 각 진공관으로 연결되는 전압 강하 회로를 채널별로 따로 만들어 보았다. 비교적 값이 비싼 국산 전해 캐패시터가 많이 들어가서 아깝다는 생각이 들었다. 테스트가 끝난 다음에 기판에서 떼어내어 재활용하기도 어려운데...

보유한 저항의 값이 다양하지 않아서 다소 이상한 모습의 회로가 되었다. 리플 제거 회로의 출력은 260V 정도이다. 2.2K옴 저항을 거치면 6N2P 부하저항에 200V 정도가 걸리게 된다. 여기까지는 아주 마음에 드는데, 43 오극관에 맞게 전압을 더 낮추는 것이 문제이다. 6.8K옴 시멘트 저항 3개를 병렬로 연결하여 2.267K옴을 만드니 드디어 120V 정도가 된다. 이보다는 약간 높은 값을 얻어야 하지만 더 이상 병렬로 이어붙일 저항이 없다.

저항의 발열을 생각하면 말이 안되는 테스트 회로이다. 2.2K옴은 겨우 2와트급이다. 3개를 병렬로 연결한 6.8K옴 시멘트 저항(각 5와트)도 매우 뜨거운데, 2와트 금속피막저항 하나가 펄펄 끓는 것은 당연하다. 대략 계산하면 2.2K옴 저항에서 약 2.2W, 시멘트 저항 3개를 병렬로 연결한 합성저항에서는 1.95W가 걸린다.  저항이 뜨거워지면, 여기에 연결된 전해 캐패시터도 덩달아 뜨거워진다. 이래서는 곤란하다. 

2.2K옴 + 1K옴 정도로 2단 구성을 하면 최종적으로 150V 가량을 뽑을 수는 있을 것이다. 문제는 적당한 저항이 없다는 것. 소리는 과거보다 더 커진 것이 확실하다.

처음부터 2차에 150V 정도 출력되는 전원 트랜스포머를 썼더라면 높은 전압을 내리기 위해서 덜 고생을 했을 것이다. 되도록 기성품으로 팔리는 일반 용도의 트랜스포머를 쓰려는 것이 취지였으니, 이런 고생을 해도 감내해야 한다.

이번 개선 작업을 통해 떼어낸 전원트랜스포머(2차 230V 120mA, 6.3V 1.6A)는 정말 무용지물이다. 진공관 앰프를 자작하면서 처음으로 주문 제작을 한 것이었는데, 전기만 통하면 정격을 넘는 것도 아닌데 요란하게 울어대니 도대체 오디오 앰프로 쓸 수가 없는 것이었다. 그래서 전단에 110V 강압트랜스(단권 트랜스라서 용량에 비해 크기는 작음)를 달아서 출력 전압을 낮게 만들어 43 앰프에 사용하였던 것이다. 겉으로 보면 알 수가 없지만, 앰프를 뒤집어 놓으면 아주 가관이다. 
전기만 흘리면 우는 전원 트랜스포머. 그림 출처는 내 블로그의 글이다(링크).
새 앰프를 만들어 보려고 알리익스프레스를 며칠 동안 뒤지다가 기존의 앰프를 개선하는 것부터 먼저 해결하기로 마음을 바꾸어 먹기를 참 잘 한 일이라고 생각한다. 이렇게 단순하고도 아름다운(!) 앰프의 기본적인 문제도 해결하지 못한 상태로 또 다른 앰프를 만든다는 것은 올바른 배움의 자세가 아니다.