2022년 2월 27일 일요일

43 오극관 앰프의 전원 트랜스 교체

돌고 돌아 결국은 원래의 모습으로 돌아갔다. 50VA급 220V:220V 트랜스를 써서 얻은 전압을 정류하면 43 오극관에 그대로 공급하기에는 너무 높다. 몇 개의 저항을 거쳐서 160V 정도로 낮추었더니 전원 트랜스에서 듣기 싫은 떨림이 발생하였다. 똑같은 전원 트랜스를 6LQ8 싱글 및 푸시풀 앰프에서는 전혀 문제 없이 쓰고 있다.

아세아 전원의 50VA 220V:220V 복권 트랜스. 3개째 구입한 것이다. 다음에 6LQ8 PP amp를 하나 더 만들게 되면 그때 사용하자. 어쩌면 6V6 싱글 앰프에 응용할 수도 있겠다.

 그래서 처음 43 싱글 앰프를 만들었던 때의 조합으로 되돌아가기로 하였다. 220V:230V + 6V 전원 트랜스에 440V + 380V:220V + 110V 단권형 트랜스(30VA)를 조합하는 것. 단권형 트랜스에서는 밑줄친 탭을 사용하였다. 이 무슨 해괴한 조합이란 말인가? 사실 220V:110V 트랜스 하나면 얼추 해결될 일이었겠지만, 그러려면 또 비용이 들게 되므로 갖고 있는 트랜스 중에서 110V를 조금 넘는 출력을 얻을 수 있는 조합을 만들어 본 것이다. 

MOSFET을 사용한 리플 제거 회로에서 12V 정도를 잡아먹기 때문에, 2차에 150V가 나오는 전원 트랜스가 있었다면 초단관(6N2P)과 출력관(43) 전부를 그런대로 만족시키는 B+ 전압을 뽑을 수 있었을 것이다. 그러나 이런 트랜스는 주문제작을 하지 않으면 살 수가 없다. 진공관 오디오용으로 만들어진 표준형 전원 트랜스는 가격도 비싸고, 이렇게 낮은 전압을 출력하지는 않는다.

1940년대에 만들어진 라디오 수신기용 출력관을 응용하여 오디오 앰프를 만드느라 나도 나름대로 고생을 많이 하였다. 변변한 회로도가 없는 상태에서 출발하여 CAD로 만든 상판을 써서 제작을 하였으니 꽤 많은 공부와 경험을 하였노라고 자부해도 좋을 것이다.

6LQ8 싱글 앰프 RCA 단자의 불량을 해결하는 것을 포함하여 2022년 2월의 오디오 DIY를 마치도록 한다. 헐거은 RCA 단자는 교체 말고는 답이 없었다. 부품은 좋은 것을 쓰자! 특히 단자류의 품질에 소홀해서는 안 된다. 직접 손으로 그 품질을 느낄 수 있고, 결과는 소리에 고스란히 반영되기 때문이다.


2022년 2월 26일 토요일

Linux용 Windows 하위 시스템을 쓰게 되다

올해에는 입문자를 위한 미생물 유전체 분석 강좌를 개설하기로 하였다. 마지막으로 강사 역할을 했던 것이 벌써 2019년 11월이니 3년이라는 세월이 흘렀다. 이런 강좌를 준비할 때마다 항상 실습 환경을 마련하는 것이 고민이다. 수강생들이 리눅스 경험을 갖고 있을까? 리눅스 환경에 어느 정도는 익숙하다는 것을 전제로 준비를 해야 할까, 아니면 아예 없다고 가정하고 준비를 해야 할까? 리눅스에 대한 흥미를 유발하게 한다면 더할 나위가 없겠지만, 시커먼 화면에 글씨만 보이는 터미널에 대해서 거부감을 갖는 사람들이 적지 않다는 것이 현실이다.

2019년 강좌에서는 KOBIC에서 마련해 준 서버(우분투)에 프로그램과 실습용 데이터를 미리 설치한 다음 수업을 진행했었다. 그러나 교육 기간이 끝나면 서버를 닫아야 하니 복습 또는 심화 학습을 하기기 어려웠다. 나는 이를 대비하여 VirtualBox ova 파일과 스크립트를 별도로 배포하여 교육 기간이 끝난 뒤에도 직접 분석 과정을 경험해 보도록 나름대로 배려를 하였다.

올해 계획하고 있는 교육은 리눅스 경험이 없는 사람을 대상으로 할 것이다. 수강생들은 각자 노트북 컴퓨터를 갖고 오게 될 것이고, 윈도우 + 인터넷 환경에서 미생물 유전체 분석 과정에 대한 맛보기를 경험하게 될 것이다. 사실 아무리 초급 수준의 생명정보학이라 해도 리눅스를 배제하고 할 수 있는 일은 그렇게 많지 않다. 교육을 마치고 돌아간 뒤, 하드웨어 사양이 좋은 컴퓨터를 구입하여 리눅스를 설치한 뒤 직접 미생물 유전체 해독 자료(fastq file)를 조립하는 일이 생길 것을 기대하는 것은 어려울 것이다. 그러나 이 일이 정말로 필요함을 느끼게 되고, 그 과정이 어떻게 돌아가는지를 최소한 교육장에서는 체험해 볼 수 있게 하는 것이 나의 목표이다.

NCBI에서 직접 내려받았거나 시퀀싱 서비스 업체에서 받은 미생물 genome assembly가 있다고 가정하자. Contig는 몇 개이고, 길이는 얼마나 되고, 예측된 유전자의 수는 어떻고... 이런 것들은 전부 리눅스에서 명령행 한 줄이면 해결이 될 일이다. 상용 genomics tool을 쓰지 않아도 얼마든지 가능하다. 이런 체험을 하기 좋은 환경은 무엇일까? 교육 기간이 끝나면 쓰지 못하는 서버에 의존할 것이 아니라, 수강생이 준비해 온 노트북 컴퓨터에서 간단히 해결할 방법은 없는 것일까? 

이번에도 Oracle VirtualBox를 생각해 보았지만, 명령행 환경만 체험하기 위한 용도로는 너무 무겁고, 다소 복잡해 보이는 세부 설정 과정이 초보자에게는 부담이 될 것이다. 그러다가 대안으로 발견한 것이 바로 Linux용 Windows 하위 시스템(Windows Subsystem for Linux, WSL)이었다. 사용이 몹시 불편하던 기존의 명령 프롬프트를 대체함과 동시에 WSL과 파웨셸까지 포함하게 된 Windows Terminal을 같이 사용하면 명령행 기반으로 돌아가는 리눅스 프로그램을 실행하는데 불편함이 없다. WSL 환경에 설치한 우분투에 Conda를 이용하여 실습에 필요한 프로그램과 데이터를 받은 뒤, 이를 tar 파일로 export하여 수강생에게 배포하면 된다. 동등한 수준의 프로그램과 데이터를 받은 VirtualBox 파일보다 크기도 훨씬 작은 것 같다. GUI를 갖춘 JAVA 프로그램은 윈도우에서 직접 돌리면 된다.


VirtualBox를 통해서 리눅스를 사용할 때 느꼈던 어색함도 없고, 상호 파일을 전송하기 위해 공유폴더를 설정할 일도 없다. 그저 윈도우에서 탐색기를 열고 주소창에 '\\wsl$'라고 치면 된다. 또는 리눅스 쪽의 명령행에서 'explorer.exe .'를 쳐도 된다.

호스트 측과 게스트 간에 서로 사용할 자원(CPU나 메모리 및 디스크 공간)에 대한 경계를 확실히 세워 두고 작동하는 VirtualBox(그것이 장점이 되고 또 필요한 순간도 있을 것임)과는 달리 윈도우 안에서 작동하므로 이질감이 적고, 시스템 자원에 대한 특별한 경계가 없으니 유연한 운용이 가능할 것이다. 물론 VirtualBox는 리눅스를 호스트로 하여 쓸 수도 있다는 것이 큰 장점이다.

이러다가는 Cygwin을 쓸 사람이 하나도 없을 것 같다. Windows 11부터는 리눅스용 GUI 프로그램도 돌릴 수 있다고 하니 말이다.

