Tabix: 엄청 빠르게 탭으로 구분된 1차원 좌표계 검색하기

시퀀싱 데이터를 대량으로 뽑아놓고 후속 통계 처리를 하자면,
각 태그가 각각 어느 영역에 들어가는지 분류하는 과정을 거쳐야 합니다.
요즘 cufflinks 같이 자동으로 요약해주는 녀석들이 많이 나와서 그냥 단순히 RNA-seq이나 ChIP-seq에서 보통 많이 하는 분석만 하는 것은 간단해졌지만, CLIP-seq처럼 어디로 분석이 튈지 모르는 다른 대부분의 경우엔 매핑된 곳이 뭔지 대량으로 정확히 분류하는 것이 중요합니다.

예를 들어 생쥐 유전체 mm9 버전에서는 Igf2가 chr7의 149,836,673부터 149,845,394까지 있습니다. 사이에 인트론이 3개있고요. 여기서 만약에 14984050 근처에 어떤 태그 하나가 매핑되었다고 치면, 이게 어느 유전자인가 물어봐서 “Igf2의 2번째 인트론입니다”하는 대답하는 것을 수천만번 매우 빠르게 할 수 있어야 하는데요. 여기서 좀 더 복잡해지는 것은 Igf2 인트론에 mmu-mir-483이 껴 있듯이, 여러 주석이 겹쳐서 시작과 끝이 순서가 달라질 수 있어서, 단순히 정렬된 리스트의 이진 탐색이나 B+ tree같은 것으로는 인덱싱이 불가능합니다.

이 문제를 해결하려면 가장 간단하게는 SQLdb에 모두 집어넣어놓고 SELECT절에 BETWEEN을 써서 검색하는 것이겠는데, UCSC Genome Browser도 이 방법으로 MySQL에 넣고 쓰고, RNA-seq 초창기 분석 도구로 유명했던 ERANGE도 sqlite에 넣고 계속 SQL 질의를 보내서 해결합니다. 간단하긴 한데 써 보면 엄청나게 느려서 분석을 구경하는 사람들에게 NGS란 이렇게 오래걸릴 정도로 대단한거구나 하는 환상을 심어주기 쉬울 정도가 되죠.

그렇지만 물론 BETWEEN을 수천만번 하는 것 보다 훨씬 빠른 인덱스를 만드는 방법이 많이 있겠죠. 대표적으로 R+treekd-tree 같은 2차원 공간 탐색이 가능한 인덱스를 쓰는 것도 있겠고요. 겹치는 녀석들끼리 링크드 리스트로 묶어서 서로 안 겹치는 클러스터로 만든 다음, 클러스터를 이진 탐색하는 방법도 있습니다. 알고리즘은 전자가 훨씬 멋있지만, 실제로 벤치마크해보면 후자가 압도적으로 빠릅니다.;; 이 방법을 구현한 파이썬 라이브러리로 Pygr의 NLMSA가 널리 쓰입니다. 그런데 이 쪽 구현은 내부적으로 깔끔하게 구현하느라 그랬는지, 메모리 소모나 인덱싱 속도, 검색 속도 모두 생각보다 훨씬 많이 먹고 느리고 그렇습니다. 특히 인덱스 들고 있는 데에 파이썬 객체를 너무 많이 만들어내서… 물론 SQL쓰는 것보단 훨씬 빠르지만 말이죠.

그러다 괜찮은 것을 발견했습니다. +_+ 어느날 samtools 새 버전이 나왔으면 받으려고 갔는데, tabix라고 처음 보는 게 올라가 있길래 궁금해서 보니까 바로 이 기능을 하는 것이더군요! 이건 좀 더 범용으로 만들어서 탭으로 구분된 모든 텍스트파일을 인덱싱 가능하도록 만들었고, 게다가 원본 그대로 gzip한 다음에 원본에서 찾아주는 것이라 여러 도구에 쉽게 녹아들어갈 수 있게 만들었네요. (압축파일 안에서 랜덤액세스가 가능하도록 하는 도구도 물론 들어있습니다.)

