Fast and flexible bacterial genomic epidemiology with PopPUNK

Genomic epidemology에 관한 최신 논문을 찾아보다가 올해 1월 24일 Genome Research에 발표된 새로운 소프트웨어 도구인  PopPUNK(Population Partitioning Using Nucleotide K-mers)를 발견하였다.

Lees JA, Harris SR, Tonkin-Hill G, Gladstone RA, Lo SW, Weiser JN, Corander J, Bentley SD, Croucher NJ. Fast and flexible bacterial genomic epidemiology with PopPUNK. Genome Research 29:1-13 (2019). doi:10.1101/gr.241455.118

출처: https://poppunk.readthedocs.io/en/latest/index.html

병원체 감시(pathogen surveillance)를 위한 유전 데이터를 분석하는 방법으로는 MLST(multilocus sequence typing)가 전통적으로 널리 쓰인다. 유사한 MLST sequence type(ST)을 갖는 균주끼리 minimum-spanning tree를 만들어서 eBURST로 시각화하는 방법이 매우 보편적인데, 몇 개의 housekeeping gene만을 대상으로 하기 때문에 해상도에는 한계가 있다. 이제는 유전체 시퀀싱을 매우 염가에 할 수 있으므로, cgMLST(core genome MLST) 또는 wgMLST(while-genome MLST)를 활용하는 방법도 바람직하다. 그러나 결국 sequence alignment라는 시간이 오래 걸리는 계산을 피할 수 없고, 유전자에 대한 정보가 필요하며, 새로운 샘플이 약간만 추가되더라도 분석을 처음부터 다시 해야 하므로 스케일업이 매우 어렵다.

PopPUNK는 이러한 단점을 극복한 typing & clustering 도구로서 여러 시각화 도구에서 사용 가능한 결과 파일을 만들어 낸다. 두 유전체 사이의 관계(p_match,k)를 가변 k-mer sketch로 계산하여 core distance(π, SNP density)와 accessory distance(a)를 축으로 하는 2차원 평면에 점으로 표시하고, 이들의 분포를 가장 잘 설명할 수 있는 모델의 피팅 및 최적화를 거쳐서 네트워크 구조를 만들어 낸다.

위의 수식에서 보였듯이 한 쌍의 서열 사이에 k-mer가 매치할 확률은 비교 대상 중 어느 하나에만 존재하는 accessory locus에 존재하지 않을 확률(p_accessory) 곱하기 mismatch가 없는 공통 core genome 서열에 존재할 확률(p_core)로 표현된다. 나도 아직은 이 수식을 완벽하게 이해한 것은 아니다. PopPUNK의 전체 과정을 설명하는 그림을 원문 사이트에서 인용해 본다.

Fig. 1. Summary of the PopPUNK algorithm. 출처 링크.

99 균주의 유전체 서열을 가지고서 PopPUNK를 돌려 가면서 감을 잡기 위해 노력하는 중이다. 아직은 모델의 최적화를 이루지 못하였다. 논문에서 예시로 든 다른 species에 비해서는 다양성이 매우 적다고 느껴진다.

microreact로 시각화한 결과.

PopPUNK 결과물을 제대로 보려면

PopPUNK가 만들어내는 다양한 interactive output은 다음의 환경에서 보면 된다.

• Microreact: Open data visualization and sharing for genomic epidemiology
• phandango: Interactive visualization of genome phylogenies
• GrapeTree: A fully interactive, tree visualization program, which supports facile manipulation of both tree layout and metadata
• Cytoscape: Network data integration, analysis and visualization in a box (네트워크 구조를 상세하게 조사하려면 cytoscape를 권장한다)

[하루에 한 R] POCP 매트릭스를 이용하여 미생물 균주를 genus 별로 군집화하기

POCP(percentage of conserved proteins)에 대한 간략한 소개는 여기를 참조하라. 결론만 이야기한다면 POCP 50%를 기준으로 하여 세균이 같은 속(genus)에 속하는지를 판별하는 것을 제안한 것이다. 50%라는 기준은 평균값을 통해서 제안된 것이라 실제의 분포는 꽤 넓다.

POCP 매트릭스를 만들면 R에서 계층적 군집화(hierarchical clustering) 분석을 통해서 같은 genus에 속하는 균주를 구분하는데 도움이 될 것으로 기대를 하였다. heatmap을 그리면 균주들 사이의 대략적인 관계를 파악할 수 있지만, 실제로 군집이 어떻게 이루어지고 그 멤버는 누구인지를 결정해 주지는 않기 때문이다.

군집 분석에 대하여 알기 쉽게 설명한 글이 있어서 먼저 소개해 본다. 거리를 계산하는 방법과 군집 알고리즘의 차이에 대해서 특히 유념하여 읽도록 한다.

군집분석(Cluster Analysis) - Amazon AWS
[R 군집분석 (Cluster Analysis)] 군집분석의 개념 및 유형

그러면 34개의 유전체 서열을 이용하여 만든 POCP 매트릭스로 작업을 해 보자. 이 매트릭스는 d3이라는 데이터 프레임에 들어 있는데, 아직 완성되지 않은 일이라서 genus 정보만 남기고 균주 정보를 살짝 가리기 위해 약간의 트릭을 썼다. 수정된 데이터 프레임은 d4로 저장하였다. paste(vector1, vector2, sep=" ")을 실행하면 vector1의 모든 원소는 같은 위치에 해당하는 vector2의 원소와 공백을 경계로 하여 연결된 상태로 반환된다. POCP 매트릭스를 dist() 함수로 처리하여 euclidean distance를 구하여 이를 기반으로 계층적 군집화를 실시한다. hclust() 함수에 넘겨주는 method를 average로 지정하면 바로 우리에게 친숙한 UPGMA(Unweighted Pair Group Method with Arithmetic mean)이 된다.

출처: R Friend

실제 코드와 실행 결과는 다음과 같다. 요즘 몰두하고 있는 미생물은 Agathobaculumn, Butyricicoccus, 그리고 Eubacterium의 일부이다.

d4 = d3
k = row.names(d3)
km = gsub(" .*$", " sp. ", k)
n = 1:34
row.names(d4) = paste(km, n, sep="")
x = hclust(dist(as.matrix(d4)),method="average")
plot(x)