어제 오후부터 WSL을 익히면서 세상이 이렇게 편해졌다는 사실에 새삼 놀라고 감탄하였다. 윈도우의 환경변수(%USERPROFILE% 또는 PowerShell 입장에서는 $env:USERPROFILE)에 대한 지식은 그동안 백지 상태였다가 이번 기회에 그 의미를 제대로 익히게 되었다.

경계가 허물어지는 세상이다. 여기에 더하여 Docker에 익숙해지면 정말 재미난 세상이 펼쳐질 것 같다.


2022년 2월 21일 월요일

FASTA file을 탭으로 구분된 파일로 펼치기(fasta2tbl 스크립트 + 줄 번호 삽입)

그렇게 중요하지 않은 일에 과도하게 집착을 하는 것은 아닌지 모르겠다. 어제 올린 글('목록을 이용하여 FASTA file에서 원하는 서열을 뽑아내기')에서 이미 awk를 이용한 방법을 소개하였었다. 이것은 특별한 문제가 없이 잘 작동하지만, description 필드를 구성하는 과정에서 사소한 오류가 있다. Description을 포함하는 서열 ID 라인을 처리할 때, awk의 필드 변수인 $1을 null 문자로 바꾸어 버린다. 그러고 나서 $0을 출력하게 만들게 하였더니 $1만 제거되고 그 뒤를 따르는 공백 문자가 description의 시작 부분에 남아 보기가 싫다. 때로는 description 없이 서열 ID만 존재하는 FASTA file도 있을 것이다.

많은 부분을 차지하지 않는 예외 사항까지 잘 처리할 수 있는 스크립트를 짜는데 의외로 많은 시간이 걸린다. 그리고 나는 이러한 마무리 과정에서 오히려 큰 희열을 느낀다. 이것이 바람직한 코딩의 자세인지는 잘 모르겠다. 보람을 상회하는 희열 추구는 자가발전의 경향이 있어서 에너지를 많이 소모하게 만드는데...

인터넷을 뒤져서 가장 완성도가 높아 보이는 awk 스크립트를 구한 다음, 기능을 확장하였다. 아래의 것을 fasta2tbl이라는 파일로 저장한 뒤 처리할 FASTA file 이름을 인수로 주어서 실행하면 된다. 결과물은 서열 ID, description(없는 경우 'NO_DESC'로 표시), 그리고 서열로 이루어진 3개의 컬럼으로 표현된다.

#!/usr/bin/awk -f
# FASTA file을 seq_ID, description, sequence의 세 컬럼으로 이루어진 tab-delimted file로 전환한다.
# description이 없는 경우 두 번째 컬럼에 'NO_DESC'을 채워 넣는다.
{
        if ($1~/^>/) {
                id = $1;
                sub(/^>/, "", id);
                if (NF == 1 )
                    desc = "NO_DESC"
                else
                    desc = $0;
                    sub(/^>[^\s]+\s/, "", desc)
                if (NR>1)
                        printf("\n%s\t%s\t", id, desc)
                else
                        printf("%s\t%s\t", id, desc)
        } else {
                printf("%s", $0)
        }
} END { printf "\n" }

# adapted from Josep Abril's Fasta2Tbl
# https://bioinformatics.stackexchange.com/questions/2649/how-to-convert-fasta-file-to-tab-delimited-file

awk는 세미콜론을 쓰는 방식이나 중괄호의 사용 등이 그렇게 엄격하지는 않은 것 같다. 문법이 헐렁하면 사용자는 게을러지는데...

다음 명령어 한 줄이면 FASTA file을 TSV file로 전환한 뒤 라인 번호(with leading zeros)를 삽입한 네 컬럼짜리 결과 파일을 얻게 될 것이다.

$ fasta2tbl input.fa | nl -n rz - > SeqWithLineNum.tab

2022년 2월 20일 일요일

목록을 이용하여 FASTA file에서 원하는 서열을 뽑아내기 - 심화편

이 주제는 나의 블로그에서도 여러 차례 다루어진 바가 있고, 구글을 검색해 보면 꽤 많은 질문과 답이 올라와 있다. 파이썬을 능숙하게 다루는 사람이라면 별로 어렵지 않게 스크립트를 짜서 원하는 바를 해결할 수 있을 것이다. 나 역시 Perl을 이용하여 이따금씩 필요한 염기서열을 multi-FASTA file에서 추출하고는 한다.

KentUtils라 불리는 Jim Kent의 유틸리티 중에서 faSomeRecords(다운로드 링크)가 이러한 목적에 딱 맞는 '이미 잘 알려진' 프로그램이다.'-exclude' 옵션을 이용하면 목록 파일(listFile)에 있는 서열을 제거하는 정반대의 동작을 한다.

$ faSomeRecords 
faSomeRecords - Extract multiple fa records
usage:
   faSomeRecords in.fa listFile out.fa
options:
   -exclude - output sequences not in the list file.

BBMap에 포함된 유틸리티 중에서 Filterbyname.sh 또한 faSomeRecords와 거의 같은 동작을 한다. 단,  fastq 파일도 다룰 수 있다는 것은 장점이다.

# By default, "filterbyname" discards reads with names in your name list, and keeps the rest. To include them and discard the others, do this:

$ filterbyname.sh in=003.fastq out=filter003.fq names=names003.txt include=t

잘 알려진 프로그램을 사용하는 방법을 알아 보았으니 이제부터는 리눅스에 기본적으로 포함된 유틸리티(awk 등)을 사용하여 다소 원시적인 방법으로 똑같은 일을 해 보자. 단, 전제 조건을 살짝 비틀어 보았다. 추출할 서열의 ID가 아니라 일련번호(1부터 시작)를 별도의 목록 파일에 갖고 있다고 가정한다.

FASTA file의 서열 ID, description 및 염기서열(줄바꿈을 제거하여 한 줄로 펼쳐야 함 - 흔히 unwrapping이라 부른다)을 탭으로 구분된 TSV 파일을 만들고, 이를 목록 파일과 비교하여 공통된 컬럼을 갖는 라인을 join으로 병합한 뒤 해당 라인만 추출하여 다시 FASTA file로 전환하는 것이 매우 일반적인 방법이다. 이러한 방법을 쓰는 경우에는 꼭 서열 일련번호가 아니라 서열 ID를 써서 원하는 서열을 뽑는 것도 무난하게 할 수 있다.

최근 알게 된 long read용 binning program인 LRBinner의 출력 파일('Bins.txt') 정보를 이용하여 입력 FASTA file에서 각 bin에 해당하는 염기서열을 별도의 파일로 뽑아내는 방법을 알아보도록 한다. Bins.txt는 입력 파일에 실린 각 염기서열이 속하는 bin 번호를 다음과 같이 보여준다. 참 건조하고 재미 없는 포맷이다. 

2
1
3
2
0
...

2번 bin에 속하는 염기서열은 1번과 4번에 해당한다. 2라는 값을 갖는 라인의 번호를 뽑아서 LineNum.txt 파일에 저장하려면 다음과 같이  awk 명령어를 실행한다.

awk '$1==2{print NR}' Bins.txt > LineNum.txt

nr이라는 유틸리티는 인수로 주어진 텍스트 파일을 라인 번호와 함께 출력한다. 당장 아무 텍스트 파일이나 가져다가 'nr file.txt'를 실행해 보라. 행 시작 부분의 공백을 제거하고 구분자를 탭으로 전환하려면 다음과 같이 약간의 수정을 거쳐야 한다. 맨 왼쪽에 라인 번호를 삽입한 텍스트 파일은 앞으로도 종종 쓸모가 있을 것이다.

nl -n ln -s $'\t' Bins.txt > BinWithLineNum.txt

이상에서 설명한 방법으로 추출해야 할 염기서열의 번호를 별도의 파일('LineNum.txt')에 저장하였다고 하자. 그러면 우선 입력물인 input.fa 파일을 TSV 파일('SeqWithNum.tab')로 펼치되, 첫 번째 컬럼에는 1에서부터 시작하는 서열 일련번호를 넣기로 한다.

