2019년 9월 20일 금요일

YouTube Preminum을 취소하고 Mendeley를 업그레이드하다

약간의 비용을 지불하고 유튜브에서 광고 없이 음악을 듣는 것은 참 즐거운 일이었다. 돈을 내면, 그만큼 얻는 즐거움과 효용이 있다. 하지만 가끔씩 선택의 순간이 오고, 더 큰 효용을 위해 작은 것을 포기하는 용기를 발휘해야 할 필요도 있다.

2018년에 5월부터 매달 8,690원을 내기 시작했을 때의 서비스 이름은 YouTube Red였다. 구글 결제 센터에 접속하면 여전히 YouTube Red라는 서비스 명목으로 돈이 빠져나감을 알 수 있다. 그런데 해지를 하였더니 YouTube Premium 멤버십을 정녕 해지하실 거냐고 이메일이 왔다. 두 서비스가 사실상 같은 것이겠지? 잘 쓰던 서비스를 해지한다고 해서 별로 아쉬울 것은 없다. 음악은 인터넷 라디오로 들으면 되니까...

그대신 유튜브에서 Mendeley의 PDF 저장 공간(Personal Space라 불림)을 기본 2GB에서 5GB로 올리는데 비용을 지불하기로 했다. 파견 근무지(회사)에서는 개인 PC에 파일을 저장하는 것이 기본적으로 금지되어 있고, PDF에 메모나 하이라이트와 같은 수정을 가하면 lock이 걸려서 사용하기에 보통 불편한 것이 아니다. 차라리 Mendeley 웹에 논문 PDF를 올린 다음 읽는 것이 더 낫다. 업로드한 논문의 성격에 맞추어서 읽을만한 최신 논문을 제안하는 기능도 꽤 괜찮다. 어차피 여기에 올리는 자료는 공개적으로 얻을 수 있는 논문이므로 회사의 기밀 유지와는 관계가 없다고 생각하여 Mendeley를 적극 사용하기로 하였다. 흠, 내가 어떤 주제의 논문을 읽는다는 것이 기밀 사항이 될 수도 있겠군.

회사에서는 Gmail 웹사이트가 열리고, Mendeley 웹사이트의 내 계정으로 PDF를 올릴 수도 있다. DropBox가 되는지는 확인해 보지 않았다. 반면 KRIBB에서는 이상의 행위를 할 수 없다. 국정원의 정책 때문인 것으로 알고 있다. 회사든 KRIBB이든 일을 하기에는 비슷한 정도로 불편한 셈이다. 다행스럽게도 KRIBB에서는 파일에 보안 잠금은 되지 않는다. SSL 인증서와 관련해서 회사 전산망에서는 너무 많은 고생을 했고 지금도 해결 가능성은 보이지 않는다.

2012년에 Solar System이라는 패키지로 총 4회 결제를 했었다. 엊그제 새로 개시한 Plus package는 여전히 매월 4.99 달러를 내면 된다. 무려 7년이 지나서 다시 유료 서비스로 전환한 셈이다.

가급적이면 모니터 화면을 통해서 자료를 읽으려 하지만, 종이에 인쇄하여 줄을 쳐 가면서 읽는 아날로그적인 습관을 아직 버리지 못했다. 내가 전자책 단말기를 사기로 결정하는 날은 언제나 올지 모르겠다. 오늘도 퇴근 후 읽으려고 논문 하나를 종이에 인쇄하여 들고 나간다.

'Quick-and-Dirty' Perl 스크립트 한 줄의 작동원리 알아보기

여러 개의 염기서열을 담고 있는 sequence.fasta에서 특정 ID에 해당하는 서열만을 뽑아내는 한줄짜리 Perl 스크립트를 설명해 보고자 한다. 추출할 서열의 ID 목록은 ids.txt에 존재한다고 가정하자.

$ perl -ne 'if(/^>(\S+)/){$c=$i{$1}}$c?print:chomp;$i{$_}=1 if @ARGV' ids.txt sequences.fasta

이 예제는 Short command lines for manupulation of FASTQ and FASTA sequence files에서 가져온 것이고, 내가 작성하여 위키 페이지에 올렸던 문서에 포함되어 있다. 솔직하게 말하자면 교육 또는 학습 목적으로는 대단히 바람직하지 않은 코드이다. 길이가 매우 짧고 어쨌든 실행은 되니까 그저 복사해서 붙여놓고 아무런 의문을 제기하지 않고 써도 된다. 물론 나는 이보다는 조금 더 긴 Perl 스크립트를 별도로 만들어서 사용한다.

마침 이 코드의 정규표현식에 관한 질문이 들어와서 블로그에 상세하게 작동 원리 전체를 설명해 보기로 하였다. 겨우 한 줄에 불과하지만 이해하기는 쉽지 않다. 한참 다른 정보를 찾아서 고민하고 테스트를 거듭하여 겨우 감을 잡을 수 있었다.

먼제 Perl 인터프리터의 명령행 옵션인 -e와 -n을 알아야 한다. -e는 잘 알려진 것으로, 스크립트를 파일로 만들지 않고 명령행에서 직접 입력하여 사용하기 위한 것이다. -n은 -e 코드를 while loop로 둘러싸는 것이다.  Perl의 명령행 옵션에 관한 상세한 설명과 예제는 Perl Command-Line Options 또는 7 of the most useful Perl command line options를 참조하기 바란다.

그러면 본격적으로 이 고약한 코드를 해부해 보자. 코드는 다음과 같이 세 개의 블록으로 나눌 수 있다. 각각을 b-I, b-II, 그리고 b-III이라 부르기로 하자. perl -n 옵션을 주었으므로 ids.txt와 sequences.fasta의 모든 줄에 대해서 다음이 순차적으로 실행되는 implicit loop를 돈다.


