'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 코드가 불러일으키는 난잡함, 그러나 그 속에서 풍겨나는 익살과 아름다움을 감상하고 싶다면 위키피디아 혹은 다음의 수상작 페이지를 방문해 보라.

• 제3회 수상작
• 제4회 수상작
• 제5회 수상작

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

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

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

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

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

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

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

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

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

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

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

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

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

화성으로 가는 탑승권

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

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


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

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

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

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

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

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

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

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

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

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

2019년 9월 13일 업데이트

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

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 사이트를 소개한다(링크).

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 (Edgar, 2010) 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도 넣어줘야 되겠다.

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

다음의 글을 참조하라.

[차이점] 개발자 vs 프로그래머 vs 컴퓨터 사이언티스트 vs 소프트웨어 엔지니어.

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

출처: Scott Hasnselman