써보니 속도도 상당히 빠릅니다. 매우 만족! 그래서 파이썬에서 쓰려면 어떻게 해야하나 보니까, 이미 펄, 파이썬, 자바 바인딩이 들어있더군용. 그러나 파이썬 바인딩이 ctypes를 쓰게 돼 있어서, 설치나 관리가 여러모로 번거롭고 속도도 (아주 약간) 느려지고.. 흐흐 예전부터 들어왔던 Cython을 배워볼 절호의 기회다! 하고 Cython으로 한참 다시 바인딩하는 작업을 해 봤는데, 아무리 찾아봐도 이터레이터 객체를 구현할 때 __next__에서 NULL을 리턴할 방법이 없어서 깔끔하게 만들어 줄 수가 없네요.. 결국은 그냥 포기하고 C 모듈로 만들었습니다. (tabix 0.2.3용 패치) 패치는 tabix 원저자에게 보냈습니다. 파이썬 2/3 겸용이죠. 호호호;

이렇게 쓸 수 있습니다~ 일단 인덱스 만들고 확인하기.

파이썬에서도 똑같이 접근하기.

원래 pygr.NLMSA쓰고 있던 곳을 요걸로 바꾸니까 거의 10배 이상 빨라졌네요. ^__^

RNA 연구에서의 시퀀싱 활용 (시작)

대용량 시퀀싱(high-throughput sequencing; 다음부터 시퀀싱)이 폭발적으로 널리 쓰이기 시작하고도 이미 한참 지나서, 네이처 주요 논문들의 그림 구성, Affymetrix 주가, 학회장의 질문들 어느 것 하나 예전과 같은 게 없어졌습니다. 그래서 이제 논문 읽을 때 유전체학이나 전사체학 용어와 개념에 익숙하지 않고서는 읽기 힘든 논문이 한 둘이 아니게 됐는데…

시퀀싱의 주요 수요처인 유전체 DNA 시퀀싱은 워낙 뻔하고 여기 저기서 늘 얘기하고 있는 것이니 넘어가고, 역동성이 넘치는(-O-) RNA 분야에서의 시퀀싱 활용에 대해서 대략 정리해 볼까 합니다. 한꺼번에 다 올리면 읽기 압박이 있으니 여러 회로 끊어서 찔끔찔끔 올리겠습니다~ (히히)

순서는 모두 9번으로 나눠서,Roblox HackBigo Live Beans HackYUGIOH DUEL LINKS HACKPokemon Duel HackRoblox HackPixel Gun 3d HackGrowtopia HackClash Royale Hackmy cafe recipes stories hackMobile Legends HackMobile Strike Hack

  1. RNA 시퀀싱 기본 지식
    1. Illumina mRNA-seq 프로토콜 대충 보기
    2. Illumina SRA 프로토콜 대충 보기
    3. 여러가지 RNA/DNA ligase
  2. RNA 프로파일링
    1. RT하는 방법에 따라 다른 것
    2. 조각내는 방법에 따라 다른 것
    3. CAGE, SAGE에서 파생된 프로파일링 기법들
    4. 작은 RNA 프로파일링
  3. 번역효율 측정하기: ribosome profiling
  4. RNA와 단백질의 상호작용 보기
    1. RIP-seq
    2. CLIP-seq (UV cross-linking and immunoprecipitation), PAR-CLIP
    3. miRNA 타겟 보기: Ago HITS-CLIP
  5. RNA 간의 상호작용 보기: CLASH
  6. RNA 잘리는 부분 알아내기: degradome sequencing
  7. RNA 2차구조 보기: Kertesz et al.과 FragSeq
  8. RNA 5′ 3′ 끝 알아내기: High-throughput 5’/3′ RACE와 PET
  9. mRNA 3′ UTR 끝 알아내기: 3P-seq

