2010년 10월 6일 수요일

[HMA] 01. Haeyoung's microarray analysis - 시작하는 글

먼저 [HMA](Haeyoung's Microarray Analysis)로 시작하는 글들은 웹을 통해 마이크로어레이 데이터 분석 방법을 소개하는 목적으로 쓴 글이 아님을 분명히 밝혀 둔다. 1998-1999년에 cDNA microarrayer와 scanner를 직접 만들어보고자 노력을 한 일이 있었으나, 실제 나만의 microarray data를 얻은 것은 작년이 처음이었다. 스탠포드대의 패트릭 브라운에 의해 이 개념이 도입되어 널리 퍼진지도 벌써 15년 가까운 세월이 흘렀다. 나는 그 사이에 오로지 미생물 genome sequencing/analysis에만 매달려 있었고, microarray 실험이 돌아가는 것을 옆에서 보기만 하였지 직접 내 손으로 해 본 일은 없었다.

이제 세월이 많이 흘러서 더 이상 컴퓨터를 이용한 기능 예측에만 머물 수는 없게 되었다. 그동안 애써 외면하고(?) 있었던 마이크로어레이 실험과 분석에 대한 경험을 직접 쌓아야 될 필요성을 절실히 느끼게 되어, 이제 나름대로 공부를 시작하면서 관련 정보들을 기록하기 위해 이 글을 시작하게 된 것이다. 나 자신도 이제 시작하는 입장일진대, 감히 다른 입문자를 위한 길잡이 노릇을 할 생각은 언감생심이로다.

자, 그러니 혹시 마이크로어레이 데이터 분석에 대한 기초를 공부하고자 웹 검색을 통해서 이 블로그를 들어오신 분은 실망하셔도 어쩔 수 없다. 그런 개인적인 목적으로 작성한 글이라면 뭐하러 웹에 공개하느냐고 물으신다면... 사무실이든 집이든 편하게 글을 작성하기 위해서 그렇게 결정했다고 어줍잖은 변명을 늘어놓으련다.

나에게 있는 것은 무엇인가? Perl에 대한 지식, 몇 편의 논문과 단행본, 그리고 몇 장의 GenePix 4000B 데이터들이다. 기본적으로 R을 이용한 데이터 분석을 시도해 볼 것이다. 필요하다면 TM4 software suite의 MultiExperiment Viewer나 CLC Genomics Workbench(상용)도 시도해 볼 것이다. 그러나 기본을 이해하는 데에는 R이 좋은 환경인것 같다.

2010년 9월 16일 목요일

Paired end sequencing vs Mate pair library sequencing

Sanger chemistry를 이용한 전통적인 샷건 시퀀싱에 익숙해 있는 내게 paired end sequencing이나 mate pair (library) sequencing이나 다를 바가 없다. 그런데 며칠 전 Illumina의 시퀀싱에 대해 검토하다가 이 두 가지 방법을 다르게 기술하고 있는 것으로 알고 깜짝 놀랐다. 라이브러리를 만드는 방법 자체가 다르기에 이를 구별할 수도 있겠지만, 나같이 전통 기술에 익숙한 사람에게는 혼동을 불러 일으킬 수도 있다. Illumina에서 이야기하는 mate pair library는 일종의 jumping library라고 하는 것이 기술적으로 더 정확할 수 있겠다.

Paired-End Sequencing - Acheving maximum coverage across the genome (Illumina)

Mate Pair Library Sequencing - Characterization genome variation (Illumina)


플라스미드에 클로닝하여 만든 라이브러리를 가지고 Sanger 방법으로 시퀀싱을 하는 경우를 생각해 보자. 클로닝 사이트 한 쪽의 벡터 서열에 붙은 프라이머로 시퀀싱을 하면 single end sequencing이고, 반대편의 프라이머로 한번 더 읽으면 이것이 바로 paired end sequencing이다. 이렇게 얻어진 한 쌍의 read를 mate라고 부른다. Paired end sequencing을 통해 얻어지는 read의 길이는 single end sequencing과 다를 바가 없다. 라이브러리를 만드는 방법도 역시 똑같다. 혹시 genomics의 석기 시대(?)에 M13 박테리오파지에 클로닝하여 ssDNA를 얻은 다음 시퀀싱을 하던 시절이 있었다면, paired end sequencing을 하려면 고생 깨나 했을 것이다. 그러면 colony picking이 아니고 plaque picking을 했었을까? M13에 클로닝한 뒤 PCR을 해서 시퀀싱을 했다는 전설같은 이야기를 들은 기억이 있는 것도 같다. 그러나 아무리 M13이라 해도 대장균을 좀먹는(?) 시스템을 가지고 대용량 시퀀싱을 하기는 좀 위험할 것이다.

Paired end sequencing이 왜 유용한가에 대해서는 이 글에서 설명하지 않아도 쉽게 알 수 있을 것이다. Resequencing이 일상화된 요즘, structural variation의 검출이라는 새로운 용도가 더해졌다는 것만 강조하자.

Genomic fragment를 뒤집어 붙여 원형을 만들어서, 수 kb 이상 멀리 떨어진 부분을 DNA 상의 연속적인 서열로 만들어 버린 뒤 이용하는 jumping library는 누가 가장 먼저 생각해 냈을까? 2005년 Science지를 장식했던 Jay Shendure의 polony sequencing에서 처음 쓰였던 것일까?

http://www.sciencemag.org/cgi/content/full/309/5741/1728

번거로와서 더 이상 거슬러 올라가지는 못하겠다. 여기서 mate-paired library란 말을 쓰고 있기는 하다. 그렇다면... mate에 해당하는(실제로는 fragment 길이만큼 떨어져 있어야 하지만) 한 쌍의 서열을 정말로 물리적으로 아주 가깝게, 짝을 이루도록(paired) 이어 놓았다는 것으로 이해해야 하나? ABI SOLiD의 기술을 소개하는 페이지에서도 mate-paired library라는 용어를 쓰고 있다. 자, 그럼 중간 결론을 내리자. 일반적인 의미의 mate pair와 mate pair(ed) library는 다르다! 그런데 위키피디아에서는 아예 Paired-end Tags라는 항목까지 있다. Paired end sequencing과 혼동할 소지가 있다...

http://en.wikipedia.org/wiki/Paired-end_Tags

정말이지 용어의 순화가 필요하겠다. 위 단락에서 빨간 색으로 표시한 것에 해당하는 라이브러리 및 이에 대한 시퀀싱을 가리키는 말을 다른 것으로 바꾸면 안될까?

내 기억이 맞다면, 원래 mate pair는 클로닝된 DNA 단편의 양 끝에서 만들어진 read의 쌍을 의미할 뿐 그 이상도 이하도 아니었다. Mate 혹은 mate pair라고만 쓰였지 Shendure에 의해서 특별한 방법을 사용해 만들어진 library로서 mate pair library라는 말이 쓰인 것 같다. 시퀀싱의 측면에서 본다면 mate를 만들어 내는 paired end sequencing과 쓰임새가 같지만.

생각난 김에 검색을 좀 해보자. 아하, "pairwise end sequencing"이 가장 역사가 오랜 용어로구나.

Pairwise end sequencing: a unified approach to genomic mapping end sequencing. Genomics 26:345-53 (1995)

이 논문에서는 어디에서도 "paired end"라는 말은 나오지 않는다. 하지만 pairwise end나 paired end나 뭐가 다른가...

2010년 9월 15일 수요일

Illumina-Solexa read 다루기

Sanger Institute에서 제안한 FASTQ format은 서열과 quality score를 하나의 파일에 축약해서 표연하는 방식으로서 high-throughout short read data를 다루기에 적합하며, 많은 de novo assembler 또는 reference mapper 프로그램이 입력 파일의 형태로서 요구하고 있다. 원래 quality score(Phred score)는 십진수 0~93으로 로 표현되지만, 이를 FASTQ format에서는 인쇄 가능한 아스키 문자로 대체하고 있다. 98, 99와 같은 높은 값의 quality score는 사용자가 편집한 값을 표현하는 용도로 쓰인다.

ASCII 32는 공백이므로, 인쇄 가능한 형태의 문자로 표현되는 ASCII 33~126이 FASTQ에서의 quality score로 쓰인다. 그러면 Phred score를 FASTQ format으로 표현하려면 어떻게 해야 하는가? 93보다 큰 값은 93으로, 이보다 작은 값은 그대로 둔 뒤 33을 더한 다음, 이 ASCII 값에 해당하는 문자로 변환하면 된다. 요약한다면 0~93까지의 Phred score를 33~126 범위의 숫자로 변환하여, 문자로 전환한다는 뜻이다.

이를 간단한 Perl code로 표현해 보자. Phred quality를 $Q, Sanger FASTQ의 값을 $q라 하면,

$q = chr(($Q<=93? $Q : 93) + 33); $Q = ord($q) - 33; Solexa/Illumina platform에서 생산되는 값은 조금 달라서, 0~62 범위의 값을 만들어 낸다. 이 경우에는 다음의 식을 적용하라. $Q = 10 * log(1 + 10 ** (ord($sq) - 64) / 10.0)) / log(10); http://maq.sourceforge.net/fastq.shtml에서 참조.