ids.txt에 seq1, seq2(실제로는 한 줄에 한 ID이므로 "seq1\n", "seq2\n"라고 생각해야 됨)라는 내용이 담겨있다고 가정하자. 이 두 줄에 대해서 루프를 돌 때는 b-I은 건너뛴다. 왜냐하면 첫 줄이 '>'로 시작하지 않기 때문이다. 뒤이어서 b-II에서는 무슨 일이 벌어지는가? A? : command 1 : command 2는 if..the..else와 매우 흡사하다. '?'는 conditional operator라 부른다(참고). A라는 조건을 판별하여 참이면 command 1을, 거짓이면 command 2를 실행하라는 뜻이다. $c?라고 물어보았으니 $c가 1이면 print 함수를 , 그렇지 않으면 chomp 함수를 실행하게 되는데, 함수의 인수가 없다고 걱정할 것은 없다. 자동으로 $_ 변수(현재 읽어들여 처리 중인 라인)을 인수로 사용하기 때문이다. chomp 함수는 인수를 명시적으로 나타내지 않으면 $_의 뒤에 붙은 "\n"을 털어낸다. 이어서 이 코드의 세번째 블록(b-III)으로 진행되어 $i{seq1}=1가 선언되고, 다음번 줄에 대해서 반복을 돌면서 $i{seq2}=1가 실행된다. 만약 b-II에서 chomp를 실행하지 않았다면 $i의 key는 "seq1\n"이 되고 만다. 이는 "seq1"과는 다른 문자열로 간주됨에 유의하라. 결과적으로 ids.txt를 이루는 두 줄에 대해서 코드 전체(b-I/II/III)가 반복되고, %i hash에는 추출해야 할 서열의 ID가 key로서 기록된다.

ids.txt에 대해서 loop를 돌았으니 다음으로는 sequences.fasta 파일을 이루는 각 줄에 대해서 b-I/II/III을 반복하여 실행할 차례가 되었다. sequences.fasta의 첫 줄은 '>sequece_id' 형태이므로 block-I에서 패턴 매칭문을 만난다. /^>(\S+)/의 의미는 '>'로 시작하는 라인에서 '>'의 바로 다음에 오는 공백이 아닌 문자 여러개에 해당하는 패턴을 찾아서 $1 변수에 넣으라는 뜻이다. if {} 블록 내부의 $c=$i{$1}는 그러면 무슨 의미인가? 만약 지금 읽어들인 서열 ID가 ids.txt에 존재하는 것이라면 %i hash의 key에 정의된 상태가 된다. 따라서 $c=1이 된다. 만약 지금 막 읽어들인 서열 ID가 ids.txt에 없다면 $c는 아직 미정의 상태이다.

다음으로는 b-II로 진행한다. 조건연산자 '?'가 $c의 참 혹은 거짓 여부를 묻는다. $c가 정의된 상태라면 '>sequence_id'를 표준출력으로 인쇄한다. 정확히 말하자면 라인의 끝을 나타내는 "\n"까지 $_에 포함된 상태이니 정상적으로 한 줄이 다 인쇄되고 다음 라인을 인쇄할 준비를 마친다. 다음으로 b-III을 만난다. $i{">sequence_id\n"}=1이 실행되지만 이는 이후의 동작에 영향을 미치지 않는다.

이번에는 sequences.fasta의 두번째 줄을 실행할 순서이다. '>'로 시작하는 줄이 아니므로 b-I은 통과한다. b-II에서는 어떻게 될까? $c는 여전히 1의 값을 갖고 있으므로 print 문이 실행된다. b-III을 만나게 되지만 $i{"ATCGACG..\n"}=1을 실행하게 되고 역시 결과에는 영향을 미치지 않는다. 즉, 추출하고자 선택된 sequence의 ID 및 실제 염기서열 라인만 계속 표준출력으로 인쇄가 된다.

그러다가 ids.txt에 정의되지 않은 ID의 서열을 드디어 만나게 된다. '>'로 시작하는 라인을 만났으니 일단 b-I으로는 들어가는데, 이에 맞는 $i{서열ID}는 존재하지 않는다. 따라서 $c는 정의되지 않는다. 따라서 이어서 진행되는 b-II에서는 chomp를 만나게 된다. 이는 마찬가지로 출력물에 영향을 미치지 않는다. 이 동작은 다시 '>'로 시작하는 줄을 만나게 될 때까지 반복된다.

솔직하게 말해서 나의 코드 분석이 100% 맞는다고 자신은 못하겠다. 겨우 한 줄에 불과한 코드를 이렇게 긴 글로 설명해야 한다면, Perl 스크립트 작성 학습용 샘플로는 적합하지 않다. 물론 사용하는 데에는 아무런 문제가 없다. 1996년부터 2000년까지 소위 '난해한 Perl 코드 콘테스트(Obfuscated Perl Contest)'가 있었다고 하니 500바이트 내외의 Perl 코드가 불러일으키는 난잡함, 그러나 그 속에서 풍겨나는 익살과 아름다움을 감상하고 싶다면 위키피디아 혹은 다음의 수상작 페이지를 방문해 보라.


2019년 9월 17일 화요일

독서 역주행 - <두 얼굴의 조선사(2016)>와 <공자가 죽어야 나라가 산다(1999)> 이어서 읽기

나의 성은 나주 정(丁)씨 '○○공(公)'파이다.

"저는 ○○ 김가입니다." ...(1)
"저는 ○○ 박씨입니다." ...(2)

다른 사람 앞에서 말을 할 때에는 (1)처럼 하는 것이 예의에 맞다고 한다. 하지만 이 블로그에서처럼 평어체로 쓰는 글에서는 (2)처럼 해도 큰 문제가 없을 것 같다. '○○공'파라고 표현한 것은 이 정보가 정확하지 않을 수도 있기 때문이다. 집안 어른께서 우리 정씨의 내력을 말씀해 주신 것을 제대로 기억하지 못해서 그런 것이 아니다. 우리 집안의 ○○공 파보(派譜)를 정리한 것은 돌아가신 할아버지의 노력이었다. 그런데 이것을 나주 정씨 '메인 스트림'에 해당하는 월헌공파에서 인정을 하지 않는 듯한 느낌을 받은 일이 있다. 그것도 처음 만나는 사람이 당사자인 내 앞에서 그런 이야기를 한다는 것은 매우 무례한 일이 아니겠는가?