로 올리도록 하겠습니다. ^__^

(다음 시간에 계속~)

후배를 찾습니다~

몇 군데에 이미 전에 올린 적이 있어서 이미 보신 분들도 있겠지만, 좀 더 많은 홍보를 (;ㅁ;) 위해서 잠잠한 블로그에도 올려봅니다. ^.^;;

제가 공부하고 있는 실험실에 새로 석사과정, 통합과정 또는 박사과정 대학원생으로 참여할 대학원생을 모집합니다.

저희 실험실은 동물 RNA의 유전자발현 조절 분야를 주로 연구합니다. 특히 microRNA의 생합성 경로와 조절 경로의 주요 단백질들의 기작을 밝힌 것으로 많이 알려져 있습니다. (지도교수님 최근 기사와 인터뷰, 조금 오래된 인터뷰 참조)

이번에 찾고 있는 대학원생(1명)은 실험실에서 하는 분야 중 유전체학과 생물정보학 쪽을 주로 하도록 뽑을 예정이고요. 따라서, 학부나 석사 전공으로 컴퓨터과학, 통계학, 생물정보학 등을 전공하신 분을 기대하고 있습니다. 원하는 경우에는 실험도 배워서 스스로 데이터 만들어서 분석하는 요즘 유행하는 스타일로 공부할 수도 있습니다.

특히 유전자 발현 조절 분야는 대규모 전사체/유전체 실험으로 계속 커지는 추세라서 대용량 데이터를 정량적으로 다뤄서 생물 연구를 할 수 있는 기회가 계속 늘고 있습니다.
실험실에서 최근에 많이 하는 high-throughput 실험들을 많이 도입해서 하고 있고, 원하는
실험을 필요하다면 웬만큼은 부담이 되더라도 할 수 있을 정도의 여유는 되기 때문에
연구할 데이터와 기회가 풍부하고요. 실험실의 주요 연구분야가 아직 미지의 영역이 매우 많은
분야라서 머리를 쥐어짜서 연구분야를 찾아 헤멜 필요도 없어서 연구 측면에서는 국내에서
대학원 생활하기에는 아주 좋은 여건이 될 것입니다~

입학전공은 생명과학부, 생물정보학 협동과정, 유전공학 협동과정, 종양생물학 협동과정 중 하나를 선택할 수 있습니다. 2011년 후기 이후 입학 가능하며 입시 이전에 협의가 되어야하고, 연구원으로 미리 일할 수도 있으므로, 되도록이면 빨리 연락을 주시는 편이 좋습니다. (보통 1년 이전에도 많이 연락하는 편입니다.)

실험실에 관한 정보는 홈페이지를 참조하시고, 문의사항이나 지원에 관련된 것은 모두 저한테 자유롭게 보내주세요. 메일 주소는 제 자기소개 페이지 맨 끝에 있습니다. ^.^; 혹시 주변에 관심 있을 것 같은 학생이 있으면 알려주세요~ (나쁜 아저씨 아니에요~;;)

생물정보 태동기의 재미있는 사실들

어느 학문 분야든 성숙하다보면 해당 분야의 역사와 철학에 대한 연구가 따라오게 된다.
학문이 생기게 된 배경과 발전 과정, 패러다임의 변화, 다른 학문에 대한 영향, 연구자들의 분야
고유적인 연구 방법을 관찰하는 것은 재미있지 않을 수가 없다.

최근에 PLoS Computational Biology에 생물정보학의 뿌리라는
기사가 올라왔다.
유전체 모델이나 RNA 2차구조 같은 것을 촘스키식 문법으로 다룬 것으로 유명한
David Searls가 쓴 생물정보학의 역사에 대한 글인데, 깊게 잘 다루었다.