ord() 함수는 인수로 주어진 문자의 ASCII 값을 반환한다. Nucleic Acids Research 38(6) 1767-71에 Sanger FASTQ file format과 Solexa/Illumina FASTQ variant에 대한 상세한 내용이 나와있다.

maq를 이용해서 Illumina-Solexa 결과를 Sanger type의 FASTQ로 바꿀 수 있다.

maq sol2sanger s_1_sequence.txt s_1_sequence.fastq

s_1_sequence.txt가 Solexa read sequene file이다.

[2012년 7월 11일에 추가 작성]

quality score를 ASCII 값으로 전환하기 위해 더해주는 기본 값이 Sanger standard FASTQ(64)와 Illumina FASTAQ(33)가 서로 달랐었지만, CASAVA(Illumina pipeline) 1.6부터인가는 64로 통일되었다. 따라서 2011년 봄 이후로 생산되는 일루미나 데이터에 대해서는 quality encoding의 호환성에 대해 고민할 필요가 없어졌다. 매우 다행스런 일이다.

2010년 5월 31일 월요일

Hypothetical Proteins

아직 논의가 완전히 끝난 것은 아니지만, 4월 말에 참석했던 NCBI Prokaryotic Annotation Workshop(NCBI PAW)은 genome project의 결과물을 수시로 등록하는 나에게는 매우 중요한 행사였다.