할아버지께서는 고향(경기도 어느 지역 A라 하자)에 남아있는 우리 조상의 비문, 우리 성씨와 연결된 지명 및 유적의 이름, 그리고 최근 기록만 남은 단편적인 우리 집안의 족보를 가지고 쥬류를 이루는 정씨의 족보에 힘겹게 연결하는 작업을 하셨노라고 내게 말씀하신 적이 있다. 내가 알기로 '○○공파'가 만들어진 것도 할아버지의 노력인 것으로 안다. 그렇다면 혹시 원래 우리 집안은 정씨가 아닐 수도 있는 것 아닌가? 조선 후기가 되면서 양반이 점점 많아지고, 전부 성과 족보를 갖게 되었다. 양반이 갑자기 생식능력이 왕성해져서 그 수가 늘었나? 절대 그럴리는 없다. 대부분은 족보를 새로 만들면서 권문세가의 성을 슬쩍한 것이 다반사일 것이다. 게대가 태반은 중국에서 들어온 사람이 자기 성씨의 시조가 되었다고 하니 정말 웃기는 노릇이다. 하긴 나주 정씨도 정덕성이란 사람이 중국에서 한국 압해도에 유배를 와서 시작했다고 전하지만, 다산 정약용은 이를 믿을 수 없다고 하였다. 보다 확실한 중시조는 고려 때의 무인 직책인 검교대장군을 지낸 정윤종(정덕성의 15세손이라고도 함)이다.

우리 조상이 원래 근본이 없는 집안을 이루고 살다가 수백년 전에 A 지역에 들어와서 丁씨 행세를 한 것일까? 기왕 그럴 것이라면 몰락한 왕족이나 당시 세도를 누리던 집안의 성을 따르지 않았을까? 아니면 티를 덜 내기 위하여? 지금으로선 우리 집안이 과거에 양반이 아니었다 하더라도 상관이 없다는 생각이다.

명절을 맞이하여 우리 나라 인구를 차지하는 성씨가 너무가 다양성이 부족하고, 대부분이 왕족이나 권문세가의 후손임을 자처하는 것이 한심하다는 생각이 들었다. 그렇게 된 원인은 어디 있을까 고민을 하다가 방송 다큐멘터리 작가 조윤민이 지은 '두 얼굴의 조선사(2016)'를 구입해서 읽게 되었고, 지금으로부터 꼭 20년 전인 1999년에 상명대 김경일 교수가 지은 '공자가 죽어야 나라가 산다'를 중고책 서점에서 구해서 단숨에 읽었다.


중국땅을 거쳐간 왕조 중에는 한나라 이래로 300년을 넘긴 것이 없다. 그러나 조선은 대한제국 선포 전까지만 따져도 500년이 넘는 역사를 자랑한다. 그런데 이것이 길다고 해서 과연 자랑할마한 일일까? 그 기나긴 역사 속에서 우리는 자주적인 근대화의 싹을 틔워 왔는가? 결코 그렇지 않다는 것이 왼쪽 책의 결론이고, 나는 그것의 근본 이유가 유학에 있다고 생각하여 오른쪽 책까지 읽게 된 것이다. 

조선의 지배세력은 그야말로 '정치 공학자'라고 생각한다. 유교라는 이념화된 사상이야말로 가장 강력한 무기였다. 흔히 유교의 기본 정신 중에 仁을 가장 중요한 것으로 생각하지만, 중국 고대문자 및 문헌 전문가인 김경일 교수에 의하면 공자가 살던 시기에 仁이라는 글자는 없었다고 한다. 오래 전 중국에서 적당히 다듬어진 사상이 어쩌다가 조선의 통치 이념이 되고, 그 틀에 갖힌 기득권 세력은 우리의 창의롭고 진취적인 기운을 죄다 틀어막은 것만 같다. 과거의 기준에 맞지 않으면 사문난적으로 몰려서 사약을 한 사발 죽 들으켜야 하는 사회에서 어찌 발전으로 생각하고 미래를 대비할 수 있단 말인가?

왼쪽 책은 조선왕조실록 등 근거를 제시해 가면서 학문적으로 기술하였고, 오른쪽 책은 일종의 에세이에 해당한다. 오른쪽 책은 당시 IMF 구제금융 시대를 맞아 답답한 현실을 개탄하면서 쓰여진 책이다. 저자는 책을 내고 나서 유림들에게 고소를 당하는 등 많은 어려움을 겪었다. 하지만 원래 저자가 의도한 제목은 훨씬 강경한 '공자를 죽여야 나라가 산다'였다나? 두 책을 읽는 내내 머리에서 김이 뿜어져 나오는 것만 같았다. 

조선에 대한 비판적인 시각을 갖고 있다고 하여 이것이 식민사관이라든가 요즘 물의를 일으키는 책 '반일 종족주의'를 옹호하고 있는 것으로 오해하지는 말라. 

다음에 읽을 책은 김경일 교수의 2013년도 저작 '유교 탄생의 비밀'이 될 것이다.

2019년 9월 16일 월요일

화성으로 가는 탑승권

NASA에서 2020년 발사 예정인 화성탐사선(MARS2020)을 기념하여 탑승권을 발급하는 이벤트를 벌이고 있다. 마감은 9월 말까지이다. 탑승권을 받은 사람들의 이름은 칩에 담아 우주선에 실린다고 한다. 사람을 보낼 수는 없는 노릇이니.

나는 9,342,924번째 탑승자이다. 탑승권의 원본 사이트는 여기에 있다.


이 이벤트에 관한 소식은 대덕넷에서 접했다. 재치있고 낭만적인 이벤트가 아닌가. 여기에 접속하여 이름과 국가명, 우편번호, 그리고 이메일 주소만 넣으면 된다.

국내만 '20만명'..."이번달까지 화성에 이름 보내세요"

잠시 상상의 나래를 펼쳐보자. 인적사항을 담은 화성 탐사선이 지구를 출발하여 가는 도중에(혹은 화성에 도착하여) 외계인에게 납치가 된다. 외계인은 탐사선이 보유한 모든 정보를 해독하여 지구인 백만명의 이메일 주소를 알아내고(그렇다, 외계인에 의한 해킹이다!) 악성코드를 심은 이메일 메시지를 뿌린다. 이메일을 열어본 사람은 갑자기 정신상태가 이상해지면서 외계인들에게 조종당한다... 어, 잘만 다듬으면 공상과학영화 시나리오 한 편이 나오겠는걸.