awk '/^>/ { id=$1; 
            sub(/^>/,"",id); 
            $1=""; 
            printf("%s%d\t%s\t%s\t",(N>0?"\n":""),N+1,id,$0);
            N++;
            next; }
          { printf("%s", $0); } 
          END {printf("\n");}' input.fa > SeqWithNum.tab

awk  명령이 좋은 점은 따옴표로 둘러싼 긴 명령어 내부에서 가독성을 위해 줄바꿈을 해도 된다는 것. 단, END와 {..} 사이에 줄바꿈을 넣으면 에러가 남을 확인하였다. BEGIN과 {..} 사이도 마찬가지일 것이라고 생각한다. 이 awk 명령어를 여기에서 설명하지는 않겠다. 첫 번째 컬럼에 일련번호를 삽입하고,  descripton에 해당하는 여러 단어를 하나의 컬럼에 들어가게 하는 등 꽤 많은 공을 들였다. 하지만 이 awk 명령어는 나의 창작은 아니다. 또한 이것이 awk를 이용하여  FASTA file을 펼치는 유일하거나 가장 능률적인(짧은?) 방법도 아니다. 가끔씩 이 awk 명령어를 바라보면서 왜 이렇게 짜여졌는지 생각을 하는 것도 좋은 공부가 될 것이다.

LineNum.txt와 SeqWithNum.tab이라는 두 개의 파일이 만들어졌다. LineNum.txt 파일에 수록된 라인 번호에 해당하는 것을  SeqWithNum.tab에서 골라낸 뒤, 이를  FASTA file로 전환하는 것이 오늘의 목표이다. SeqWithNum.tab 파일의 첫 번째 컬럼에 1로부터 시작하는 라인 번호가 들어 있어서 점검하기는 편하다. 

오늘의 최종 목표를 달성하는 방법도 몇 가지가 있을 것이다. 우선 첫 번째 방법을 소개한다. 사실 이 방법을 쓰는 경우라면 SeqWithNum.tab 파일에 라인 번호를 별도의 컬럼으로 넣지 않아도 된다. 자세한 설명은 생략한다.

awk 'FNR==NR{a[$1];next}(FNR in a){print}' LineNum.txt SeqWithNum.tab

표준 출력으로 SeqWithNum.tab에서 조건에 맞는(즉 LineNum.txt에서 라인 번호를 지정한) 라인이 나올 것이다. 따라서 이를 FASTA file로 전환하여 저장하려면 파이프('|') 기호 뒤에 다음 명령어를 넣어야 한다(A).

awk -F"\t" '{printf(">%s %s\n%s\n",$2,$3,$4)}' | seqret fasta::stdin fasta:outfile.fa

왼쪽의 awk 명령어만을 사용해도 FASTA 기본 골격으로 출력은 된다. 그러나 서열 영역은 한 줄로 펼쳐진 형태가 되므로 이를 EMBOSS의 seqret로 처리하여 wrapping을 실시하였다.

최종 목표를 달성하는 두 번째 방법을 소개한다. 그것은 바로 sort와 join 명령어를 이용하는 것인데, 검색을 해 보면 두 텍스트 파일에서 특정 컬럼을 서로 비교하여 같은 값을 갖는 라인(혹은 선별된 컬럼)을 추출하는 일반적인 솔루션으로 소개하고 있다. 예를 들어서 두 개의 파일 file_A와 file_B가 있고, 각각의 첫 번째 컬럼이 같은 값을 갖는 경우 file_A의 1, 2번째 컬럼과 file_B의 4, 5번째 컬럼을 추출하 싶다면 다음과 같이 해야 한다. 컬럼은 탭 문자로 구분되어 있다고 가정한다.

join -t $'\t' -1 1 -2 1 -o 1.1,1.2,2.3,2.5 <(sort -t $'\t' -k1,1 file_A) <(sort -t $'\t' -k1,1 file_B)

생각보다 꽤 복잡하다. 리눅스 shell에서 컬럼 구분자를 탭문자로 설정하는 것(-t $'\t')도 그렇고, file_A와 file_B가 join field(같은 값을 갖는지 비교할 컬럼)에 대해서 이미 sort가 되어 있다 해도 여기에서 보인 process substitution에 의해서  join  명령어 내에서 sort를 다시 하지 않으면 경고문을 발생한다. '--nocheck-order' 옵션을 주면 될 것 같으나 제대로 작동을 하지 않아서 몹시 마음에 들지 않는다. sort나 join 모두 사용법을 정확히 알지 않으면 낭패를 겪기 쉬운 명령어이다.

따라서 LineNum.txt와 SeqWithNum.tab의 두 파일이 준비된 경우에 각 파일의 첫 번째 컬럼을 비교하여 같은 값을 갖는 라인을 SeqWithNum.tab에서 추출하여 FASTA file로 되돌리는 '두 번째의 복잡한' 방법은 다음과 같다.

join -t $'\t' -1 1 -2 1 <(sort -t $'\t' -k1,1 LineNum.txt) <(sort -t $'\t' -k1,1 SeqWithNum.tab)

이 명령어 뒤에 파이프 기호를 쓴 뒤 위에서 소개한 (A) 명령어를 써 넣어야 FASTA file로 결과가 저장된다. 이 방법은 아직 완벽하지 않다. process substitution 내에서 일반 문자의 순서로 소팅했기 때문에 라인 넘버로는 2, 3, 15 순서가 아니라 15, 2, 3...의 순서로 출력된다. sort 명령어의 단독 실행에서는 -k1n,1 또는 -nk1,1 옵션을 제공하여 숫자 기준의 정렬이 되지만, 이상하게도 join 명령어 안에서 <(sort ...) 형태로 실행하면 sort가 되지 않았다는 오류 메시지가 나온다. 따라서 단순하게 -k1,1  또는 -k1 옵션을 주어서 문자 순서로 sort를 해야만 한다.

오늘 학습의 교훈:

  • Don't reinvent the wheel. 창의력을 시험하고 싶다면 바퀴를 다시 발명해도 좋다.
  • awk에서 가끔은 printf() 함수를 사용할 필요가 있다.
  • sort 명령의 -k(key)는 시작과 끝을 지정하는 키 2개를 지정하는 것이 기본이 되어야 한다.
  • 예전에 쓴 글 Awk를 이용한 간단한 텍스트 파일의 조작(join)을 가끔 들추어 보면서 awk  명령어의 기본을 잊지 말도록 하자.

 

Numeric sort 문제가 없는 join 방법

라인 번호를 1, 2, 3이 아니라 00001, 00002, 00003과 같이 자릿수를 꽉 채운 0('leading zero')이 있는 것으로 바꾸면 해결이 된다. 단, seq_id, description 및 서열의 3개 컬럼으로 이루어진 TSV 파일을 먼저 만든 뒤  nl 명령을 사용하여 첫 컬럼에 라인 번호를 넣는다. nl 명령은 라인 번호를 오른쪽으로 정렬한 뒤 leading zero를 넣는 옵션('-n rz')이 있다. 전체 명령어는 다음과 같다.
 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
# 3개 컬럼으로 구성된 파일을 만든다.          
$ awk '/^>/ { id=$1; 
             sub(/^>/,"",id); 
             $1=""; 
         # 아래 코드는 seq_id, description, sequence의 3개 컬럼을 만든다,
             printf("%s%s\t%s\t",(N>0?"\n":""),id,$0);
         # 아래 코드는 1로부터 시작하는 번호를 첫 번째 컬럼에 넣는다(총 4개 컬럼).
         #   printf("%s%d\t%s\t%s\t",(N>0?"\n":""),N+1,id,$0);
             N++;
             next; }
           { printf("%s", $0); } 
           END {printf("\n");}' input.fa > Seq.tab
           
# 라인 번호를 첫 컬럼으로 삽입한다.
$ nl -n rz -s $'\t' Seq.tab > SeqWithNumLeadingZeros.tab

# LineNum.txt의 숫자에 leading zero를 채운다. 
# 그러려면 SeqWithNumLeadingZeros.tab에서 자릿수가 얼마나 필요한지를 알아내야 한다.       
$ for n in $(cat LineNum.txt)
> do
> printf "%06d\n" $n >> LineNumLeadingZeros.txt
> done