그런데 이렇게 그린 그림에서 세로축의 height는 우리에게 별다른 느낌을 주지 않는다. POCP 값(퍼센트)을 100에서 뺀 수치를 나만의 거리로 사용하면 안 되는 것일까? ANI도 마찬가지 성격의 지표가 된다. 이렇게 자체적으로 만든 매트릭스를 그대로 디스턴스로 삼으려면 as.distance() 함수를 쓰면 된다.

mydist = 100 - d4
x = hclust(as.dist(as.matrix(mydist)),method="average")
plot(x)
rect.hclust(x, border="red",h=50)

덴드로그램의 구조는 약간 달라진다. 그러나 여기에서는 (100 - POCP)가 height로 나타나므로, 일정 기준을 만족시키는 클러스터 주변에 사각형 경계를 씌울 수 있다. 동일 클러스터를 이루는 균주들은 POCP가 50을 넘는 것들이니 같은 genus에 속하는 것으로 볼 수 있다.

각 클러스터를 이루는 멤버들의 이름을 확인하고 싶다면? 트리 형태의 데이터를 데이터 그룹으로 잘라내려면 cutree() 함수를 쓰면 된다. cuttree가 아님에 유의하라. height = 50을 기준으로 잘라내는 코드와 결과를 같이 살펴보도록 하자. 클러스터의 수를 지정하여 잘라낼 수도 있다(k = 원하는 수).

> groups = cutree(x2,h=50)
> table(groups)
groups
 1  2  3  4 
31  1  1  1 
> names(groups[groups==1])
 [1] "Butyricicoccus sp. 1"     "Agathobaculum sp. 2"     
 [3] "Clostridia sp. 5"         "Butyricicoccus sp. 7"    
 [5] "Butyricicoccus sp. 8"     "Agathobaculum sp. 9"     
 [7] "Butyricicoccus sp. 10"    "Agathobaculum sp. 11" 

 ....(후략)

벡터의 원소를 가리키는 이름(names)이 groups 벡터를 다루는데 얼마나 요긴한지를 알 수 있다.

[2019 오디오 자작] 6N2P + 43 싱글 엔디드 앰프의 전원부 수정

전원트랜스에서 인출한 110V 교류를 다이오드 브리지로 정류하여 평활회로를 거친 뒤 초단과 출력관에 전부 연결하면 출력관의 캐소드-애노드 사이에 85V 이상을 걸기가 어렵다. 전원회로쪽의 저항을 약간 작은 것으로 하면 몇 볼트를 올릴 수 있지만 험이 커진다.

원래대로 DC-DC boost converter를 쓰면 어떨까? 6N2P 초단 회로의 B전원도 여기에서 한꺼번에 공급한 이후로는 스피커에서 간헐적으로 잡음이 들리기 시작했다. '웅`' 혹은 '윙-'의 중간쯤 되는 소리인데, 입력을 끊고 잠시 기다리면 사라진다. 이래서야 음악 감상을 하기가 곤란하다. 잠시 트랜스를 사용한 전원회로로 외도를 했던 것도 이러한 잡음 때문이었다. 6N2P 한 알이 들어가는 프리앰프만을 전원트랜스로 구동할 때에는 문제가 없었다. 물론 여기에는 많은 편법이 동원되었다. 출력이 6.3V(히터) 및 230V라서 정류를 한 뒤에는 도저히 43 오극관의 애노드에 맞출 수가 없었기 때문이다. 갖고 있는 캐피시터의 내압도 250V가 최고이고, 저항을 연결하여 낮추자니 적당한 수치의 것이 없다. 발열도 심할 것이 자명하다.

컨버터 출력에 RC 필터를 달면 어떨까? 아내의 푸념을 뒤로 하고 주말 작업에 돌입하였다. 아래 회로도에서 파랑색 곡선으로 둘러친 것이 어제 새로 추가한 필터에 해당한다. 컷오프 주파수는 15.39Hz로 계산된다. 만약 220R과 100uF 조합이라면 7.32Hz가 되어 더욱 양호한 직류 전원이 얻어질 것이다. 그러나 부품통에는 47uF 250V 캐패시터가 전부이다. 아주 초기에는 컨버터 출력에 47uF 전해 캐패시터 하나를 병렬로 연결하여 사용했던 적이 있었다.

확정된 6N2P + 43 싱글 엔디드 앰프의 회로도. 파랑색으로 표시한 RC 필터를 전원부에 추가함으로써 잡음을 해결하였다.

컨터버 모듈의 가변저항을 살며서 돌려서 출력관의 캐소드-애노드 사이에 130V 이상이 걸리게 하였다. 소리는 어떠한가? 매우 양호하고 음량도 높아졌다. 무신호 상태에서 스피커에 귀를 바싹 대고 볼륨 놉을 최대로 올리면 약간의 잡음이 들릴 수준이다. 실용적으로 문제가 없다.

만약 다음 기회에 다시 전원 트랜스를 쓰게 된다면 2차 110V 출력을 배전압 정류한 다음 저항으로 적절히 강압하여 사용하면 될 것이다.

어제 작업을 하는 동안 작은 사고가 하나 일어났다. 2016년 5월 구입한 목인두가 앰프 샤시(케익 틀)를 살짝 스치면서 불꽃이 일어나는가 싶더니 차단기가 내려가고 집안의 모든 전기가 끊겼다. 테스터로 찍어보니 인두의 히터는 단선이 되었다. 도대체 왜? 누전 때문이 아닌가 싶다. 접지되지 않은 멀티탭을 쓰다가 최근에 접지가 있는 것으로 바꾸었는데, 같은 멀티탭에는 43 앰프가 연결된 상태이다. 여기에 꽂힌 12V SMPS 어댑터는 출력의 마이너스 극과 220V 전원의 접지가 내부적으로 또한 연결이 된 상태이다. 아마도 납땜인두의 팁 부분에 누전이 일어나서 앰프의 새시 -> SMPS 어댑터의 접지를 통해 전기가 흘렀고, 차단기가 이를 감지하여 자동으로 끊긴 것으로 생각된다. 누설전류는 앰프쪽 SMPS가 아니라 인두 파워 케이블의 접지를 통해서 흘러나갔는지도 모른다. 분명히 눈 앞에서 '딱' 소리와 함께 불꽃이 튀었는데 너무나 순식간이라서 정확히 어디인지를 기억하지 못하겠다.