실제로 우주선에 실리는 정보는 이름 정도일 것으로 생각된다. 덜렁 영문으로 쓰인 이름만 가지고 엉뚱한 일을 하기는 어려울 것이다.

이러한 이벤트를 우리나라에서 할 수 있을까? 아마도 개인정보보호법에 의거하여 동의서를 받아야만 될 것이다. 우리나라에서 이메일 주소는 개인정보로 취급하고 있다. 다른 정보를 포함하면 개인을 식별할 수 있다고 판단하기 때문이다. 따라서 화성 탐사선에 내 이름을 실어보내자는 낭만적인 이벤트를 우리나라에서는 하기 어렵다. 자유로운 상상력을 틀 안에 가두는 규제가 우리에게는 너무 많다.

2019년 9월 12일 목요일

6LQ8 PP 앰프의 소소한 개선 - low-pass filter를 달기

6LQ8 진공관은 원래 오디오용 주파수의 증폭에 쓰기 위해서 개발된 것이 아니다. 가뜩이나 증폭률이 높은 진공관이라서 오디오용 앰프를 만들면 발진(oscillation)을 일으키기 쉽다고 한다. 발진이라고 해서 반드시 '삐~'하는 소리가 난다는 뜻은 아니다. 저항이나 캐패시터와 같은 수동 부품을 적당히 덧입혀서 발진을 잡는 방법이 잘 알려져 있지만, 어차피 오디오용 관이 아닌 것으로 진지하게 음악 감상을 하기 위한 앰프를 만드는 것을 싫어하는 사람들도 있다. 이것은 선택과 취향의 문제이다.

다음은 제이앨범에서 권장하는 6LQ8 푸시풀 앰프용의 '검증된' 발진 대책이다. 초단(전압증폭단)의 플레이트에 연결된 부하 저항(R5)에 저항+캐패시터로 이루어진 바이패스 회로를 붙이는 것이다. 이어지는 회로에는 낮은 주파수를 보내는 셈이니 일종의 low-pass filter라고 불러도 무방하겠다. 높은 주파수를 차단하는 역효과가 없지는 않겠으나 가청 주파수 영역에서는 상관이 없을 것이다. 이렇게까지 하여 오디오용이 아닌 진공관을 사용하여 오디오 앰프를 만드는 것에 대하여 부정적인 시각을 가진 사람이 없지는 않다. 이런 노력을 들일 바에야 아직도 쓸만하게 많이 남아있는 오디오용 진공관을 쓰는 것이 낫다는 의미가 되겠다.


부품을 붙일 자리가 없어서 PCB를 뒤집어서 패턴면에 부품을 납땜하였다. 중간쯤에 보이는 103 메탈필름캐패시터(IC114 링크)가 이번에 새로 붙인 부품이다. 한쪽 리드선에는 1K 저항을 직렬로 연결하였다.



신호를 연결하지 않고 볼륨 노브를 최대로 돌렸을 때 들리는 잡음은 확실이 줄어들었다. 그러면 가청 주파수의 높은 쪽 대역에서 손해를 보았을까? 파형 발생기와 오실로스코프로 측정을 하지 않으면 알기가 어렵다. 뿐만 아니라 바이패스 회로를 연결하기 전후의 상황을 비교해야 차이를 발견할 수 있을 것이다. 단지 음악 신호를 연결해서는 고음 신호를 깎아먹어서 음질이 더 나빠졌는지를 알 수가 없었다. 약간 둔탁해진 것인가? '그럴 것'이라는 선입견이 감각을 왜곡하는지도 모르겠다. 만약 측정을 통해서 가청 주파수 대에서는 변화가 없다는 것이 입증되면, 아마도 그 순간부터는 갑지가 소리가 좋게 느껴질지도 모르는 것이다. 마치 플라시보 효과처럼 말이다.

2019년 9월 13일 업데이트

필터를 붙인 이후 소리가 어딘가 모르게 둔탁해졌다는 느낌이 들기 시작했다. 고음을 너무 깎아버려서 그런 것이 아닐까. 앰프를 다시 작업대에 올리고 한쪽 채널만 필터를 제거하고 신호를 연결하여 좌우를 비교하니 너무나 확연한 차이가 났다. 이런 식의 잡음 제거 효과라면 없느니만 못하다는 생각에 양 채널에서 필터를 전부 제거하였다. 정상적인 negative feedback 회로 이외의 것을 달아서 잡음 등을 제거하려는 시도는 상당히 조심해야 되겠다.

2019년 9월 9일 월요일

Ubuntu에서 hear string이 작동하지 않는 이유

CentOS(bash)에서는 잘 돌던 스크립트가 Ubuntu에서는 에러가 발생하였다. 왜 그런 것일까?

$ cat test
echo $(cut -d_ -f1 <<< "123_456")
$ sh test
test: 1: test: Syntax error: redirection unexpected
$ bash test
123

이유를 알아보자. 우분투의 기본 shell은 bash가 아니고 dash라고 한다. 이 환경에서는 here document 또는 here string이 잘 작동하지 않는다고 하였다.

[Ubuntu/Linux] #!/bin/sh에 대한 간단한 이야기

따라서 위에서 소개한 test 스크립트의 shebang line에 #!/bin/sh를 삽입해 봐야 소용이 없다. 우분투에서는 /bin/sh가 dash를 가리키기 때문이다. 이 한줄짜리 스크립트를 제대로 돌리려면 첫줄에 #!/bin/bash를 삽입해야 한다. 그러나 실행할 때 sh test라고 하면 도루묵이다. 실행 권한을 주어서 ./test라고 실행해야 한다. 당연한 이야기이다.


AWK를 이용한 fasta file splitter

여러 염기서열로 이루어진 multi-fasta file을 쪼개에서 각 서열이 하나씩의 파일에 담기게 하는 작업은 매우 흔히 벌어지는 일이다. 나도 10줄 조금 넘는 Perl 스크립트를 만들어서 fasta file을 분리하는데 종종 쓰고는 한다.