철학적인 생각은 글에 남겨두고, 의외로 모르고 지나가기엔 너무 아쉬웠을 만한 재미있는 사실 몇 가지만 추려보면,

  • 컴퓨터를 생물 연구에 처음으로 쓴 사람은 너무 뻔해서 약간은 재미없게도(?) Ronald Fisher인데, EDSAC을 개발한 Wilkes와 Wheeler가 직접 작업을 돌려주었다. (1950년)
  • 소개가 필요없는 Alan Turing은 말년에 주로 발생학 연구를 했으며 (1952년~), 역시 정보이론과 논리회로의 창시자격인 Claude Shannon은 심지어 박사학위를 계산유전학에 대한 연구로 받았다. (1940년)
  • 빅뱅이론으로 유명한 이론물리학자 George Gamow와 Monte Carlo 시뮬레이션으로 유명한 이론물리학자 Nicholas Metropolis는 유전코드의 상세한 기전이 밝혀지기 전에, 서열의 통계적 분석과 시뮬레이션으로 유전코드의 이론적 특성 연구를 했는데 이 연구가 거의 역사 최초의 생물정보학 연구로 보통 받아들여진다. (1954년)
  • 역시 초기에 컴퓨터를 가장 널리 사용한 것은 결정학자들이었는데, 1952년에 이미 EDSAC으로 계산한 논문이 나왔다.
  • 또 다른 생물정보학의 주세부분야 중 하나인 계통분류계산은 1957년에 처음 시작되었다. 요즘 화학유전체학에서 거의 표준처럼 쓰이는 타니모토 계수는 1960년에 IBM의 수학자인 타니모토가 세균 분류를 위해 개발했다.

보통 어디서 트렌드따라 뚝 떨어진 신생융합듣보잡 취급을 많이 받는 생물정보학이지만 의외로 뿌리는 깊다. +_+

시퀀싱 데이터에서 3′ 어댑터 서열 제거

small RNA 시퀀싱에서는 리드보다 RNA가 더 짧아서, 5′ 끝부터 읽을 경우에는 3′ 어댑터 시퀀스가 나오고, 프라이머를 뒤집어서 3′ 끝부터 읽으면 5′ 시퀀스가 나올 수 밖에 없다. small RNA가 아니더라도 CLIP에서는 보통 바인딩 사이트를 정확히 알기위해 짧게 쳐내서 시퀀싱하는 경향이 있어서, 보통 30nt 안쪽으로 들어오는 편이라 시퀀싱한 뒤에 어댑터 제거가 꼭 필요하다.

그렇지만 역시나 PCR 오류, 시퀀싱 오류, 어댑터 불량 등등 수많은 잡음때문에 역시 단순 문자열 비교로는 잘 안 통한다. 그래서 정규식을 쓰기도 하는데 영 속도가 만족스럽지 못하고, 모든 자리에서 어댑터 시퀀스랑 비교해서 미스매치를 세는 등의 방법(HTSeq 패키지)을 쓰기도 하는데, 갭을 전혀 허용하지 않아서 어댑터 합성 품질이 안 좋은 경우는 놓치는 것이 너무 많아서 결과를 보면 답답~하다.

최근 많이 쓰이는 방법으로 Needleman-WunschSmith-Waterman을 섞어서 어댑터의 5′ 끝에게는 지역정렬처럼 아무데서나 시작하게 하고, 3′ 끝에는 전체정렬처럼 끝까지 가게 하는 것이 있다. 그런데, 소프트웨어는 홈페이지에 오거나, 저자에게 말하면 준다고 다들 써 놓고서는 정작 홈페이지에 가면 아무 것도 없고, 메일 보내면 묵묵부답이라, 1년 넘게 정규식으로 불쌍하게 쓰다가, 결국 큰 마음먹고 토요일 밤을 투자했다. +_+

소스코드 받기 (파이썬용 C 확장 모듈)