만약 접지가 없는 멀티탭을 계속 사용했더라면 이러한 '사고'는 일어나지 않았을지도 모른다. 그러나 이런 인두를 계속 쓰다간 부품이 망가지거나, 더욱 운이 없었더라면 내가 감전으로 다쳤을지도 모른다. 세대용 분전반에 설치된 차단기는 과전류가 흐를 때 차단을 시키는 목적도 있지만, 더욱 중요한 것은 누전을 감지하여 차단하는 것이다. 사실 이것을 이해한 것도 비교적 최근의 일이다.

남아있는 소용량 납땜인두로는 단자 등 큰 부위의 납땜이 어렵다. 안전을 위해서 다음에는 세라믹 인두팁이 장착된 것을 구입하도록 하자.

작업을 마친 상태. B전원용 컨버터는 섀시 위로 올렸다. 볼트 구멍의 흔적이 어지럽다.

LED 파일럿 램프를 달았다.

만약 이 앰프에 다시 손을 대게 된다면, 그것은 새로운 섀시에 완전히 새로 제작하기로 결심을 한 이후가 될 것이다.

[하루에 한 R] pair 형태의 데이터를 matrix로 전환할 때 주의할 점

오늘 쓰려는 주제는 R 활용 기법과 직접적인 관련은 없음을 미리 밝힌다.

한 쌍의 유전체 사이에 어떤 계산을 하여 수치를 얻었다고 가정하자. 기본적으로 얻어지는 결과물은 다음과 같이 pair 형태일 것이다.

친절한 프로그램이라면 matrix 형태로 전환한 파일을 제공해 줄 것이고, 이를 R에서 데이터 프레임으로 읽어들여서 여러 가지의 후속 분석을 실시할 수 있다. 만약 pairwise data만 위 그림과 같이 제공한다면 R에서 reshape 패키지를 이용하여 matrix 형태로 전환하면 된다.

그런데 여기에서 주의할 점이 하나 있다. 바로 자기 자신에 대한 데이터를 포함하고 있느냐를 고려해야 한다. 위에서 보인 그림 중 첫번째 것은 dRep의 결과물이다. FASTA 파일 여러개를 제공하면 프로그램이 알아서 모든 쌍에 대한 ANI를 계산하는데, 이때 자기 자신에 대한 것도 포함한다. 계산값을 보여주는 첫 줄이 그러한 것으로서 alignment coverage와 ani 전부 1,0으로 나타난다.

이것이 전부가 아니다. A->B, B->A 각 경우에 대한 값이 전부 존재하는가? 그렇다면 동일한가? 결과가 없는 것은 어떻게 할 것인가? 의외로 생각할 것이 많다. 그리고 각 경우에 대해서 처리하는 R 코드가 조금씩 달라진다.

위에 보인 그림에서 아래의 것은 POCP.sh를 사용하여 percentage of conserved protein을 계산한 것이다. 이 스크립트에서는 유전체 서열 한 쌍과 스레드 수만 인수로 제공해야 한다. 따라서 query와 subject 서열의 쌍을 별도로 계산하여 넣어 주어야 한다. 나는 Perl의 Algorithm::Combinatorics 모듈을 사용했으므로 당연히 자기 자신에 대한 계산을 하지 않게 만들었다(샘플 코드는 여기를 참조). 계산을 하나마나 POCP = 100%가 나올 것이기 때문이다. 하지만 이를 고려하지 않고 R의 reshape 패키지를 그냥 돌리면 매트릭스 형태가 좀 이상해진다. 전부 1(=100%)로 채워지는 대각선에 대한 배려가 없기 때문이다. 10개의 유전체 서열을 전부 비교하여도 R로 읽어들어 전환한 다음에는 10x10의 매트릭스에서 하나가 부족해진다. query = subject에 해당하는 대각선도 제대로 존재하지 않는다.

이 문제를 해결하려면 pairwise 결과 파일의 뒷부분에 자기 자신에 대한 분석값을 넣는 것이 좋다. A, B, C,...J의 10개 유전체에 대해서 POCP를 계산한다면 다음의 10줄을 넣기만 하면 된다. 자기 자신에 대한 POCP 값은 당연히 100%이므로 이를 세번째 컬럼에 입력해 넣는다.

A A 100.0
B B 100.0
C C 100.0
..
J J 100.0

이렇게 하여 reshape의 cast() 함수를 돌려 처리를 하면 10x10의 매트릭스, 정확히 말하자면 데이터 프레임(df라 하자)이 생길 것이다. A->B, B->A 두 경우에 대해서 한쪽으로만 계산을 하였으니 데이터 프레임으로 전환된 뒤의 구조를 확인해 보면 대각선은 전부 100이고 나머지 절만만 채워진 상태일 것이다.

NA를 전부 0으로 치환한 다음 df와 t(df)를 더하면 대각선 절반에 대한 값이 채워진다. 그러나 대각선은 전부 200이 될 것이다. 이는 또 적절히 처리하면 된다.

오늘의 교훈은 무엇인가? 데이터 프레임을 한번쯤은 엑셀 등에서 열어보고 어떤 상태인지를 눈으로 확인하는 것이 좋다.

2019년 3월 14일 업데이트

Algorithm::Combinatorics의 combination_with_repetition() 서브루틴을 사용하면 자기 자신에 대한 검색 결과도 만들 수 있다. 이것을 사용하여 POPC.sh를 실행하는 Perl 코드는 다음과 같다.

#!/usr/bin/perl

use strict; use warnings;
use Algorithm::Combinatorics qw(combinations_with_repetition);

open LIST, 'list';
my @temp = ();
while (<LIST>) {
    chomp;
    push @temp, $_;
}
my $faa_ref = \@temp;

my $iter = combinations_with_repetition($faa_ref, 2);

while (my $c = $iter->next) {
    my $out = 'pocp' . '_' . $c->[0] . '_vs_' . $c->[1];
    print "Running for $out...\n";
    system("POCP.sh $c->[0] $c->[1] 24 > $out");
}

[2019 오디오 자작] 6N2P 프리앰프와 43 전력증폭회로의 통합

오늘의 제목에서 '통합'이라 함은 두 개의 스테이지(초단 + 전력증폭단)를 하나의 섀시에 올리고 동일한 전원을 공급하게 만들었다는 뜻이다. 진공관 앰플리파이어에서 '초단'은 여러 의미를 동시에 갖는다. 아주 제대로 만들어진 장비라면 전치증폭기(프리앰플리파이어), 드라이브 회로, 위상반전회로(푸시풀 앰프를 위해)를 통틀어서 초단이라고 부를 수 있다. 각 세부적인 기능은 별도의 진공관이 해결하는 것이 일반적이나 내가 만드는 앰프는 단 하나의 3극관이 초단의 일을 맡는다. 푸시풀 앰프가 아니므로 위상반전단은 존재하지 않는다.