awk를 쓰면 달랑 코드 한 줄로 fasta file splitter를 만들 수 있지 않을까? 구글에서 검색을 해 보니 이에 대한 많은 질문과 답이 있었다. 이를 참조하여 다음과 같은 awk 명령어를 만들어 보았다.  awk의 명령에 내에서 shell과 유사하게 파일로 직접 인쇄결과를 리다이렉션할 수 있다는 것을 몰랐다. 파일핸들을 정의하고 여닫는 일을 하지 않아도 되니 Perl보다 더 간단한 코드가 되었다. 아무런 실수가 없을 것으로...믿는다.

$ awk '/^>/{s=$0; sub(/^>/, "", $1); F=$1".fa"; print s > F; next}{print >> F}' input.fasta

awk를 쓰면서 종종 혼동하는 것은 sub() 함수가 치환된 문자열을 반환하지는 않는다는 것이다. 즉 var = sub(regexp, replacement, target) 형태로 사용하면 안된다.

코드를 간단히 분석해 보자. 서열 ID(>seq_id)가 있는 줄을 만나면, 첫번째 컬럼($0)에서 '>'를 제거한 뒤 .fa를 붙여서 출력할 파일의 이름을 만들어 변수 F에 저장한다. 그런데 sub() 함수에서 직접 $1을 건드리게 되므로, 나중에 원본 라인을 그대로 출력하기 위하여 명령어 맨 처음에 s라는 변수에 미리 서열 ID 라인을 저장해 놓는 꼼수(s=$0)를 사용했다. 그 다음에 유의할 것은 '/조건/{조건에 맞는 줄에 대하여; next} {조건에 맞지 않는 나머지 줄에 대해 실행}' 명령을 이해하는 것이다. 여기에서 next가 매우 중요한 역할을 담당한다. 이에 대해서는 예전의 포스팅에서 두번 정도 다룬 것 같다.

이 방법을 이용하면 SRR2162225(Gopalakrishnan et al., 2018)의 fastq 파일을 각 샘플별로 쪼갤 수도 있다. 이 fastq file은 다음과 같은 형식의 sequence ID를 갖고 있다. 내용상으로는 이미 demultiplexing이 된 것에 해당하지만 하나의 fastq 파일에 전부 수록되어 있다.

@ERR2162225.10 Wargo.109865_330283 length=253

두번째 필드를 추출하여 빨간색으로 표시된 것을 파일 이름으로 하는 fastq 파일을 만들어 보자. 현 디렉토리에 43개의 fastq 파일이 생겼다.

$ awk '/^@ERR2162225/{s=$0; sub(/_.+$/, "", $2); F=$2".fastq"; print s > F; next}{print >> F}' ERR2162225.fastq
$ ls
ERR2162225.fastq    Wargo.121296.fastq  Wargo.122637.fastq  Wargo.124418.fastq  Wargo.130511.fastq  Wargo.132266.fastq
Wargo.109865.fastq  Wargo.121409.fastq  Wargo.122639.fastq  Wargo.124419.fastq  Wargo.130801.fastq  Wargo.132268.fastq
Wargo.115690.fastq  Wargo.121585.fastq  Wargo.122724.fastq  Wargo.124421.fastq  Wargo.131051.fastq  Wargo.132538.fastq
Wargo.116774.fastq  Wargo.121891.fastq  Wargo.123009.fastq  Wargo.126747.fastq  Wargo.131300.fastq  Wargo.133270.fastq
Wargo.117349.fastq  Wargo.121892.fastq  Wargo.123197.fastq  Wargo.126802.fastq  Wargo.131302.fastq
Wargo.118369.fastq  Wargo.122114.fastq  Wargo.123305.fastq  Wargo.127261.fastq  Wargo.131512.fastq
Wargo.118370.fastq  Wargo.122410.fastq  Wargo.123704.fastq  Wargo.127262.fastq  Wargo.131541.fastq
Wargo.120662.fastq  Wargo.122412.fastq  Wargo.123986.fastq  Wargo.130510.fastq  Wargo.132265.fastq

다른 사례를 하나 더 들어보자. 다음은 vsearch --derep_fulllength로 dereplication을 거친 염기서열 파일의 일부를 보이고 있다.

>Mock.1;size=492
TACGTAGGTGGCAAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGTAGGCGGTTTTTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGAAAACTTGAGTGCAGAAGAGGA
AAGTGGAATTCCATGTGTAGCGGTGAAATGCGCAGAGATATGGAGGAACACCAGTGGCGAAGGCGACTTTCTGGTCTGTAACTGACGCTGATGTGCGAAAGCGTGGGGATCAAACAGG
>Mock.2;size=345
TACAGAGGGTGCGAGCGTTAATCGGATTTACTGGGCGTAAAGCGTGCGTAGGCGGCTTATTAAGTCGGATGTGAAATCCCCGAGCTTAACTTGGGAATTGCATTCGATACTGGTGAGCTAGAGTATGGGAGAGGA
TGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGATGGCGAAGGCAGCCATCTGGCCTAATACTGACGCTGAGGTACGAAAGCATGGGGAGCAAACAGG

size=뒤에 나오는 모든 숫자를 더하면 dereplication 작업 전 FASTA file에 포함된 전체 서열의 수를 알 수 있다. 오늘 익힌 awk의 sub() 함수를 이용하면 간단하다.

$ awk 'BEGIN{s=0}/^>/{sub(/^>.*size=/, "", $0);s+=$0}END{print s}' input.fasta
4268

별로 설명할 필요도 없다. 서열 ID 라인에 description이 있다면 $0을 $1로 바꾸면 그만이다. 아니, 처음부터 code를 $1로 고쳐 놓으면 description이 존재하는지 여부에 상관없이 잘 돌아간다.

Bioinformatics on-liner

One-liner 애호가를 위하여 Stephen Turner의 GitHub 사이트를 소개한다(링크).

2019년 9월 6일 금요일

Robert Edgar, 멋진 사업가? 과학자? 프로그래머? 소프트웨어 엔지니어?

16S rRNA amplicon sequencing에 기반하는 메타게놈 분석 방법(주로 QIIME)을 요즘 집중적으로 공부하는 과정에서 QIIME documentation 웹사이트 말고도 하루에도 몇번씩 방문하는 곳이 있다. 바로 http://drive.com/이다. 이 사이트의 운영자는 Robert Edgar이다.