실제로는 affine gap penalty를 써서 행렬이 3개지만, 그냥 linear gap penalty를 쓴 경우라고 하고 행렬을 예를 들어 보면,

  • 시퀀스 리드: CCAGTCCA
  • 어댑터 시퀀스: CCAG
  • 점수: match 2, mismatch -3, gap penalty -3

3' 어댑터 떼기Watch movie online The Transporter Refueled (2015)

이렇게 하면 짠~

구글 리더에서 Papers로 URL열기

대학원생의 친구, Papers!
공부하는 척 할 때 참 좋은 녀석이지만, RSS기능이 없어서 새 글을 받아보려면 반드시 나머지
반쪽이 필요한 어정쩡한 녀석이죠~

NetNewsWire에 연결해서
보면 구글리더와 동기화해서 돌아다니면서 버스에서도 볼 수 있고 참 좋긴 한데.. 문제는
구글리더가 최근 글 10개만 주는 바람에, AOP feed가 없거나 글이 한꺼번에 잔뜩
올라오는 PNASPLoS ONE같은 경우에는 거의 글을 대부분 놓쳐버려서
가끔 구글리더 들어가서 확인해야해서 여러모로 귀찮습니다.

결국 그냥 구글리더로 싹 읽으면 깔끔하고 좋기는 하지만, Papers에 논문 가져다 넣으려고 긁어 넣고 붙이고
하기가 귀찮아서 결국은 NNW에 남아있다가, 오늘 마음잡고 오랜만에 GreaseMonkey로 정리해 봤습니다~

Google Reader에서 c 눌러서 지금 글 Papers열기 (그리스몽키)

설치하시면 C 누르면 Papers에서 글이 열립니다~

BSD에서 타임머신처럼 백업 관리하기

엊그제 큰일날 뻔했습니다. 연구실 홈페이지 고치다가 www디렉토리를 복사한 다음에 지운다는 것이, 제 홈디렉토리에서 지워버리는 바람에 제 홈페이지가 몽땅 날아가 버린 것. –;
다행히도 작년 5월 백업도 있었고, 그 이후로는 업데이트를 거의 하나도 안 해서 간신히 살릴 수 있었지만 초등학교 때 만들었던 프로그램을 고등학교 다닐 때 모두 날린 이후로 최대 사건이 될 뻔 했습니다.;

데스크탑은 타임머신으로 잘 백업하고 있었는데 작년에 서버 옮기고 설정하기 귀찮다고 작업서버 홈 디렉토리 백업을 설정해 두지 않았는데요. 그래서 이참에 rsync로 타임머신하고 거의 똑같이 쓸 수 있다는 얘길 들은 생각이 나서, 설정했습니다. 잘 되더군용! ^_^ 안심

물론 파일은 대부분 하드링크 덕에 용량을 크게 차지하지는 않지만, 디렉토리 구조가 계속 복사되는 통에 생각보다는 용량을 꽤 차지했습니다. 타임머신은 최근 하루는 1시간 간격, 1주일은 하루 간격, 1달은 1주일 간격 이런 식으로 오래된 시간일 수록 띄엄띄엄 저장할 수 있도록 돼 있는데요. 그래서 이걸 어떻게 해야하나 검색을 좀 해 보다가 적당한 것이 안 보여서 간단하게 만들어 봤습니다.

소스파일 다운로드 (파이썬 2용)

음 마음이 편안하군요.~

근황

거의 한 해 동안 글을 안 썼습니다. 바쁜 일도 많았지만 안 쓰다 보니 안 들어오고, 안 들어오니 안 쓰고 순환의 연속으로… 크크.
그래도 여지껏 RSS 구독을 남겨두신 분들께 혹시 궁금하시면 요즘 어떻게 살고 있는지
알려드리려고 근황을 남겨둡니다.

셀카질
올림푸스 E-P1 산 기념으로 시험 셀카;