해묵은 논점 중의 하나가 기능을 모르는 단백질의 이름을 어떻게 붙일 것인가 하는 점이다. 그러면 hypothetical protein은 도대체 무엇인가? 전통적으로 이 명칭은 다음의 세 가지 경우(엄밀히 말하면 모두 다른 의미)에 대해서 쓰여 왔다.

1. A protein whose existance is in question
2. A protein whose function/process is unknown or not predictable
3. Both of above

문법적으로 가장 적합한 것을 고른다면, 1이 되겠다. 그러나 NCBI PAW를 통해서 다음과 같은 합의를 이끌어 낼 수 있었다. 첫째, 단백질의 이름은 그 단백질이 존재함을 주장하는 기구는 아니다. 마찬가지로, 단백질의 이름이란 그 단백질의 실험적 특성 결정 수준을 주장하는 기구도 아니다.

Hypothetical protein, uncharacterized protein, protein of unknown function, protein of unassigned function의 네가지 중에서 protein of unassigned function이 가장 유력하다.

potential, precursor, conserved, unique, protein of unknown function, novel, fragment, similar to...

자, 그렇다면 gene prediction program으로 예측은 되었으나 다른 증거는 하나도 없는 단백질, 그리고 기능을 모르는 단백질에 대해 상동성이 있는 단백질을 서로 구별할 것인가? 그동안 널리 써 오던 명칭으로 부른다면 전자는 hypothetical protein이요, 후자는 conserved hypothetical protein 정도가 되겠다. AutoFACT 식으로 이야기한다면, 전자는 Unclassified, 후자는 Unassigned protein이다. 이를 모두 뭉뚱그려서 protein of unassigned function으로 해도 좋을까?

2010년 4월 12일 월요일

gsAssembler 2.0.0을 이용한 hybrid assembly

Sanger read(screened fasta file and quality file)내부에 각 read의 정보를 넣는 것이 관건이다. 우선 다음과 같은 스크립트를 만들어서 detPair.pl이라 이름을 짓는다. 여기서는 *.f1 또는 *.r1이라는 형식의 read tag을 쓰게 되어있다. 나는 주로 St. Louis convention(*.b1, *.g1)을 사용하므로 이를 약간 수정해야 했다.