Robert Edgar(출처)
1982년 런던 대학에서 입자물리학으로 박사학위를 받았고 IT 분야의 비즈니스 경험도 많지만 생명정보학 분야에서는 MUSCLE이나 USEARCH와 같은 유명한 프로그램의 개발자로 잘 알려져 있다. 그는 2001년부터 스스로를 특별한 소속 없이 일하는 독립과학자(independent scientist)로 칭하면서 주로 전산생물학과 생물정보학 분야에서 많은 사람들에게 도움을 주고 있다.

그가 어떤 사람인지는 개인적으로 전혀 모르지만, 대학이나 기업에 몸담고 있지 않으면서 이렇게 프리랜서와 같이 일하면서 꾸준히 프로그램을 개발하여 배포하고 논문을 발표한다는 것이 대단하다고 느껴진다. 그의 홈페이지를 방문하면 최근 발표한 논문의 목록이 보이는데 저자는 전부 혼자다! 1975년에 대학에 입학하여 지금 교수로 재직한 전형적인 한국인 과학자가 있다고 하자. 이제 정년 퇴임이 얼마 남지 않은 나이일 것이다. 과연 동료나 후배(선배는 전부 은퇴했을 것이니), 또는 대학원생의 참여 없이 논문에 저자로 이름을 올릴 수가 있을까?

R. Edgar의 USEARCH는 64비트 버전을 이용하려면 돈을 주고 구입해야 한다. 하지만 32비트 버전은 무료이다. 무료라고 해서 소스 코드가 공개되어 있지는 않다. 프로그램의 기본 알고리즘은 논문이나 drive5.com 웹사이트에 소개가 되어 있으나 소스 코드는 비공개이니 사람들이 불만을 느낄 수도 있다. 그런 취지에서 VSEARCH가 3년 전에 나왔고, 웹을 검색해 보니 아주 훌륭한 16S rRNA 서열 분석용 파이프라인이 공개된 것도 발견할 수 있었다(링크). 이것 역시 공부하기에 아주 좋은 자료이다.

VSEARCH의 논문 초록에서는 USEARCH를 살짝 '디스'하고 있다. 인용해 볼까?

 [VSEARCH] is designed as an alternative to the widely used USEARCH tool () for which the source code is not publicly available, algorithm details are only rudimentarily described, and only a memory-confined 32-bit version is freely available for academic use.

그래도 USEARCH는 학계에서 충분히 제 역할을 하고 있지 않은가?

나는 생명정보학을 체계적으로 배우질 않아서 이 학문의 역사에 대해서는 매우 단편적인 지식밖에 갖고 있지 못하다. 나의 짧은 경험으로 존경하고 싶은 생명정보학자, 또는 이 분야의 소프트웨어 개발자는 이런 분들이 있다. 이 목록은 절대로 생명정보학 명예의 전당 같은 것도 아니고, 이 분야를 모두 망라한 것도 아니다. 누구나 흔히 생각할 수 있는 David Lipman 같은 사람은 뺐다. 내가 현재 하고 있는 일에 대해서 직접적인 영향을 미친 사람들이다.

Margaret Oakley Dayhoff
Phil Green
David Gordon(CONSED의 개발자)
Jim Kent(흠, 정작 나는 BLAT를 거의 쓰지 않는데)
Sean Eddy
Lincoln Stein(BioPerl하면 떠오르는 사람)
Art Delcher(MUMmer 개발자)
Ewan Birney
Eugene Myers 등.

여기에 R. Edgar도 넣어줘야 되겠다.

개발자, 프로그래머, 컴퓨터 사이언티스트, 소프트웨어 엔지니어의 차이는?

다음의 글을 참조하라.


그리고 이 글에서 인용한 원문 및 그림(링크).

출처: Scott Hasnselman

2019년 9월 4일 수요일

QIIME 1을 위한 sequence data 준비 과정 - demultiplexed paried-end fastq의 복잡한 사정

약 일주일에 걸쳐서 QIIME 1을 집중적으로 공부하면서 느낀 것은 sequencing data 파일로부터 seqs.fna를 만드는 과정이 전체 수고의 1/3 정도는 되는 것 같다는 점이다. 더군다나 아직 실제 16S rRNA 데이터를 생산한 경험이 없어서 SRA에서 자료를 받아서 연습을 하니 그 과정이 더욱 어렵게 느껴졌다. 왜냐하면 서열 ID 라인의 형식이 실제 MiSeq에서 생산된 직후의 fastq 파일 모습을 그대로 따르고 있지 않으며, 메타데이터 역시 제대로 마련되지 않았기 때문이다. 하지만 그 덕분에 SRA-Toolkit의 prefetch와 fastq-dump의 기능을 제대로 이해하게 되었으며, NCBI 웹사이트에서 run accession의 목록 파일을 손쉽게 입수하는 방법을 알게 되는 성과는 있었다.

QIIME 1에서는  Roche/454 pyrosequencing data 또는 barcode read가 별도로 주어지는 Illumina multiplexed single end data가 가장 다루기 편한 것 같다. 바코드 서열과 샘플의 연결 관계는 metadata mapping file(map.tsv)에서 주어진다. 이는 간단히 mapping file이라고도 부른다. 그러나 2019년 9월 초순 현재 국내 시퀀싱 서비스 업체에서는 16S rRNA 시퀀싱 데이터를 이러한 형태로 고객에게 전해주는 것 같지는 않다. 결과 파일인 seqs.fna에서 각 염기서열의 ID는 (샘플이름)_(일련번호)의 형식을 갖는다.


$ split_libraries_fastq.py -o slout/ -i forward_reads.fastq.gz -b barcodes.fastq.gz -m map.tsv