43이라는 오극관의 작동 전압이 워낙 낮으므로 초단관도 이에 맞추어야 했고, 더 큰 문제는 두 진공관의 히터 전압이 매우 다르다는 것. 43는 25V, 6N2P는 6.3V가 필요하다. 대충 24V를 연결해도 43 오극관 작동에는 문제가 없다. 그렇다면 요즘 널리 쓰이는 스텝다운 컨버터를 쓰면 24V에서 6.3V를 만들 수 있지 않겠는가?

두 개의 SMPS 어댑터를 이용한 6N2P+43 진공관 싱글 앰프의 전원부 구성. 어댑터의 용량은 각각 12V 5A, 24V 2A이다. 이론적으로는 용량이 충분한 24V SMPS 하나만 있으면 이를 승압하여 B+ 전원도 공급 가능하다. 그러나 내가 구한 24V 고용량 SMPS는 이렇게 연결했더니 잡음이 심했다. SMPS의 문제인지 컨버터의 문제인지는 잘 모르겠다.

LM2596칩을 사용한 이 컨버터 모듈은 최대 1.5A의 전류를 흘릴 수 있으며 IC114에서 3500원에 팔린다(링크). LED 표시창에서는 입력과 출력 전압을 전부 표시하므로(버튼 스위치를 눌러 전환) 매우 쓰기에 편리하다. 어제의 퇴근 후 프로젝트는 이 모듈을 사용하여 앰프를 다시 완성하는 것이었다.

얼추 완성이 된 이 앰프는 사실 대단히 위험하다. 왜냐하면 만능기판을 뒤집어서 만든 6N2P 드라이브 회로에는 위에 노출된 배선을 따라 155V의 직류가 흐르기 때문이다. 일반적인 진공관 앰프에 비해서는 상당히 낮은 편이지만 그래도 조심해야 된다!

너무 가공을 많이 하는 바람에 케익틀은 잘못 뚫은 구멍들로 엉망이 되었다. 소리는 별다른 잡음 없이 잘 난다. 12V 어댑터를 꽂는 소켓의 내경이 잘 안맞는지 케이블을 당기면 접촉 불량이 약간 일어난다. 흔히들 간과하는 것이 전원 어댑터 잭의 내경을 맞추는 일이다. 가장 많이 쓰이는 어댑터에서는 내경 2.1mm와 2.5mm의 두 가지 규격이 있으니 주의가 필요하다.

전통적인 전원장치를 쓰고 싶다면 다음 회로도를 쓰면 된다. 반도체(다이오드)를 정류기로 채용하였으니 완벽히 전통적인 장치가 아니라고 볼 수도 있다. ~55V에 43 5극관 두 개의 히터를 직렬로 연결하는 아이디어는 아직 구현해 보지 않았다. 다만 B+ 전원을 이것으로 연결하면 스피커에서 약간의 험이 들리는 것은 사실이다. 초크 코일을 쓰면 나아지기는 하겠지만.

진공관 히터에 전원을 공급하여 가열하는 것을 '점화'한다는 표현을 많이 쓰게 된다. 아마도 일본쪽 기술 서적에서 쓰던 낱말을 그대로 가져다 쓰는 것이 아닐까 한다. 백열전구를 켤 때, 점등을 한다는 말은 쓰지만 우리나라에서는 점화한다고는 쓰지 않는다.

왜 점등·점화에는 한자 '점(點)'을 쓰는 것일까? 심지에 점 찍듯이 불을 붙인다고 해서 이런 한자어가 생겨났다고 한다(링크).

[BASH] cut 명령을 이용하여 파일 이름을 간략하게 정리하기

cut - remove sections from each line of files (print selected parts of lines from each file to standard out)

커맨드 라인에서 man cut이라고 치면 나오는 설명이다. CSV 파일처럼 각 라인이 일정한 구분자(delimiter)로 나뉜 'section'의 집합으로 이루어진 경우, 이를 나누어서 작업하는 데에 아주 적당한 명령이다. 구분자를 제공하지 않으면 라인 내의 위치를 기준으로 잘라내는 것도 가능하다.

cut을 사용하면 NCBI에서 다운로드한 유전체 정보 파일의 이름을 간략히 하는 일이 매우 쉬워진다. RefSeq에서 받은 FASTA 파일의 이름은 다음과 같이 여러 필드가 밑줄('_')로 구분된 형태를 하고 있다.

GCF_000152245.2_ASM15224v2_genomic.fna
GCF_000242675.1_Euba_infi_F0142_V1_genomic.fna
GCF_900100105.1_IMG-taxon_2593339210_annotated_assembly_genomic.fna