# join 실행(sort가 필요하지 않음)
$ join -t $'\t' LineNumLeadingZeros.txt SeqWithNumLeadingZeros.tab |\
    awk -F"\t" '{printf(">%s %s\n%s\n",$2,$3,$4)}' |\
    seqret -auto fasta::stdin fasta:outfile.fa

SeqKit의 fx2tab이라는 명령어를 쓰면 FASTA file을 TSV file로 전환할 수 있다. 물론 내가 원하는 형태와는 약간 다르다. 

어쩌면 나는 바퀴를 새로 발명하려고 애를 쓰고 있었는지도 모른다...

How to convert fasta file to tab delimited file

1
2
3
4
from Bio import SeqIO

for record in SeqIO.parse('input.fa', 'fasta'):
    print('{}\t{}\t{}'.format(record.id, record.description, record.seq))                                                                          

2022년 2월 16일 수요일

PyTorch를 이용하여 long read binner(LRBinner) 실행하기 - 일단은 실패!

Nanopore read 혹은 이를 대상으로 만들어진 조립물에 대하여 직접적으로 사용할 수 있는 metagenomics binning tool이 없을까? 지금까지 알려진 대부분의 도구는 short read를 대상으로 하였다. 이를 long read 혹은 그로부터 만들어진 assembly에 적용하기에는 무리가 있다. 왜냐하면 거의 필수 정보로 여겨지는 coverage information을 long read에 대해서 얻는 것은 - 불가능한 것은 아니지만 - 적합하지 않기 때문이다.

비교적 초창기에 개발된 소프트웨어인 MyCC를 잠시 테스트해 본 일이 있다. Coverage information을 주지 않고도 썩 잘 작동한다고 알려져 있지만, 요즘에는 사용자들의 관심 밖으로 밀려난 것만 같은 느낌이 든다. 계속 검색을 해 보니 MetaBCC-LR이라는 것이 눈에 뜨였다. 논문을 찾아서 읽어 보았다. k-mer와 trinucleotide composition 등 다양한 자료를 수집하여 bin 수립을 위한 모델을 만들어서 read를 분류한 뒤, 이를 별도로 조립하면('partitioned assembly') 샘플 속에 섞인 다양한 미생물의 유전체를 나누어서 재구성할 수 있다는 것이다.

이를 테스트하려고 소스를 컴파일하면서 잘 안되는 부분을 해결하려 애쓰다가 MetaBCC-LR의 저자 4명 중에서 3명이 모여 만들어낸 LRBinner라는 것까지 발견하게 되었다. 이는 MetaBCC-LR의 개선판으로서, read뿐만 아니라 contig도 입력물로 사용할 수 있다고 하였다. 개별 연구실에서 nanopore 장비를 직접 구입하여 응용의 폭을 점점 넓혀가는 요즘 시대에 아주 적절한 도구가 아닐 수 없다.

GitHub를 참조하여 conda에 환경을 만든 뒤 파이썬 소스를 클론하여 빌드를 하였다. 개발자가 제공하는 샘플 데이터를 이용하여 테스트 실행을 하는데 '--cuda'라는 옵션이 말썽을 부린다.

"raise AssertionError("Torch not compiled with CUDA enabled")"

그러면 이 옵션을 빼고 실행하면 되지 않겠나... 성공적으로 실행을 완료하였다. 

자, 그러면 최근 구입한 GeForce RTX 3080 장착 컴퓨터에 이를 설치하면 GPU를 이용한 계산이 되지 않을까? GitHub에서 알려준 방법 그대로 새 컴퓨터에 LRBinner를 설치하여 테스트 데이터를 분석해 보았다. 역시 같은 문제가 발생하였다. nvidia-cuda-toolkit 패키지가 설치되지 않았다는 사실을 발견하고 apt로 이를 설치한 뒤 LRBinner를 깔았으나 역시 마찬가지의 문제가 발생하였다. Guppy basecaller 실행에는 nvidia-driver-510을 설치하는 것으로 충분하였기에(이때 자동으로 cuda driver도 설치되는 것으로 알고 있음) 그때는 잘 몰랐었다.

PyTorch는 요즘 인기를 얻고 있는 오픈소스 머신 러닝 라이브러리로서, GPU 사용이 가능하다. 사실 Pytorch라는 것이 있다는 것도 오늘 처음 알았다! 윈도우에 PyTorch 설치, GPU 설정, 자세하게라는 글을 참조하여 다음과 같이 실행해 보았다.