작년 2월 말에 대전에서 졸업하고, 서울로 이사했습니다. 요즘은 낙성대역 근처에 삽니다. 보기보다 꽤 살 만한 동네입니다. 언덕이 많아 눈 쌓이면 매우 곤란한 점만 빼면. 처음 몇 달 간은 아침에 바삐 출근하는 사람들을 보기가 좀 갑갑했죠. 대전에선 하늘이 넓은 곳에서 여유롭게 돌아다녔는데 아무래도 서울은 좀 달라요.

직접(?) 구운 RNA 쿠키~
직접(?) 구운 RNA 초코쿠키

작년 9월부터는 일하고 있던 연구실에서 박사과정을 시작했습니다. 2009년 2학기 시작이니까, 보통 4~5년 정도 한다고 생각하면 2013년이나 2014년 정도까지는 계속 학생입니다. (히히) 새로 옮긴 연구실은 microRNA라는 생분자를 연구하는 곳입니다. 실험실 분위기가 아주 좋아서 매일 출근하는 게 즐겁습니다. 지도교수님은 잘 모르는 분야의 얘기라도 호기심을 가지고 항상 관심있게 들어주시고, 학문적으로도 놀랍도록 식견이 있으시지만 인품도 모두가 내가 연구책임자가 되어도 저렇게 할 수 있을까 싶을 정도로 좋으셔서, 일할 수 있는 가장 좋은 환경에서 행복하게 일하고 있습니다. ^_^

실험실 설정샷
일하는 척 설정샷. ㅋㅋ;

연구실에서 생물정보를 전공한 사람은 혼자라서, 여러 프로젝트에서 나오는 대용량 자료 처리는 대부분 맡아서 하고 있습니다. 그 덕에, 많은 사람들과 많은 연구주제에 동시에 참여할 수 있어서, 다양한 분야를 배울 기회가 됐습니다. 혼자하는 주력 연구로는 새로운 작은 RNA 발견을 위한 유전체학적 분석을 하고 있습니다. 돈도 많이 들고 유독물질과 방사능도 많이 접하게 돼서 원래 하던 일처럼 편하지는 않지만, 이제 실험도 어느 정도는 적응이 돼서 재미있게 하고 있습니다. 이히

내 자리내 실험대
공부하는 자리와 실험대. 사진 찍을 때는 HHKB였지만 지금은 Filco 쓰고 있어요. 파이펫은 길슨.

그리고 토요일마다 실험실원 여러 명에게 프로그래밍을 가르쳐드리고 있습니다. 요즘 생물 데이터가 워낙 대용량화되는 추세라서 뭘 하려면 정보의 흐름을 다루는 능력이 필요해서 다들 흥미를 가지고 참여하고 있습니다. 파이썬으로 하고 있는데, 전에 C++을 잠시 배운 적이 있었던 사람들이 “프로그래밍이 이렇게 재미있는 것인 줄 몰랐어요!”라고 합니다. 역시 파이썬! ㅋㅋ;

연구실 밖 풍경
연구실 자리에서 보이는 바깥 풍경. 쭈욱 오르막이라 모두 가까이 보여서 좋음 +_+

주 5일제 하던 곳에서 주 6일제인 곳으로 옮기다보니, 사실 대전에서 서울에 왔어도 오히려 사람을 만날 수 있는 시간이 더 줄었습니다. 그래서 거의 모임이 있어도 못 가고, 점점 사회에서 떨어져서 산에 사는 사람이 되는 느낌이.. >_<.

저도 작년 예약판매할 때 아이폰을 샀습니다. 원래 아이팟 터치를 늘 친구처럼 데리고 다녀서 소지품 수를 줄이는 효과도 좋지만, 거의 생활을 완전히 바꿔놓은 놀라운 기계네요. 화장실에서 책을 안 읽게 되었다는 슬픈 단점도 있습니다만.. 크… 몇년간 IT 관련해서는 뭘 배운 것이 없었는데, 아이폰 프로그래밍도 조금씩 해 보고 있습니다. 꼭 뭐 만들어 봐야지! 히히