굵은 글씨로 표현된 assembly accession 뒤의 정보는 assembly ID로서 제출자가 정하는 것이라서 특별히 정해진 규약이 없다. 그래서 rename 명령어나 $VAR{//문자열} 방법으로 처리하기에 좋지 않다. cut을 사용하면 밑줄을 delimiter로 삼아서 첫번째 및 두번째 섹션을 취한 뒤, 여기에 .fna 확장자를 붙이는 것으로 새 파일의 이름을 정하면 된다. cut을 사용하면 이러한 작업을 매우 간단하게 할 수 있다.

$ ls | while read f
> do
> cp $f $(echo $f | cut -d'_' -f 1,2).fna
> done

더 이상의 무슨 설명이 필요한가? 더욱 많은 사례는 Linux Cut Command With Samples를 방문해 보라.

cut 명령어의 사례를 보면 -f1,2 -f3과 같이 옵션과 지정한 값 사이에 아무런 공백을 넣지 않는 것을 보게 된다. 옵션을 길게 쓰려면 --fields=LIST와 같이 등호('=') 좌우로 공백을 넣지 않는 것이 맞다. 그러나 짧게 쓰는 경우라면 -f 1처럼 사이에 공백을 넣는 것이 자연스럽지 않을까? 실제로 작동을 시켜 보면 -f1과 -f 1 어느 것으로 하든 상관은 전혀 없다. 어쩌면 이러한 관행은 옵션을 인수(여기에서는 작업 대상 파일명)와 명확히 구분하기 위함인지도 모른다. 옵션과 인수(argument)는 다름에 유의해아 한다. 옵션은 명령의 실제 동작을 세부적으로 조절하는 것이고, 인수는 작업이 이루어지는 대상이다.

POCP(percentage of conserved proteins): 세균의 속(genus) 구분을 위한 기준

두 미생물 균주가 같은 종(species)에 속하는지를 확인하는 지표로서 아직까지도 가장 널리 쓰이는 것은 16S rRNA 유전자의 서열 유사도(sequence similarity)이다. 두 균주의 16S rRNA gene sequence similarity가 97% 미만이면 이는 서로 다른 종에 속하고(역은 항상 성립하는 것은 아님), 95% 미만이면 서로 다른 속에 속한다는 것이다(Stackebrandt, E. & Goebel, B. M. (1994). Taxonomic note: a place for DNA-DNA reassociation and 16S rRNA sequence analysis in the present species definition in bacteriology. Int J Syst Bacteriol 44, 846– 849). 종 구별을 위한 서열 유사도 cut-off는 그 이후에 98.7%까지 올라갔고(Stackebrandt, E. & Ebers, J. (2006). Taxonomic parameters revisited: tarnished gold standards. Microbiol Today 33, 152–155), 2014년에는 98.65%가 기준값으로 제시되었다.

그러나 2015년 IJSEM에 실린 논문 "Cautionary tale of using 16S rRNA gene sequence similarity values in identification of human-associated bacterial species"에 의하면 inter-species sequence similarity 범위 95 & 98.7% threshold를 준수(?)하는 실제 사례는 그렇게 많지 않다고 한다. 이 논문의 저자들은 따라서 genus마다 고유하게 적용할 수 있는 sequence similarity의 범위를 사용할 것을 제안하였다.

유전체의 시대를 맞은 요즘 세균의 종 구분 지표로써 가장 널리 쓰이는 것은 두말할 나위 없이 ANI(average nucleotide identity)이다. 그러나 이 값은 속(genus)를 판별하기에는 적합하지 않다. Qin 등은 2014년에 50%의 POCP(percentage of conserved proteins)가 속 구분의 지표가 됨을 제안하였다.

J Bacteriol. 2014 Jun;196(12):2210-5. doi: 10.1128/JB.01688-14. Epub 2014 Apr 4.
A proposed genus boundary for the prokaryotes based on genomic insights. PMID 24706738

그림 출처: 링크. 같거나 다른 속 사이에서 16S rRNA gene identity나 ANI의 분포는 상당 부분 서로 겹친다.

그림 출처: 링크.

The percentage of conserved proteins (POCP) between two genomes was calculated as [(C1 + C2)/(T1 + T2)] ·  100%, where C1 and C2 represent the conserved number of proteins in the two genomes being compared, respectively, and T1 and T2 represent the total number of proteins in the two genomes being compared, respectively.

이 논문에서는 두 단백질을 서로 BLASTP로 검색했을 때 다음의 조건을 만족시키면 conserved protein인 것으로 간주하였다.

• an E value of less than 1e−5
• a sequence identity of more than 40%
• an alignable region of the query protein sequence of more than 50%.

그러나 이 논문에서는 POCP 계산에 필요한 스크립트를 제공하지 않는다. Harris 등은 이에 따라 "Phylogenomics and comparative genomics of Lactobacillus salivarius, a mammalian gut commensal (링크)" 논문에서 POCP 매트릭스 산출에 사용했던 스크립트를 공개하였다. 실제의 스크립트는 figshare 사이트에 있다(DOI: dx.doi.org/10.6084/m9.figshare.4577953.v1). 봉사정신이 투철한 개발자가 이를 잘 손을 보아서 bioconda나 GitHub 같은 곳에서 배포하면 좋을 것이다.

[2019 오디오 자작] 6N2P 프리앰프 실험-III, 돌고 돌아 다시 43번 오극관에게...

이번에 만든 6N2P 프리앰프가 반도체 앰프와 연결하여 좋은 성능을 냈으니 43번 오극관을 사용한 전력증폭회로에 붙여서 쓰지 못할 이유가 없다. 그 동안 43번 앰프 회로는 OP amp를 이용한 드라이브단을 사용하고 있었다. 어제의 실험 회로는 다음 그림과 같다. 반도체를 전부 몰아내고 전부 진공관으로만 구성된 오디오 앰프를 추구하는 사람도 있다. 이 또한 확고한 철학을 기반으로 하는 목표이지만 전원부에서는 반도체 소자를 쓰지 않을 도리가 없다. 적어도 내 기술 수준으로는 그렇다.

이 회로도에서는 초단관(6N2P)와 출력관(43)에 각각 별도의 전원을 쓰는 것으로 그렸지만 꼭 그럴 필요는 없다. SMPS + DC-DC boost converter로 두 종류의 진공관에 B+를 전부 공급하니 가장 깨끗한 소리가 났다. 그 반대로 트랜스와 정류/평활회로를 사용하여 공급하면 볼륨을 최대로 할 경우 약간의 험이 들린다.

세상에 이렇게 엉성한 앰프는 없을 것이다. 기성품이 아니라 자작이라는 것이 더 이상 변명이 되지는 못한다. 그렇지만 내 귀에 듣기에 소리는 만족할만한 수준이다. 음량 조절용 가변저항을 드라이브단과 출력단 사이에 넣은 것이 좋은 선택이었던 것으로 생각된다. 43번 오극관의 바로 앞에 가변저항이 위치하여 그리드와 접지를 연결하고 있으므로 별도의 그리드 리크 저항은 필요하지 않다.

프리앰프부와 전원부.

43번 오극관 전력증폭회로를 연결하였다. 소리는 이만하면 OK.

6N2P 프리앰프에 별도의 섀시에 수납하지 않고 이 상태 그대로 유지하면서 실험용으로만 사용할 것인가? 혹은 43번 오극관 출력단과 하나의 앰프로 구성할 것인가? 만약 그러하다면 어떤 섀시를 만들까? 나무 도마 위에 늘어놓을까? 케이크 틀에서 도마에 이르기까지 주방와 관련된 물품보다 더 좋은 아이디어는 떠오르지 않는다.

[하루에 한 R] 리스트 생기초

data = c("a", "b", "c", "d", "e")라는 벡터가 있다고 가정하자. list(data)와 list("a", "b", "c", "d", "e")는 전혀 다른 구조의 리스트를 만들어 낸다.

data = c("a", "b", "c", "d", "e")
d1 = list(data)
d2 = as.list(data)
d3 = list("a", "b", "c", "d", "e")

as.list() 함수를 쓰지 않으면 인수로 주어진 벡터를 하나의 원소로 갖는 리스트를 만든다. 벡터의 개별 항목이 리스트의 원소가 되게 만들려면 as.list() 함수를 쓰거나, 혹은 list() 함수의 인수를 개별적으로 전부 나열해야 한다.

d2와 d3 리스트는 같은 동일한 구조이다.

d2(=d3) 리스트의 원소에는 아직 이름이 지어지지 않았다. 이름 정보를 수록한 벡터를 nato = c("Alpha","Bravo","Charlie","Delta","Echo")라 하면 names() 함수를 사용하여 이를 리스트의 각 원소에 부여하면 된다.

리스트를 생성하는 단계에서 이름을 부여할 수도 있다. setNames() 함수가 이러한 때에 쓰인다. 첫번째 인자가 벡터가 아니라 리스트임을 명시적으로 보여야 한다.

만약 실수로 as.list() 함수를 쓰지 않았다면 어떻게 될까? 그러면 단순한 named vector가 된다.

edit() 함수를 이용하여 실제로 이들 데이터 구조가 어떻게 선언되는지를 알아보자. 

만약 리스트의 원소 중에 벡터가 포함되어 있다면 모양이 좀 더 복잡할 것이다.

[2019 오디오 자작] 6N2P 프리앰프 실험-II

작년에 만들어 두었던 프리앰프 보드의 수정에 착수하였다. 입력측에 위치하던 음량 조절 전위차계를 출력쪽으로 옮기고 대신 grid leak 및 grid stopper resistor를 붙이는 것이 가장 큰 수정사항이었다. 출력측에 연결한 반도체 앰프의 전단에 이미 전위차계가 붙어 있어서 오늘의 실험에서는 커플링 캐패시터를 반도체 앰프의 입력에 직결하였다.

오늘의 실험에 사용한 6N2P 프리앰프 회로. 전원장치 회로도는 여기에 있다.

위에서 보인 '발로 그린 회로도'는 계산이나 검증을 거친 것이 전혀 아니다. 인터넷에 돌아다니는 심플한 프리앰프 회로를 적당히 짜깁기한 것에 불과하다. 다만 비교적 낮은 애노드 전압을 쓰는 회로를 좀 더 참고하였다. Triode / Pentode Loadline Simulator 사이트에서 힌트를 얻을 수 있을 것이다. 그리드 바이어스 전압의 측정치는 -0.65V 정도였다.

만능기판의 배선 모습을 사진으로 남겼다. 자세히 살펴보면 히터 점화용 전선의 연결 상태가 그다지 좋지 않다. 와이어 스트리퍼로 전선(연선)의 피복을 벗길 때 가장 짜증스러운 것은 심선 몇 가닥이 이와 같이 끊어질 때이다. 땜납으로 덧칠을 좀 더 해야 되겠다. 히팅 건이 없어서 라이터로 수축튜브를 가열해야 하니 항상 검게 그을음이 묻는 것도 마음에 들지 않는다. 해결 방법이 있기는 하다. 검정 수축튜브를 쓰면 된다...

페놀만능기판에 배선하는 한 이보다 더 간단하게 할 방법이 있을까? 회로 자체가 워낙 간단하기 때문에 하드 와이어링이든 고급 에폭시 기판이든 이런 수수한 배선과 비교하여도 음질에는 차이가 없을 것이라고 생각하였다.

반도체 앰프에 연결하여 소리를 들어 보았다. 반도체 앰프 자체의 게인이 워낙 높아서 음량 조절용 전위차계의 놉을 조금만 돌려도 매우 큰 소리가 났다. 전위차계가 프리앰프의 전단에 있을 때보다는 음량을 조절하기가 훨씬 수월하였고, 수정하기 전의 회로와 비교하면 훨씬 양호한 작동 상태를 보여 주었다. 히터 전원의 가상 중점을 접지에 연결했을 때 잡음이 줄어드는 효과는 명백하였다.

오늘 실험의 결과에 점수를 매긴다면 83점 정도? 만약 낙제점이 나온다면 다 거두어 넣고 앞으로 모든 오디오 자작 활동을 접을 생각까지 하던 차였다. '이번이 마지막'이라는 비장한 각오로 일을 벌이면 이처럼 실낱 같은 희망이 보인다는 것이 문제이다. 그래서 아직 유치한 수준의 이 취미를 손에서 놓지 못하는 것이다.

[하루에 한 R] read.table() 함수는 왜 #로 시작하는 라인을 읽지 못하는가?

데이터를 텍스트 형태로 수록한 자료 파일에서 주석(comment) 라인은 보통 첫번째 글자를 '#'으로 시작하는 경우가 많다. Shell 혹은 Perl script에서 주석을 달 때 쓰는 관행을 그대로 답습하는 것으로 보인다.

자료 파일의 중간에 주석을 넣을 수도 있지만, CSV(comma-separated values) 파일에서는 첫 줄에 컬럼의 이름에 해당하는 정보를 담기도 한다. antiSMASH database에서 다운로드한 파일의 사례를 보자. 확장자는 csv지만 실제 구분자는 탭이다. 첫 줄은 컬럼의 이름 정보를 싣는 동시에 '#'로 시작함으로서 주석임을 명시하였다.

이것을 read.table("file.csv",sep="\t",header=T)로 읽으면 어떻게 될까? '#Genus', 'Species', 'Strain'...이 컬럼 이름으로 들어가는 일은 절대로 벌어지지 않는다. 첫줄이 '#'로 시작하면 comment로 인식하여 읽어들이는 과정에서 그냥 지나가기 때문이다.

이를 제대로 읽으려면 다음과 같은 명령을 써야 한다.

data = read.table("file.csv",sep="\t",header=T,comment.char="")

read.table() 함수의 기본 동작은 comment.char="#"이다. comment.char에 지정할 수 있는 값은 당연히 단일 문자여야 한다. ""를 지정하면 모든 라인을 주석으로 해석하지 않는다. 읽어들인 파일이 데이터프레임으로 잘 전환되었는지를 확인해 보자.

왜 첫번째 컬럼의 이름이 X.Genus인가? 원본 파일의 값은 #Genus인데 특수문자가 있어서 이를 그대로 컬럼 이름으로 쓰기에 곤란하기 때문이다. 원본 파일의 #Genus를 AGenus, 1Genus, _Genus, 등으로 바꾸어 보면서 R이 이를 어떻게 적절히 컬럼 이름으로 변환하는지를 확인해 보라.

read.csv() 함수는 comment.char=""임을 기억해 두자. 주석이 없는 순수한 데이터값으로만 이루어진 파일을 읽어들이는 상황이 대부분일 것이라고 가정하여 동작한다는 의미가 아니겠는가?

함수의 용법을 정확히 하는 것 못지않게 데이터 파일을 텍스트 편집기로 한번쯤은 열어서 훑어보면서 헤더는 어떻게 구성되어 있는지, 중간에 주석줄은 없는지 등을 확인해 보는 것이 현명할 것이다.

[2019 오디오 자작] 6N2P 프리앰프 실험

모든 전자장치 DIY의 기본은 전원회로를 잘 만드는 것에서 시작한다. 1N4007(위키피디아) 네 개를 평면 상의 사각형 모양으로 엮어서 만든 정류용 브리지를 어떻게 매만져야 보기에도 좋고 공간도 적게 차지할지 고민한 끝에 다음과 같은 최종 결과물을 얻었다. 왼쪽이 교류 입력, 오른쪽이 정류된 직류(정확히 말하자면 맥류) 출력에 해당한다. 직육면체를 구성하는 평행한 네 개의 모서리에 다이오드가 하나씩 위치하게 만든 것이다.

만능기판에 브리지 정류회로를 구성해도 되지만 생각을 게을리하면 예쁜 배선 패턴이 나오질 않는다. 잠시 검색을 통하여 모범이 될만한 패턴을 찾아보았다. 다음의 것이 아주 적당하다.

출처: http://hkpinvent.blogspot.com/2013/

알리익스프레스에서 구입한 뒤 아직 만들지 않은 정류회로 보드의 패턴도 옮겨서 그려 보았다. 위상수학의 관점에서는 위의 것이나 아래의 것이나 똑같은 그림일 것이다.

2019년 설 연휴의 끝자락에 이루어진 이번 DIY에서는 6N2P를 사용하여 만들었던 프리앰프부 회로가 과연 쓸만한 물건인지를 알아보자는 것이 목표였다. 전원부는 보유한 부품을 최대한 활용하기로 하였다. 43 pentode 싱글 앰프에서 12AU7을 꽂아서 드라이브 회로로 썼다가 갑자기(?) 발생한 극심한 잡음으로 인하여 퇴출이 된 상태였다. 12AU7은 그 후로 몇 번 점검을 했지만 상태가 좀 이상한 것 같아서 갖고 있던 4개의 6N2P 중 하나를 쓰기로 했다.

6N2P를 사용한 프리앰프 회로 중 가장 모범적(?)인 것은 다음과 같다.

출처: How to use 6N1P and 6N2P Russian tube

이 프로젝트를 구상하면서 처음에 참조한 회로도는 Matt Renaud의 4S universal tube preamplifier for 12A*7 tubes이다. 여기에서는 캐소드 바이패스 캐패시터가 선택으로 되어 있고 음량 조절용 폿은 출력부에 위치한다. 위에서 보인 6N2P 회로와 비교하면 바이어스 저항의 값이 가장 많이 다르다.

최근에 참조한 회로도는 다음과 같다("an example of a 20 dB audio gain stage"). B+ 전압은 50-150V를 권장하였다.

출처: Hackaday

6N2P를 비교적 낮은 B+ 전압으로 구동하려면 나머지 부품의 수치는 어떻게 해야 하는가? 6N2P tube pre-tone low volt & low cost라는 글을 참조하면 약간의 아이디어를 얻을 수 있다. 출력측은 커플링 캐패시터에서 그대로 뽑아내는 것이 아니라 몇 가지의 수동 부품이 여러개 엮여진 다소 복잡한 형태이다. 뭔가 이유가 있을 것이다.

출처: 6N2P tube pre-tone low volt & low cost

그러나 부품을 교체하는 일은 하지 않았다. 계산과 실험에 의해서 최적의 부품값을 결정해야 하는데, 너무 주먹구구로 일을 하고 있다! B+ 전압, 더욱 정확히 말하자면 anode voltage(Va or Vak)가 낮아지면 그리드 바이어스 전압을 높여야(0에 가까와지도록) 하므로 당연히 캐소드 저항치도 작아져야 할 것이다(맞나??).

전원회로를 직접 만들어야 하는데 마침 내가 갖고 있는 캐패시터 중 내압이 가장 높은 것은 250V가 고작이다. 어떻게 할 것인가? 전원트랜스(1차 220V, 2차 230V + 6.3V, 링크)의 2차에 강압용 단권트랜스(1차 440V-380V, 2차 220V-110V, 링크의 사진에서 오른쪽)을 뒤집어 연결하여 115V 정도를 얻어내는 것으로 결정하였다(아래 회로도). 다이오드 브리지의 모습은 이 글의 도입부에서 소개하였다. 저항과 캐패시터는 하드와이어링으로 대충 연결하였고, 최종단의 0.1uF 캐패시터는 다른 회로에서 떼어낸 Dain MPX 275VAC(metallized polypropylene film capacitor)를 재활용하였다. 갖고 있던 세라믹 캐패시터는 내압이 낮아서 쓸 수가 없었다.

대충 만든 앰프에서는 대충 소리가 난다.. 아직은 그렇게 만족스러운 수준의 소리가 아니다. 앞으로 이 실험용 프리앰프가 삼극관을 이용한 전압증폭회로(기본 중의 기본!)를 공부하는데 좋은 교재가 될 것으로 믿는다.

[하루에 한 R] Key/value 형태의 자료를 이용한 치환 - named vector 이야기

Perl에서는 hash라는 자료형을 이용하여 key/value 형태의 자료를 다룬다.

$color_of{apple} = 'red';
# initialize a hash
my %color_of = (
    "apple"  => "red",
    "orange" => "orange",
    "grape"  => "purple",
);

R에서는 이와 유사한 자료 구조로서 list가 있지만 hash보다는 훨씬 복잡한 구조의 데이터를 담을 수 있다. 예를 들어서 list_data = list("red","green",c(1,2,4), TRUE,13.43)의 사례에서 볼 수 있듯이 문자열, 숫자, 벡터, 논리값 등 어떤 타입의 원소든지 담을 수 있다. 물론 Perl의 hash에서도 참조(reference)를 이용하면 어떤 자료형이든 scalar로 전환되므로 hash의 값으로 저장할 수 있다. Perl의 hash에서는 각 원소에 접근하려면 key를 사용해야 되지만, list는 훨씬 다양한 방법으로 인덱싱하여 접근하거나 혹은 슬라이싱이 가능하다.

이번 글에서는 R에서 key/value 목록을 사용하여 간단하게 데이터를 변환하는 방법을 알아보고자 한다. 알파벳 소문자 a, b, c, d, e를 NATO 음성 문자(NATO phonetic alphabet) alpha, bravo, charlie, deltam echo로 변환하는 아래의 사례에서는 list를 쓰지 않고도 named vector를 사용하고 있다. Perl의 사고 체계와는 약간 다르다는 것을 이해해야 한다. Perl에서는 $hash{key} = value의 형태로 자료를 정의하므로 특정 key에 해당하는 value를 확인하는 과정이 복잡하지 않다. 그러나 R의 named vector에서는 names(var)[var=="key"] 구문을 사용해야 key에 대응하는 value가 반환된다. 언뜻 난해해 보인다.

> abcde = c("a","b","c","d","e")
> nato.abcde = c("Alpha","Bravo","Charlie","Delta","Echo")
> names(abcde) = nato.abcde
> abcde
  Alpha   Bravo Charlie   Delta    Echo 
    "a"     "b"     "c"     "d"     "e" 
> names(abcde)[abcde=="a"]
[1] "Alpha"
> data = c("e","e","a","b","c","d")
> for(i in data){print(names(abcde)[abcde==i])}
[1] "Echo"
[1] "Echo"
[1] "Alpha"
[1] "Bravo"
[1] "Charlie"
[1] "Delta"

내용상으로 value에 해당하는 것이 names() 함수의 값으로 반환되는 것이 좀 어색하다. 왠지 names()로 지정 또는 반환되는 값은 key에 해당하는 것이 더 자연스러울 것이라는 생각이 든다. 다시 말해서 "a"를 입력했을 때 "Alpha"를 얻고자 한다면, 다음 중 어느 것이 더 자연스럽겠는가? 두 경우에서 myVar 벡터의 구조는 서로 정반대이다.

1. result = myVar["a"] 
2. result = names(myVar)[myVar=="a"]

두말할 나위도 없이 (1)번이 더 직관적이다. 하지만 이 글의 윗부분에서 든 사례에서는 (2)의 방식을 택하고 있다. 어떤 방식이 더 낫다고는 할 수 없다. 다만 어떤 벡터 원소의 이름은 그 원소의 속성 중 하나라고 이해하면 (2)의 방법도 수긍이 간다. 위에 보인 R 코드에서 abcde를 입력했을 때 화면으로 출력되는 자료의 모습을 한번 음미해 보자. Alpha, Bravo, Charlie...는 마치 데이터프레임의 column label과 흡사하다. 이러니 혼동을 불러 일으키는 것도 당연하다.

좀 더 현명한 방법은 아래와 같이 list를 사용하는 것이다. setNames() 함수를 사용하면 리스트를 생성하면서 동시에 이름을 부여할 수 있다(괄호 안에서 as.list() 함수를 쓰지 않으면 단순한 named vector가 된다). 여기에서는 names()가 가리키는 방향이 위의 named vector의 사례와는 반대임을 유의해야 한다. 즉, key에 해당하는 것을 가리킨다.

> myList = setNames(as.list(nato.abcde),abcde)

# myList = list(a="Alpha",b="Bravo",c="Charlie",d="Delta",e="Echo")
# add a new key/value element
> myList$f = "Foxtrot"
> myList[["h"]] = "golf"
> names(myList)[names(myList)=="h"] = "g"
> for (i in ls(myList)) {print(myList[[i]])}
[1] "Alpha"
[1] "Bravo"
[1] "Charlie"
[1] "Delta"
[1] "Echo"
[1] "Foxtrot"
[1] "golf"
> for (i in data) {print(myList[[i]])}
[1] "Echo"
[1] "Echo"
[1] "Alpha"
[1] "Bravo"
[1] "Charlie"
[1] "Delta"

'펄(Perl)'과 '알(R)'의 짤막한 비교언어학이었다.

[업데이트 2019-02-18] 실용적인 사례

위에서 보인 예문에서는 for 반복문 안에서 단지 print를 이용하여 값을 화면에 인쇄한 것에 지나지 않는다. 실제로는 data = c("e","e","a","b","c","d")라는 벡터가 주어졌을 때 각 원소를 이에 대응하는 "Echo"    "Echo"    "Alpha"   "Bravo"   "Charlie" "Delta"로 바꾸어 벡터로 반환해 주어야 비로소 쓸모가 있다. Named vector를 사용하여 간단하게 처리하는 코드를 소개해 본다. Named vector 'abcde'는 위에서 정의한 것을 그대로 사용한다. 변환 작업을 할 벡터 data가 주어졌을 때 여기에 원소와 동일한 이름을 붙이는 것이 핵심이다.

> abcde = c("a","b","c","d","e")

> nato.abcde = c("Alpha","Bravo","Charlie","Delta","Echo")
> names(abcde) = nato.abcde

> data = c("e","e","a","b","c","d")
> names(data) = data
> for (i in data) {
+ temp = names(abcde)[abcde==i]
+ names(data)[names(data)==i] = temp
+ }
> names(data)
[1] "Echo"    "Echo"    "Alpha"   "Bravo"   "Charlie" "Delta"  

위의 방법에서는 data vector의 이름을 변경하였다. 이름과 값이 동일하다면, 값을 바꾸어도 된다. 맨 마지막 행에서는 unname() 함수를 사용하여 data 벡터의 이름을 제거해야 원하는 결과가 나온다.

> data = c("e","e","a","b","c","d")
> names(data) = data
> for (i in data) {
+ temp = names(abcde)[abcde==i]
+ data[names(data)==i] = temp
+ }
> unname(data)
[1] "Echo"    "Echo"    "Alpha"   "Bravo"   "Charlie" "Delta"  

named vector에서는 이름이 같아도 문제가 되지 않는다. 아마 data frame에서는 곤란할 것이다.