Paired end read라면 join_paired_ends.py 스크립트로 연결한 뒤 추가적으로 바코드 정보를 업데이트해야 하므로 살짝 작업이 번거로워진다. 바코드 업데이트가 필요한 이유는 모든 read pair가 전부 성공적으로 연결되는 것은 아니기 때문이다. 지난주에 천랩에 고객에게 제공되는 16S rRNA amplicon sequencing data의 타입에 대해 문의하니 MiSeq 한 레인 단위의 복수의 샘플을 받아서 시퀀싱을 한 다음 전부 demultiplexing을 해서 고객에게 보내준다고 하였다. 이러한 데이터는 QIIME 1에서 오히려 다루기가 더 불편하다. QIIME 2의 문법으로 말하자면 이는 Casava 1.8 paired-end demultiplexed fastq 파일에 해당한다.

Demultiplexed paired end 데이터를 이용한 QIIME 분석의 좋은 사례는 2016년도 EDAMAME 워크샵(위키GitHub)에서 제공한 교육 자료에 있다. 이 자료를 만든 사람들 중에는 한국 사람이 틀림없는 이름(Sang-Hoon Lee)이 있어 반갑게 느껴지기도 했다. 여기에서 샘플 파일을 다운로드하여 압축을 풀면 V4 region을 증폭하여 MiSeq으로 양쪽에서 150 bp를 읽은 파일 54쌍이 나온다(C01D01F_sub.fastq, C01D01R_sub.fastq...). QIIME 1 script인 join_paired_ends.py는 한번에 한 쌍의 파일만 처리하므로, 파일이름 베이스를 수록한 list.txt를 만들어서 각 파일에 대하여 순차적으로 join_paired_ends.py을 실행한 다음 연결이 된 sequence만 가져다가 사용하는 스크립트인  Merged_Reads_Script.sh를 사용하는 방법을 소개하였다. 이 스크립트를 실행하면 몇 단계를 거쳐서 Merged_Reads라는 디렉토리에 52개의 fasta 파일이 만들어진다. 다음 단계에서는 add_qiime_labels.py를 사용하여 combined_seqs.fna라는 최종 결과물을 만든다. add_qiime_labels.py는 mapping file을 필요로 한다. 원래 이 파일에는 BarcodeSequence 컬럼이 있어서 각 read에 부가된 바코드로부터 sequence를 구분하는 중요한 역할을 하지만, 이 경우에는 이미 demultiplexing이 되어 있어서 이는 중요하지 않다. 단 InputFileName 컬럼에 있는 fasta file 이름이 각 sequence 파일을 분류하는 기준이 된다. 이 경우에는 필수 컬럼인 BarcodeSequence 컬럼은 공백으로 비워져 있어도 된다.

나는 Merged_Reads_Script.sh를 응용한 merge_paired_fasta.sh 스크립트를 만들어 보았다. 별도의 list.txt 파일이 없이도 쌍을 이룬 fastq 파일이 있으면 작동하며, 최종 결과물로서 fasta + fastq를 같이 만들어내는 것이다. Merged_Reads_Script.sh에서는 merged fastq를 지우도록 만들어 놓았다.

combined_seqs.fna 또는 seqs.fna 파일 내에서 각 서열의 번호는 새로운 샘플에 대해서 무조건 0또는 1부터 시작한다고 생각했는데 그건 착각이었다. 위에서 설명한 방법으로 만들어진 combined_seqs.fna의 서열 일련번호는 _1에서 시작하여(_0으로 시작하게 만들 수도 있음; 전산인들이 좋아하는 방법) 샘플이 바뀌어도 계속 증가한다. 사실 각 read가 구별만 되면 충분한 것이다. 만약 샘플마다 1번부터 시작하게 만들면 일련번호만으로 sequence를 구별하기 어려워지므로 이는 당연한 것이다.

Demultiplexed paird end data를 다루는 방법은 이상의 것이 가장 권장할 만한 것으로 보인다. 이것 외에도 오늘 오후를 종일 투자하여 QIIME 스크립트 자체에 충실한 방법을 만들어내기는 했는데 하도 신경쓸 것이 많아서 과연 이것을 블로그에 기록하는 것이 현명한 일인지 아닌지 잘 판단이 서지 않는다. 간단히 소개하자면 여기에서는 multiple_join_paired_ends.pymultiple_split_libraries.py를 사용한다. 하나의 디렉토리 안에 샘플별로 이름이 붙은 read pair가 전부 들어있는 경우(아마 대부분의 경우일 것임)와, 샘플별로 구분된 디렉토리 안에 한 쌍씩의 파일이 들어있고 이것 전체가 하나의 인풋 디렉토리에 들어있는 경우에 따라서 옵션이 약간 달라진다. 후자의 경우라면 fwd 및 rev read fastq file이 같은 이름을 갖고 있을 수도 있다. 샘플별로 이름이 붙은 디렉토리에 들어있으니 당장은 문제가 되지 않는다.

EDAMAME 데이터를 예로 들어서 설명해 본다. 54쌍의 read는 전부 샘플 이름을 기준으로 이름이 지어졌으니 하나의 인풋 디렉토리(예: Fastq)에 넣어서 작업을 해도 문제가 되지 않는다. --include_input_dir_path, --include_input_dir_path --remove_filepath_in_name 옵션을 주고 실행한 다음 결과 디렉토리가 어떻게 다르게 만들어지는지 한번 확인해 보라. --remove_filepath_in_name는 단독으로 쓰이지 못하고 항상 --include_input_dir_path와 같이 쓰여야 한다.

$ multiple_join_paired_ends.py -i Fastq --read1_indicator F_sub --read2_indicator R_sub -o output_folder
$ ls output_folder
C01D01F_sub  C01D02F_sub  C01D03F_sub  C02D01F_sub  C02D02F_sub  log_20190904164417.txt
$ ls output_folder/C01D01F_sub/
fastqjoin.join.fastq  fastqjoin.un1.fastq  fastqjoin.un2.fastq

각 샘플에 대하여 결과물에 대한 서브디렉토리가 생겼고 그 하위에 연결(join)된 것과 그렇지 않은 것들에 대한 fastq file이 생겼다. 여기서부터 주의를 해야 한다. 연결되지 않은 파일은 일단 제거를 해야 하고, 결과 디렉토리며에서 _sub는 없애는 것이 여러모로 편하다. 이것은 알아서 해결을 했다고 가정하고, 다음 단계인 multple_split_libraries_fastq.py로 넘어가자.