스탠포드 로댕 미술관
스탠포드 로댕 미술관 앞에서 실험실 동생과 동상 따라하기~

올해 초엔 처음으로 미국에 갔습니다. 아 말로만 듣다가 직접 가 보니 음식도 맞고 사람들도 괜찮고 좋네요. 다음에 또 가려면 열심히 연구해서 학회에 뭔가를 꼭 내야겠어요. (동기유발 효과가 불끈불끈) 콜로라도에 있는 키스톤 리조트에서 관련분야 학회에 참가한 뒤 샌 프란시스코에 갔습니다. 샌 프란시스코는 항상 날씨가 좋고 좀처럼 비가 오지 않는다는데, 제가 갔던 기간 중엔 폭풍우가 몰아치더군요.; 스탠포드 로댕 미술관. 그동안 미술관에 여러 번 갔지만 감동을 느끼기는 처음이었습니다. 작가가 그리면서 느꼈던 감정을 그대로 느끼는 듯한 감동! 이히히 다른 미술 전시회도 다시 가야 겠어요.

이제 봄 입니다. 돌아올게요. 종종 또 보아요!

군집이룬 자료의 비모수 두 표본 차이 검정 (C 구현)

두 표본 집단의 차이가 의미가 있나 알아볼 때 보통
t-검정을 많이 씁니다.
t-검정에서는 값의 분포가 정규분포라는 가정이 있어서, 자료의
분포를 모르거나 정규분포로 바꾸기 매우 난감한 경우에는 t-검정을
할 수 없어서 비모수 검정을 사용하는데요, 이쪽으로 가장 인기있는
검정법은 아무래도 Mann-Whitney U test입니다.

어느날 MW-U로 재미나게~♪ 통계를 하다가~ 아니!! 결과가 엄청 편향되어 나오는 것입니다!
비모수라더니!! 마이크로RNA
마이크로어레이로 나온
결과를 분석하고 있었는데, 마이크로RNA가 염색체에서 떼로 몰려 있어서
한 놈이 올라오면 같이 우루루 올라오는 특성이 있다보니 군집이 엄청 큰
놈들에 의해 전체가 흔들려서, 아 무슨 좋은 방법이 없을까! 하다가
군집을 이룬 자료들에 쓸 수 있게 변형한 것
(Rosner and Grove, 1999,
Datta and Satten, 2005,
Haataja et al., 2009)
들을 발견했습니다. +_+

오…. 다 괜찮아 보이지만, 결국은 코드가 공개돼 있는
Datta와 Satten의 것으 로 해서 일단 결과를 뽑았습니다. 매우 만족스럽군요! ^_^

그런데 문제는 순수 R로 구현되어 있다보니
속도가 엄청나게 느려서, 자료가 별로 크지도 않은데 거의 2시간씩 돌려야
돼서 시험삼아 돌려볼 때도 마음을 크게 먹고 돌려야 해서, 바로바로
결과를 확인하는 재미가 없었습니다.

그래서 쭉 벼르고 있다가, 주말에 여자친구가 기말고사 공부해야 해서, 같이 공부하는
척 하느라 짬이 난 김에
C로 새로 코딩했습니다 (다운로드).
R이나 SciPy같은데서 쉽게 쓸 수 있게 외부 의존성을 없애려고, Z-통계치에서
p-value구하는 부분은 빼서 의존성을 줄였습니다. 그냥 libm만 있으면 됩니다.
(파이썬에서는 scipy.stats.norm.pdf를 써서 Z에서 p-value를 구할 수 있습니다.)
대략 돌려봤더니 2시간 걸리던게 30초로 줄었군요! 이히히 ^__^*
간단하게 컴파일한 다음에 파이썬에서 쓰려면 이렇게 쓸 수 있습니다~

으음.. R에 붙이는 것은… 아직 잘 몰라서… (나중에..;)