그 다음에 할 일은 phred와 cross_match를 통해 screened sequence file을 만드는 것이다. edit_dir에서 phredPhrap 한번만 타이프하면 되지만, 공부를 하자는 의미에서 순서대로 적어 보도록 한다. Sanger ab1 file들이 chromat_dir에 있다고 가정하자.
phred -id chromat_dir -pd phd_dir
(determineReadTypes.perl은 별로 필요하지 않다)
phd2fasta -id phd_dir -os sequence.fasta -oq sequence.fasta.qual
cross_match -minmatch 12 -penalty -2 -minscore 20 -screen sequence.fasta
여기까지 하면 sequence.fasta.screen과 sequence.fasta.screen.qual 파일이 생길 것이다. 이제 다음과 같이 한다.
detPairs.pl sequence.fasta.screen > sanger.acc
fnafile -o sanger.fna -af sanger.acc sequence.fasta.screen
이제 만들어진 sanger.fasta와 sanger.qual을 gsAssembler에 넣어주면 된다.
다른 사용자를 통해 물어보니 Sanger paired end read를 넣어 주면 scaffold는 생성되지만, mate의 정보(clone insert size & direction)를 assembly에 직접 활용하지는 않는 것 같다. 만일 그렇다면, 라이브러리별로 크기에 대한 정보를 넣어 주어야 하나 그런 배려는 없다.
(참고자료: GS FLX software manual 5.5.4)

2010년 3월 29일 월요일

MIRA를 이용한 454+Sanger hybrid assembly

454에 기반한 fungal genome project를 하나 진행하고 있다. Genome size는 9 Mb 남짓이지만, 이것도 eukaryotic genome이라 세균에서는 접하지 못했던 약간의 고난이도의 repeat에 직면하게 되었다. GS assembler가 현재로서는 가장 만만한(?) 도구이지만, repeat에 의해 발생한 misassembly가 눈에 뜨이는 것은 사실이다. 정확도를 높이기 위한 방편으로서 mira를 활용해 보기로 한다.

mira 버전 3이 나오면서 설치 방법도 약간 까다로와졌다. boost C++ library(>=1.35)라는 것을 깔아야 하는데, CentOS 5에는 1.33이 설치되어 있으며 이보다 높은 버전의 rpm package는 아직 없다. 어찌어찌하여 boos 1.42의 소스를 구해서 설치하는데 이것도 쉽지가 않다. 압축만 풀면 단지 헤더 파일만 설치될 뿐이다. 나에게 필요한 것은 boost library binary이다. 설치 디렉토리에 가서 다음의 것을 해 주어야 한다.

$ ./bootstrap.sh --prefix=path/to/installation/prefix
$ ./bjam install

라이브러리 바이너리는 /usr/local/lib으로 하였다. $LD_LIBRARY_PATH를 설정해 주어야 mira의 빌드가 비로소 가능해졌다. mira는 매뉴얼 파일도 전부 tex 소스로 되어 있어서 이것 역시 애를 먹었다. make docs를 해야 man, html 등의 파일이 나오는데, latex2man이라는 유틸이 있어야 한다. 이것은 소스 rpm을 받아다가 설치했다. CenOS로 전환한 이후로 항상 yum으로 대충 패키지를 설치하고는 했는데, 소스 rpm을 풀어서 컴파일한 것은 이번이 처음이다. rpm -ivh <패키지명>을 실행하면 /usr/src/redhat/SOURCES에 tar.gz 파일이 생성된다. 이것을 풀어서 설명 대로 빌드하고 설치하면 된다.

Sanger/454 hybrid assembly

454 data만으로 작업을 하는 것은 비교적 쉽다. 하지만 Sanger data를 같이 섞어서 합체하는 경우 약간의 수고를 들여야 한다. Quality file을 수반한 FASTA sequence file, 그리고 NCBI TRACEINFO XML 파일이 필요하다. 최소한 vector clip은 되어 있어야 한다. 내가 쓰는 방식은 phred/cross_match를 통해 screening을 실시한 다음, PHD file을 참조하여 XML 파일 내에 clip quality left/right 정보를 넣는 것이다. 단, PHD file 내에 trim: -1 -1로 정의된 것은 clip quality position이 무의미하므로, 아예 서열 파일 뭉치에서 제거해 버리는 것이 낫다.

다음의 순서를 따르면 된다.

1. phred/cross_match를 실시하여 screened fasta file 및 qual file을 만든다.

2. trim range가 최소 100 bp인 read만을 골라내어, 즉 successful read만으로 구성된 screened fasta file 및 quality file을 만들어라. 어떻게 하냐고? 이럴 때 바로 BioPerl을 쓰는 거다.

3. traceinfo 파일을 만들어라.