$ multiple_split_libraries_fastq.py -i output_folder/ -m sampleid_by_file -o final 
$ head -n 1 final/seqs.fna 
>fastqjoin.join.fastq_0 HWI-M03127:41:ACE13:1:2109:11596:14331 1:N:0:GGAGACAAGGGA orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0

$ multiple_split_libraries_fastq.py -i output_folder/ -m sampleid_by_file -o final2 --include_input_dir_path
$ head -n 1 final2/seqs.fna
>C01D01fastqjoin.join.fastq_0 HWI-M03127:41:ACE13:1:2109:11596:14331 1:N:0:GGAGACAAGGGA orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0

$ multiple_split_libraries_fastq.py -i output_folder/ -m sampleid_by_file -o final3 --include_input_dir_path --remove_filepath_in_name
$ head -n 1 final3/seqs.fna
>C01D01_0 HWI-M03127:41:ACE13:1:2109:11596:14331 1:N:0:GGAGACAAGGGA orig_bc=AAAAAAAAAAAA new_bc=AAAAAAAAAAAA bc_diffs=0

--include_input_dir_path --remove_filepath_in_name 옵션을 전부 준 세번째 것이 가장 바람직한 서열 ID를 만들어 내었음을 알 수 있다. 이것으로 pick_open_reference_otus.py를 아무런 에러 없이 실행하였다.

split_libraries_fastq.py를 무식하게 실행하여 소기의 목적을 달성할 수도 있다. 다음은 QIIME의 스크립트 설명에 나오는 가장 마지막 예제이다.

$ split_libraries_fastq.py -i lane1_read1.fastq.gz,lane2_read1.fastq.gz --sample_ids my.sample.1,my.sample.2 -o slout_not_multiplexed_q20/ -q 19 --barcode_type 'not-barcoded'

Fastq 형태로 join을 마친 파일을 한 곳에 가져다 놓고 그 이름을 콤마로 연결하여 -i file1,file2,file3...으로, 그리고 순서를 맞춘 샘플 이름을 --sample_ids sample1,sample2,sample3...으로 제공하면 될 것이다. Shell script를 이용하여 이는 어렵지 않게 구현해 보았다. 단, 바코드는 의미가 없으므로 --barcode_type 'not_barcoded' 옵션을 준다. 상당히 번거로왔지만 테스트는 잘 되었다.

마지막으로 가장 무식한 방법을 생각해 보았다. fastq 파일에서 barcode 판독 정보를 뽑아내는 것이다(multiple_extract_barcodes.py). 그 다음에 이들을 concatenation하여 도로 multiplexed fastq로 전환한 후 split_libraries_fastq.py를 실행하는 것. 원본 파일에 가까운 것을 재구성할 수는 있지만 너무 미련한 방법이다. 그리고 이미 demultiplexing이 된 파일을 입수해 보면 바코드 정보가 fastq 파일 내에 남아있지 않은 경우가 많다. 이론적으로는 안 될 이유가 없지만 시도하지는 말자.

2019년 9월 12일 업데이트

나의 위키 사이트(genoglobe.kr/kribb)에 일루미나 16S rRNA 시퀀싱 데이터를 처리하는 방법을 체계적으로 설명한 글을 올렸다.

2019년 9월 3일 화요일

간단한 AWK 활용 팁 - 컬럼 수의 총 합계 구하기 등

QIIME에서 사용하기에 좋은 매우 간단한 awk 활용 예제를 소개해 본다. OTU map file은 similarity에 의해 묶인 서열의 클러스터를 표현한 tab-delimited text file이다. 한 줄은 하나의 클러스터이며 맨 왼쪽은 클러스터의 ID, 나머지 컬럼은 클러스터를 이루는 각 서열의 ID이다.


총 몇 개의 서열이 클러스터를 이루는 것일까? awk가 이를 간단하게 계산해 준다.

$ awk 'BEGIN{s=0}{s+=NF-1}END{print s}' final_otu_map.txt

NF는 awk의 내장 함수로서 컬럼의 수를 의미한다. 첫번째 컬럼은 클러스터 ID므로 NF에서 하나를 뺀 것을 더해나간 뒤(변수 s에 저장), 파일의 모든 라인에 대한 처리가 끝나면 인쇄를 하면 된다.

다음은 조건에 맞는 줄은 편집을 하고, 그렇지 않은 줄은 그대로 출력하는 방법을 알아본다. 작년에 'awk를 이용한 같단한 텍스트 파일 조작(join)'이라는 제목의 글로 그 기본 원리에 대해서 이미 포스팅을 한 적이 있었다. 다음은 Gopalakrishnan의 2018년도 Science 논문에서 다룬 16S rRNA sequencing 결과 중 하나인 ERR2162226의 일부를 표시한 것이다.


엄밀히 말하자면 이 파일은 MiSeq가 생산한 그대로의 파일이 아니다. Wargo.R2d2_45029를 서열 ID로 전환하는 것이 오늘의 목표이다. '/1'은 제거해야 한다.

Awk를 사용하면 다음과 같이 명령행 한 줄로 끝난다.


$ awk '/^@ERR2162226/{printf "%s%s\n", "@", substr($2,1,length($2)-2);next}{print}' ERR2162226.fastq > ERR2162226-mod.fastq

awk '조건{명령A; next}{명령 B}' 텍스트파일

이 명령을 이해하는 것이 핵심이다. 조건에 맞으면 A를 실행한 뒤 next를 통해서 뒷쪽의 {명령B} 블록을 건너뛰고, 조건에 맞지 않으면 곧바로 {명령B}를 실행한다. if..else에 해당하는 문법이 아니다!. 만약 next가 없으면, 조건에 맞는 줄에서 명령A를 실행한 뒤, 그 줄에 대해서 또, 명령 B를 실행하게 된다. 이것을 가지고서 fasta로 전환한 뒤(split_libraries_fastq.py의 -q 19 옵션을 적용하여 quality에 의해 필터링하는 것을 똑같이 구현하려면 쉽지는 않음) pick_otus.py 명령을 실행하여 실습을 진행하면 된다.

awk의 printf 함수와 substr 함수의 용법도 이번에 공부를 해 두면 도움이 될 것이다.