$ conda create -n lrbinner -y python=3.7 numpy scipy seaborn h5py tabulate hdbscan gcc openmp tqdm biopython # 일단은 pytorch를 포함시키지 않음
$ conda activate lrbrenner
(lsbinner) $ conda install pytorch torchvision cudatoolkit=10.1 -c pytorch
$ python
Python 3.7.12 | packaged by conda-forge | (default, Oct 26 2021, 06:08:21) 
[GCC 9.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import torch
>>> torch.cuda.get_device_name(0) 
'NVIDIA GeForce RTX 3090'
>>> torch.cuda.is_available() 
True
>>> torch.__version__
'1.4.0'
>>> quit()
(lrbinner) $ git clone https://github.com/anuradhawick/LRBinner.git
(lrbinner) $ cd LRBinner; python setup.py build

LRBinner가 잘 돌아갈 수 있게 PyTorch가 깔린 것 같다. 그런데 이번에는 다른 에러가 난다.

RuntimeError: cuDNN error: CUDNN_STATUS_EXECUTION_FAILED

으으으... 이건 또 뭔가? 어쩔 도리가 없이 '--cuda' 옵션을 제거한 뒤 다시 LRBinner를 실행하여 결과를 얻을 수 있었다.

기왕 장착한 GPU를 이용하여 더 많은 일을 해 보고자 하였는데 처음부터 완벽할 수는 없을 것이다. 일단 guppy는 잘 돌아가고 있고, 비록 GPU를 쓰는 상황은 못 되었지만 long read(or assembly)를 대상으로 하는 최신의 metagenomics binning tool을 쓰게 되었으니 완벽한 실패는 아닌 셈이다.


2022년 2월 13일 일요일

RTC1302를 이용한 아두이노 시계 코드

Michael Miller("Makuna")의 Rtc 라이브러리에 포함된 예제(DS1302_Simple)를 활용한 아두이노 시계의 코드를 싣는다. 택트 스위치를 누르면 인터럽트가 발동하여 12시-24시 표시를 전환하게 만드느라 이곳 저곳을 참조하여 코드 조각을 가져다가 일단은 돌아가게 만들어 놓았다.

// CONNECTIONS:
// DS1302 CLK/SCLK --> 7
// DS1302 DAT/IO --> 6
// DS1302 RST/CE --> 5
// DS1302 VCC --> 3.3v - 5v
// DS1302 GND --> GND

#include <ThreeWire.h>  
#include <RtcDS1302.h>
#include <LiquidCrystal.h> 
#define debounceTime 75 // <<- Set debounce Time (unit ms)

int interruptPin1 = 3; // 2 or 3
volatile int state = LOW;

char *dayArray[] = {
  "Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat"
};

ThreeWire myWire(6,7,5); // IO, SCLK, CE
RtcDS1302<ThreeWire> Rtc(myWire);
LiquidCrystal lcd(8, 9, 10, 11, 12, 13); //RS, EN, data(4,5,6,7)

void setup () 
{
    Serial.begin(9600);

    Serial.print("Compiled: ");
    Serial.print(__DATE__);
    Serial.print("\t");
    Serial.println(__TIME__);

    lcd.begin(16, 2);
    lcd.clear();

    pinMode(interruptPin1, INPUT_PULLUP); // pulled down using a 10K resistor
    attachInterrupt(digitalPinToInterrupt(interruptPin1), buttonPushed, FALLING);
    
    Rtc.Begin();

    RtcDateTime compiled = RtcDateTime(__DATE__, __TIME__);

    if (!Rtc.IsDateTimeValid()) 
    {
        // Common Causes:
        //    1) first time you ran and the device wasn't running yet
        //    2) the battery on the device is low or even missing

        Serial.println("RTC lost confidence in the DateTime!");
        Rtc.SetDateTime(compiled);
    }
    if (Rtc.GetIsWriteProtected())
    {
        Serial.println("RTC was write protected, enabling writing now");
        Rtc.SetIsWriteProtected(false);
    }
    if (!Rtc.GetIsRunning())
    {
        Serial.println("RTC was not actively running, starting now");
        Rtc.SetIsRunning(true);
    }

    RtcDateTime now = Rtc.GetDateTime();

    if (now < compiled) 
    {
        Serial.println("RTC is older than compile time!  (Updating DateTime)");
        Rtc.SetDateTime(compiled);
    }
    else if (now > compiled) 
    {
        Serial.println("RTC is newer than compile time. (this is expected)");
    }
    else if (now == compiled) 
    {
        Serial.println("RTC is the same as compile time! (not expected but all is fine)");
    }
}

void loop () 
{
    RtcDateTime now = Rtc.GetDateTime();
    printDateTime(now);

    if (!now.IsValid())
    {
        // Common Causes:
        //    1) the battery on the device is low or even missing and the power line was disconnected
        //Serial.println("RTC lost confidence in the DateTime!?!");
    }

    delay(1000); // one seconds
}

void printDateTime(const RtcDateTime& dt)
{
    char daystring[16];
    char timestring[16];
    String flag = "AM";
    String stateStr = "24-hr system";
    short DOW = dt.DayOfWeek();
    
    sprintf(daystring, " %4u-%02u-%02u %s", dt.Year(), dt.Month(), dt.Day(), dayArray[DOW]);

    short Hour = dt.Hour();

    if ( state ) { // 24-hr system
      lcd.clear();
      sprintf(timestring, "    %02u:%02u:%02u", dt.Hour(), dt.Minute(), dt.Second());
    } else {  // 12-hr system with AM/PM designation
       if (dt.Hour() > 12) {
         Hour = dt.Hour() - 12;
         flag = "PM";
       }     
       stateStr = "12-hr system with AM/PM designation";
       sprintf(timestring, "  %2d:%02u:%02u <%s>", Hour, dt.Minute(), dt.Second(), flag.c_str());
    }
    Serial.print("State ");
    Serial.print(state);
    Serial.print(": "); 
    Serial.println(stateStr);
    lcd.setCursor(0, 0);
    lcd.print(timestring);
    lcd.setCursor(0, 1);
    lcd.print(daystring);
}

void buttonPushed() {
    static unsigned long lastTime = 0;
    unsigned long Now = millis();
    if((Now - lastTime) > debounceTime)
    {
        state=!state;
        Serial.print("Button pressed: state changed to ");
        Serial.println(state);
        RtcDateTime now = Rtc.GetDateTime();
        printDateTime(now);
    }
    lastTime = Now;
}

C++ 함수 선언부에서 참조를 통하여 파라미터를 전송하는 방법을 겨우 이해하는 수준이다. 따라서 추가로 장착한 버튼을 눌러 시간 설정을 고치도록 개선한 버전 2 코드를 만들려면 아직 갈 길이 멀다. 그 다음 목표는 알람 기능 추가하기.

12시-24시 전환은 3번 버튼이 담당한다.

버튼을 이용한 시 변경 기능은 Accurate Clock Just Using an Arduino을 참조할 예정이다.  paulsb가 공개한 이 프로젝트에서는 그러나 RTC 모듈이라는 외부 하드웨어를 전혀 사용하지 않기 때문에 전원을 새로 연결할 때마다 시각을 새롭게 맞추어야 한다. Makuna의 Rtc 라이브러리를 가져다가 RTC 모듈에 데이터를 써 넣고 읽어오는 기능을 만들어 넣는 일이 필요하다. paulsb의 코드는 시간과 관련한 특별한 객체를 전혀 사용하지 않아서 시간 증가 및 표현을 전부 계산을 통해 구현한다. 공부 목적으로 살펴볼 가치가 있다.

LCD 또는 7-세그먼트 표시기를 이용한 시계 만들기는 아두이노 입문자의 프로젝트로 널리 인식되고 있다. 매우 단순한 목표라고 생각할 수도 있으나 시간이라는 데이터 구조가 의외로 복잡하여 공부를 할 것이 많다. 1970년 1월 1일 0시라는 유닉스 세계의 epoch time이 정해진 이유(사실 명확한 것은 없음) 등에 대해서도 흥미로운 이야깃거리가 많을 것이다.

2022년 2월 11일 금요일

숙원사업 - 드디어 GPU가 장착된 컴퓨터를 구입하여 나노포어 시퀀싱 준비를 완료하다

이제는 다른 연구팀의 장비를 빌리거나 CPU에서 극도로 느린 guppy basecalling을 하면서 며칠이고 기다릴 필요가 없게 되었다. NVIDIA GeForce RTX 3090이 장착된 데스크탑 컴퓨터를 구매하였기 때문이다. 메모리는 64 기가바이트로 맞추어서 주문하였다. 우분투 20.04의 설치용 USB 매체를 늘 갖고 다니기 때문에 아주 쉽게 리눅스부터 설치를 하였다. 

모니터 뒤편에는 같은 케이스에 만들어진 AMD Ryzen 5950X가 놓였다. AMD 머신에 비하면 말할 수 없이 조용하다! 

무슨 유흥업소 간판도 아니고 이렇게 요란하게 그래픽 카드를 만들다니...

NVIDIA 드라이버와 CUDA toolkit을 수동으로 설치하려다가 이상하게 꼬여서 고생을 하였다. 우분투를 싹 새로 설치한 다음, ubuntu-drivers devices 명령이 제안하는 드라이버를 찾아서 설치하니 아주 쉽게 CUDA까지 해결이 되었다.

$ ubuntu-drivers devices
== /sys/devices/pci0000:00/0000:00:01.0/0000:01:00.0 ==
modalias : pci:v000010DEd00002204sv00001458sd00004043bc03sc00i00
vendor   : NVIDIA Corporation
driver   : nvidia-driver-470 - distro non-free
driver   : nvidia-driver-510 - distro non-free recommended
driver   : nvidia-driver-470-server - distro non-free
driver   : xserver-xorg-video-nouveau - distro free builtin
$ sudo ubuntu-drivers autoinstall
# or sudo apt install nvidia-driver-510
$ sudo shutdown -h reboot

설치 상태를 확인해 보았다. "NVIDIA GeForce RTX 3090"이라는 명칭이 잘 보인다.

$ nvidia-smi  --query | grep "Product Name"
    Product Name                          : NVIDIA GeForce RTX 3090
hyjeong@ionic:~$ nvidia-smi
Sat Feb 12 20:50:46 2022       
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 510.47.03    Driver Version: 510.47.03    CUDA Version: 11.6     |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp  Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|                               |                      |               MIG M. |
|===============================+======================+======================|
|   0  NVIDIA GeForce ...  Off  | 00000000:01:00.0 Off |                  N/A |
|  0%   32C    P8    15W / 370W |     68MiB / 24576MiB |      0%      Default |
|                               |                      |                  N/A |
+-------------------------------+----------------------+----------------------+
                                                                               
+-----------------------------------------------------------------------------+
| Processes:                                                                  |
|  GPU   GI   CI        PID   Type   Process name                  GPU Memory |
|        ID   ID                                                   Usage      |
|=============================================================================|
|    0   N/A  N/A      1098      G   /usr/lib/xorg/Xorg                 56MiB |
|    0   N/A  N/A      1240      G   /usr/bin/gnome-shell                9MiB |
+-----------------------------------------------------------------------------+

다음으로는 나노포어 커뮤니티 사이트에 접속하여 MinKNOW와 Guppy GPU 버전을 설치하였다. 오랜만에 커뮤니티 사이트를 들어갔더니 Firefox(리눅스)와 Chrome(윈도우)에서 보이는 창이 달라서 소프트웨어 다운로드 위치를 찾느라 애를 먹었다. 더욱 혼동스러웠던 것은 'MinION Software' - 'MinIT Software' - 'MinKNOW Stand Alone GUI'라고 표현되어 있다는 점이다. 셋 중에서 가장 앞에 있는 것이 MinION Mk1B를 위한 것인데, 내 기억으로 전에는 분명히 MinKNOW라고 되어 있어서 어느 것을 골라야 하는지 주저할 필요가 없었다.

MinION Mk1B를 연결하여 인식이 됨을 확인하였다.

다음주에는 작년에 fast basecall로 생산해 둔 시퀀싱 결과를 guppy GPU version으로 다시 처리할 수 있을 것이다. 연구비에 여유가 있어서 이렇게 장비를 마련하게 된 것을 감사히 생각하며...

Torsten Seemann의 nullarbor 테스트

호주 관광청의 눌라보 평원 관광 안내(링크)
앞으로 96 킬로미터를 더 달려야 캥거루를 만날 수 있다고? 단 한번도 굽은 길을 만나지 않고? 

nullarbor를 conda로 한꺼번에 설치하는 것은 쉽지 않았다. 어제 저녁부터 명령을 걸어 놓았는데 도무지 진도가 나가질 않는다. 소스는 git clone으로 복사해 놓고, 필요한 패키지를 conda install로 하나씩 설치하였다. 약간의 오타가 있지만 애교로 봐 주자. 예를 들어 newick-utils라는 dependency가 있는데, 이는 newick_utils를 잘못 적은 것이었다.

Nullarbor GitHub 웹사이트에서 기술한 모든 dependency(링크)를 수작업으로 전부 설치하였다. 새로운 환경을 만들지는 않았고, 이전에 만들어 둔 snippy 환경에 하나씩 더해 나간 다음 환경 이름을 아예 바꾸어 버렸다. 이름은 prokaryotic whole genome analysis의 약자인 pwga로 하였다. Conda는 rename이라는 유틸리티를 제공하지 않기 때문에, 이미 설치한 환경의 이름을 바꾸려면 조금 까다롭다(참고의 글).

패키지를 하나씩 설치하다가 맨 마지막인 centrifuge 단계에서 에러가 발생했다. 파이썬 버전이 맞지 않는다는 것이었다. Centrifuge는 꽤 오래전에 만들어진 것이라서 데이터베이스 업데이트를 제외하면 더 이상 개선이 이루어지지 않는 것 같다. 이미 설치된 파이썬 버전을 3.7에서 3.6으로 내려서 해결을 했다. Cetrifuge의 개발자인 김대환 교수는 미국에서 활동하는 한국인 생명정보학자이다. TopHat의 개발자라고 하면 더 빨리 기억할 수 있을 것이다. 구글 스칼라 정보를 보니 인용 횟수 역시 대단하다. 아마 Steven Salzberg가 그의 지도교수였던 것 같다.

필요한 것을 다 설치했는데 이제는 fq, fa라는 실행파일이 없다고 한다. 이건 또 뭔가? Conda로 배포되는 패키지도 아니고... 크! nullarbor 배포본의 bin 디렉토리에 있는 자체 스크립트였다. 필요한 것들이 다 설치되었는지 확인하려면 'nullarbor.pl --check'를 실행하면 된다. 'You deserve a medal!'이라는 익살맞은 메시지와 함께 준비가 다 되었음을 알려주었다. 앞으로 pwga 환경에서 미생물 유전체를 다루는 대부분의 일을 해도 될만큼 많은 프로그램이 한꺼번에 설치되었다.

테스트 실행을 해 보았다. Roary 때문에 최소 4개 isolate의 시퀀싱 자료가 필요하다고 한다. 탭으로 구분된 샘플 설명 파일을 만든 뒤 실행 준비를 시키면 이런 메시지가 나온다.

Run script를 별도로 만들어 주는 것이 아니라 'nice make all...' 명령어를 쓰라는 것이다. make는 프로그램 빌드를 할 때나 쓰는 것 아니었나? 뭔가 있어 보이는 실행 방법이라는 생각이 든다.

최종적으로 만들어진 html 리포트는 Tormes의 그것보다 좀 더 전문적으로 보인다. 단, 구식 Kraken(DB: minikraken_20171019_8GB.tgz)를 이용한 k-mer based species identification에서는 결과가 나오지 않았다. Kraken 2 및 centrifuge를 이용하도록 설정을 바꾸어 보도록 하자. Assembler 등 사용자가 선택할 수 있는 것들이 있는데 이를 어디에서 건드려야 하는지는 확실하지 않다. 어디 보자... 아, Makefile에서 plugin이라는 이름으로 바꿀 수가 있다. 다음에 Makefile의 일부를 보였다. "TAXONER"의 오른쪽 항목을 다른 것으로 바꾸면 원하는 프로그램을 쓰게 되는 것이다.

NAME := TEST
GCODE := 11
MIN_CTG_LEN := 500
PUBLISH_DIR := /home/hyjeong/public_html/MDU
ASSEMBLER := cpus=$(CPUS) opts="" /data/apps/nullarbor/bin/../plugins/assembler/skesa.sh
TREEBUILDER := cpus=$(CPUS) opts="" /data/apps/nullarbor/bin/../plugins/treebuilder/iqtree_fast.sh
TAXONER := cpus=$(CPUS) opts="" /data/apps/nullarbor/bin/../plugins/taxoner/centrifuge.sh
ANNOTATOR := cpus=$(CPUS) opts="" /data/apps/nullarbor/bin/../plugins/annotator/prokka_fast.sh
NW_DISPLAY := nw_display -S -s -w 1024 -l 'font-size:12;font-family:sans-serif;' -i 'opacity:0' -b 'opacity:0' -v 16
SNIPPY := snippy --force
SNIPPYCORE := snippy-core
ROARY := roary -v
ABRICATE := abricate --threads $(CPUS)

내가 nullarbor를 설치한 것은 centrifuge를 좀 쉽게 써 보려는 이유도 있다. 마치 zga 내부에서 checkm을 이용할 수 있는 것처럼 말이다.
 

쿠팡에서 전자 부품을 사면 이런 일이 생긴다

아두이노를 이용한 시계 만들기 프로젝트를 진행하다가 시각 설정(또는 알람 설정) 기능을 추가하고 싶어서 택트 스위치를 이용하기로 하였다. 아두이노 스타터 키트에는 4개의 택트 스위치가 들어 있었는데, 라즈베리 파이(볼루미오 3 설치)에 셧다운용 스위치를 만들어 주느라 하나를 쓰는 바람에 수량이 부족한 상태가 되었다. 별 생각을 하지 않고 쿠팡에서 택트 스위치를 10개 주문하여 어제 받았는데... 딱 파리 한 마리만한 스위치가 봉투에서 나온다. 재래식 화장실 주변에서 노니는 약간 큰 사이즈의...

엥? 뭐가 이렇게 작아?

DIP 패키지의 IC와 동등한 크기였다.
IR2153(데이터시트)이 오랜만에 사진에 찍혔다. 600V를 다룰 수 있는 self-oscillating half-bridge driver 칩이다. 이것을 보고 있노라니 갑자기 전압 부스터를 만들고 싶어졌다. 전에 실험을 할 때에는 상용전원 220V를 직접 정류하여 다루느라 칩을 종종 망가뜨렸었다. 차라리 요즘 쉽게 구할 수 있는 DC 12V 전원을 입력으로 이용하여 250V 정도를 만들게 하면 어떨까? Elektor magazine에서 진공관 앰프를 위한 직류 컨버터 자작 기사를 발견하여 여기에 소개한다. 히터를 제외하고 50와트 출력을 뽑을 수 있다면, 늘 꿈꾸어 오는 6V6 싱글 정도는 구동을 너끈히 할 수 있을 것이다. 회로도 등의 상세 자료를 보려면 로그인이 필요하다.

Fairly easy to build DC-DC converter for valve amplifiers using through hole components (just one SO-8 IC). Galvanically isolated high voltage output. The transformer must be created by oneself and its turns ratio sets output voltage. Input voltage is a 5A/ 12VDC AC adapter, output power is 50 W max.

오른쪽의 것이 이번에 구입한 것.

에구구... 쓰지 못할 물건은 아니지만 이렇게 작으니 다루기가 나쁘다. 온라인 쇼핑으로 모든 것을 다 해결할 수 있을 것만 같지만 현실은 그렇지 않다. 직접 눈으로 보는 것이 중요하다. 영상으로 보는 것은 주변 맥락을 제거해 버리므로 현실을 그대로 반영하지 않는다.

2022년 2월 15일 업데이트

National Valve Museum 웹사이트에서 DC-DC booster module을 이용한 흥미로운 앰프 제작 기사를 보았다. "What about driving a valve amplifier from an old laptop power supply? Finding G A French's push-pull EL91 amplifier set me thinking." 닉시(Nixie)관을 이용한 자작을 위하여 높은 DC를 만들어 내는 모듈이 팔리고 있으니, 진공관 앰프가 요구하는 공급 전력이 높지 않다면 이를 이용해도 된다는 것이다. 




2022년 2월 10일 목요일

박테리아 유전체 조립 및 분석용 파이프라인, Nullarbor와 ZGA

나는 아직도 논문을 종이에 인쇄하여 줄을 쳐 가면서 읽는 것을 좋아한다. 모니터 화면이나 태블릿 PC에 PDF 파일을 띄워 놓고 읽는 것이 더 편하다는 사람도 있고, 나 역시 앞으로는 종이 인쇄물보다 컴퓨터 화면을 이용하는 일도 점점 더 많아질 것이다. Mendeley에는 어디서든 볼 수 있도록 클라우드 서버에 동기화한 PDF 파일이 수북하게 존재한다지만, 최근에 읽고 중요하다고 생각한 논문을 다시 찾으려면 결국 인쇄하여 놓은 논문 더미를 뒤지는 것이 더 빠르다. 

박테리아 유래의 raw sequencing 자료를 조립하고 기본적인 분석을 실시하는 공개 소프트웨어로서 TORMES를 즐겨 사용하는 편이다. 단, 명령어를 날리기 위해 메타데이터 파일을 만드는 일이 조금 귀찮고, 서로 관련성이 없는 미생물 자료를 한번에 처리하는 것은 넌센스이다(왜냐하면 pan genome analysis를 기본으로 수행하므로). KRAKEN DB 때문에 설치 용량도 조금은 큰 편이다. 이보다 더 심플하고 가벼운 파이프라인이 없는지 최근 자료를 뒤져보니 바이오아카이브에 ZGA라는 소프트웨어가 소개된 것을 발견하게 되었다(논문, GitHub). 논문의 제목은 "ZGA: a flexible pipeline for read processing, de novo assembly and annotation of prokaryotic genomes"이다. 제목에 잘 나타나 있듯이 딱 유전체 주석화까지만 해 준다. pan genome analysis, core/accessory gene-base tree buidling, MLST, VF/AMR prediction 등은 다루지 않는 가벼운 파이프라인이다. 

개발자인 러시아의 Korzhenkov는 많은 미생물 유전체 조립 파이프라인이 명멸하였으나 현재까지도 계속 업데이트가 지속되는 것은 TORMES와 nullarbor뿐이라고 하였다. 2015~2016년 무렵에 A5-miseq이라는 일루미나 전용의 조립 파이프라인을 썼던 기억이 난다. 현재 인용 횟수는 814회. Nullarbor가 뭘까? 이 블로그에서 여러 차례 소개했던 호주의 과학자 Torsten Seemann의 작품이었다(GitHub). 

Nullarbor는 원래 호주 남부의 평원 지대를 뜻한다. arbor(나무)가 없어서(null) 이런 이름이 만들어진 것인지? 호주 관광청에서 소개한 눌라보 평원 관련 기사를 보라. 자동차로 여기를 여행하는데 무려 6일이 걸린다고 한다. 오늘 검색을 해 보기 전까지는 이 단어가 무엇을 뜻하는지 전혀 알지 못했다. 스테이크 전문점으로만 알고 있었던 아웃백(Outback)이라는 낱말도 '호주의 건조한 내륙부에 사막을 중심으로 뻗어있는 넓고 인구가 희박한 지역'이라고 하니, 호주의 독특한 자연환경을 가리키는 말들이 이것 말도고 더 있지 않을까?

Nullarbor라는 영단어에 대한 나의 추론은 정확하였다. Torsten Seemann의 설명을 인용하자면 다음과 같다.
The Nullarbor is a huge treeless plain that spans the area between south-west and south-east Australia. It comes from the Latin "nullus" (no) and "arbor" (tree), or "no trees". As this software will generate a tree, there is an element of Australian irony in the name.
ZGA의 장점은 일반적인 소비자용 PC(modern consumer PC)를 초과하는 시스템을 요구하지 않는다는 것, 그리고 long read를 다룰 수 있다는 점이다. 항생제 내성 유전자를 검색하는 기능은 없다. 유전자 주석화에는 역시 Seemann의 유명한 작품인 Prokka를 쓰는 것이 아니라 일본에서 개발한 DFAST(DDBJ Fast Annotation and Submission Tool)을 사용한다. 인용횟수가 8천회 가까이 되는 대중적인 Prokka를 택하지 않은 것은 이 응용 프로그램이 Perl을 기반으로 하기 때문이란다. Perl은 이제 버전 7을 향해 바뀌고 있어서(이것도 오늘 처음 알았음!) 앞으로 의존성 문제가 생길 수 있기에 아예 DFAST를 쓴 것이라 한다. BBTools를 적극적으로 사용하는 것도 바람직하다고 본다. Read 전처리를 위한 도구 역시 trimmomatic에서 fastp로 대체되었다.

다소 성가신 계산 작업을 하려고 머리를 쥐어뜯고 있다가 문득 BBTools 패키지에 이미 그러한 기능을 수행하는 스크립트가 들어 있음을 알게 되었을 때, 잠시 작업을 멈춘 다음 비로소 고개를 숙이고 숙연한 마음을 갖게 된다. MetaBAT 2를 이용하는 Best Binning Practices를 공부하다가 알게 된 jgi_summarize_bam_contig_depths 스크립트는 BBTools에 포함되어 있다.

오늘 제목으로 삼은 소프트웨어는 전부 bacterial isolate genome의 조립과 분석을 위한 것이다. 그러나 요즘 푹 빠져 있는 주제는 shotgun 방식의 metagenome assembly이다. 특히 Nanopore sequencing 자료 단독으로 어떻게 하면 좋은 분석물을 만들게 될지 눈이 빠지게 조사를 하고 있다. 일단 조립 측면에서는 metaFlye와 medaka가 무난하지만, circlator fixstart를 돌려 보면 여전히 dnaA gene은 시퀀싱 에러에 의해서 두세개로 조각이 난 상태이다. 결국은 일루미나로 추가 시퀀싱을 한 뒤 pilon을 돌려야 하는 것인지...

앗, 잠시 착각을 했었다. medaka 교정 전의 assembly를 가지고 fixstart를 하고 있었다. 교종 후의 것으로 circlator fixstart를 했더니 완벽하게 dnaA gene을 기준으로 시작 위치가 조정되었다.

Assembling microbial genomes from complex metagenomic samples using long nanopore reads. 원본 링크

ZGA 설치 후 테스트하기

Python 3.7을 써야 모든 dependency를 conda 환경에서 에러 없이 설치할 수 있다고 한다. 나는 다음의 명령을 이용하여 zga environment를 설치하였다.

$ conda create -n zga zga python=3.7

123설치 후 /opt/miniconda3/env/zga의 용량을 확인해 보니 3.1 GB 정도가 된다. TORMES는 8 GB의 Kraken DB를 설치해야 하고, nullarbor는 여기에 Centriguge DB(이것도 역시 8 GB)를 더해야 한다. 결과 리포트는 nullarbor가 더 '프로페셔널'해 보인다(샘플). ZGA는 가장 심플하고, TORMES는 그 중간 정도?

ZGA는 read quality check로부터 genome annotation에 이르는 6가지 과정 중 시작과 끝을 임의로 지정할 수 있어서 매우 편리하다. 예를 들어 이미 조립을 해 놓은 유전체 염기서열을 투입하여 CheckM에 의한 quality assessment가 가능하다.

오늘은 nullarbor까지 설치해 놓고서 세 가지 프로그램을 서로 비교해 가면서 용도에 맞추어 쓰는 연습을 해야 되겠다.

2022년 2월 7일 월요일

따옴표 안에 필드 분리자가 있는 텍스트 파일을 awk로 파싱하려면?

콤마로 구분된 필드로 구성된 텍스트 파일을 awk로 파싱하여 특정 컬럼만을 추출하고자 한다. 각 컬럼은 드는 다음과 같이 큰따옴표로 구분된 상태라고 가정하자. 

"column_1","column_2","column_3",...

그런데 하필이면 어떤 컬럼은 영문 이름 정보를 담고 있어서 그 안에 콤마가 또 들어가 있는 상태이다.

"Jeong, Haeyoung","other_column",...

단순히 'awk -F,' 명령을 이용해서는 이러한 상황을 잘 처리하지 못한다. 마이크로소프트 엑셀에서는 어떻게 대응하는지 아직 확인해 보지 않았다. 따옴표가 둘러싼 내부에서는 필드 구분자(field delimiter)와 동일한 문자를 만나더라도 건드리지 않게 하는 지혜가 필요하다. 이는 널리 논의된 문제이고 이미 해결 방안도 잘 알려져 있다. Starck Overflow 웹사이트에서도 중복된 질문과 답으로 표시되어 있기 때문이다.

    How to make awk ignore the field delimiter inside double quotes? [duplicate]

    예를 들어서 1번과 4번 컬럼만을 추출하려면 이렇게 하면 된다. 내가 처음에 참고하여 사용했던 방법이다.

    $ awk '{while(match($0,/("[^"]+",|[^,]*,|([^,]+$))/,a)){
         $0=substr($0,RSTART+RLENGTH);b[++x]=a[0]}
         print b[1] b[4];x=0}' sequences.csv > accession.csv.mod
    

    그런데 이 코드는 이해하기 난해하다. 답변 중에서 가장 상단에 있는 것이 가장 깔끔하고 보기에도 좋다. GNU awk의 '-vFPAT' 옵션을 쓰면 된다고 한다. 다음과 같이 명령어를 구성하는 같은 결과를 얻는다.

    $ awk -vFPAT='([^,]*)|("[^"]+")' -vOFS=, '{print $1,$4}' file
    

    콤마를 다른 문자로 바꾸면 다양한 형태의 텍스트 파일에 대응할 수 있다.

    2022년 2월 1일 화요일

    다이소 5천 원짜리 와이파이 수신기를 리눅스(CentOS 7.9)에 설치해 보다

    직장에서는 전산망보안 때문에 특정 웹사이트에서 파일을 다운로드하기 어려운 경우가 가끔 있다. 이를 해결하는 방법을 나름대로 연구하여 문서로 정리하여 놓았다. 그러나 443번 포트를 이용하여 파일을 다운로드하도록 파이썬 코드 안에 설정된 상태는 내 실력으로 해결하기가 어려웠다. 그래서 연구소에서 쓰는 것과 유사한 CentOS 환경을 집에 마련해 두고 부득이한 경우에 이용하고 있다.

    가끔 업무를 위해 꼭 접속을 해야 하는 웹사이트임에도 불구하고 국가정보원에서 공공기관으로 하여금 접속을 차단하도록 강제하는 IP address 목록에 등재되는 바람에 일을 못하는 경우도 있다. 그런 상황이라면 집에서 다운로드하여 가져와야만 한다. 근무지의 전산망 보안 담당자에게 구체적으로 문의하지 않고서는 알 수가 없다! 흑흑...

    이를 위해서 마련한 낡은 컴퓨터와 인터넷 공유기 사이의 거리가 멀어서 유선으로 연결하기가 매우 불편하였다. 다이소에 가면 매우 싼 가격의 와이파이 수신기가 있다고 하여 테스트를 해 보기로 했다. 5천 원짜리 제품이 두 가지가 있었는데, 드라이버 CD가 들어있지 않은 것으로 여겨지는 작은 패키지의 제품을 구입하였다. 겉면에 인쇄된 윈도우/리눅스/Mac OS X 지원이라는 말만 믿고서. 무선 마우스 연결을 위한 블루투스 동글 정도의 크기로 과연 제대로 수신이나 될런지?

    리눅스에서 자동 인식 또는 드라이버 설치가 될 수준의 물건이 아니다. 윈도우라면 또모르겠지만. 포장 뒷면을 보니 이런 문구가 있었다.

    드라이버 다운로드: http://www.freshfine.co.kr 접속 후 SUPPORT > FAQ

    웹사이트에 들어가서 와이파이 랜카드 150 Mbps 설치파일이라는 제목의 글을 찾은 뒤 다운로드 링크(.rar 파일)를 클릭하였다. 어이쿠, 일일 트래픽을 초과하여 전송이 불가능하단다. 웹사이트 관리자는 다음날이나 되어야 이를 해결할 것이 자명하다. 다른 방법을 알아보기로 하였다. 나중에 CentOS에서 쓸 수 있도록 컴파일된 rar/unrar 유틸리티(링크)를 설치하여 공식 웹사이트에서 제공하는 드라이버를 받아서 열어 보았더니 오직 윈도우용 드라이버만 들어 있었다.

    드라이버 파일 이름(RTL8188EUS ETV WIFI ADPATER.rar)을 참조하여 구글을 검색해 보았다. REALTEK 공식 웹사이트의 정보(링크)로도 리눅스 사용자에게 도움이 될 만한 것은 없었다. 최종적으로 GitHub에서 가장 적당해 보이는 것을 찾았다. 사용자 커뮤니티에서 자체적으로 리눅스용 드라이버를 만들어서 공개하는 것 같다.

    https://github.com/aircrack-ng/rtl8188eus

    README 파일을 살펴 보았다.

    Like https://github.com/cccooo/rtl8812au-centos-7.6, forked from aircrack-ng/rtl8188eus and modified for CentOS 7.9 as CentOS Kernel 3.10 contains many code from 4.x

    rtl8188eus v5.3.9

    Realtek rtl8188eus & rtl8188eu & rtl8188etv WiFi drivers

    현재 설치된 kernel은 3.10.0-1160.53.1.el7.x86_64이니 호환성에는 문제가 없을 것이다. 드라이버 소스를 git로 클로닝한 뒤 make를 했더니 에러가 발생하였다. kernel-devel 패키지가 없어서 발생한 문제였다. 필요한 패키지를 설치한 뒤 다시 'make && make install'을 실행한 뒤 재부팅을 하였다. 

    실제 컴파일 과정은 이보다 복잡하였다. 설치된 최신 커널의 버전과 같은 kernel-devel로는 make가 잘 되지 않아서 이전 버전(3.10.0-1160.45.1.el7)을 이용하니 모듈 빌드가 잘 되었다. 그런데 이를 쓰려면 부팅할 때에 커널을 이에 맞추어야 한다는 불편함이 있었다. 그래서 원래의 최신 커널로 부팅하여 다시 모듀을 빌드하니 이번에는 잘 진행이 되는 것이다. 재부팅을 몇 번이고 거듭하면서 최신 커널에 대하여 잘 작동하는 것을 확인하였다.

    그런데 와이파이 네트워크가 보이질 않는다. 다시 구글링을 해 보니 이더넷 케이블을 꽂은 상태에서는 안된다나? 이 힌트에 대한 원문 링크는 다시 찾기가 어렵다. 

    CentOS 7에 대한 지원은 2024년에 공식 종료된다고 한다. 내가 사용하는 리눅스 컴퓨터 중에 Ubuntu가 깔린 것이 점점 더 많아지는 것은 지극히 자연스럽다. 그렇다 해도 앞으로 CentOS Stream이 무엇인지 이해하는 것이 남은 숙제이다.

    [Red Hat Blog] FAQ: CentOS Stream Updates

    또는 Rocky Linux라는 것이 대안이 될 수도 있다.


    2022년 5월 17일 업데이트

    우분투 22.04 LTS 'Jammy Jellyfish'를 설치하였더니 이 와이파이 수신기가 잘 작동하였다. 사용할 장치는 다 꽂아 놓은 상태에서 리눅스를 설치하는 것은 상식이다.