NGS를 위한 생물정보 프레임워크의 필요성

최근 차세대염기서열결정법(NGS)이 각 생물학 연구에서 주요한 도구로 활용됨에 따라 생물정보 분석 서비스를 위한 하드웨어 및 소프트웨어 일체를 NGS에 맞출 필요가 있습니다. NGS는 기존 염기서열결정법에 비해 짧은 서열을 대규모로 생산하므로 이에 맞춰 메모리 활용이나 병렬 계산법이 적극적으로 적용된 소프트웨어가 개발되고 있고, 따라서, 최적의 성능을 발휘하기 위한 하드웨어의 구성도 함께 필요한 것입니다.

최근 발간된 보고서(The Global Outlook for Next Generation Sequencing: Usage, Platform Drivers & Workflow, October 31, 2011. BioInformatics, LLC. See on YouTube)에 따르면 NGS를 이용하는 연구자들의 이슈 중에서 가장 어려운 문제로, 분석 소프트웨어의 성능 개선으로 조사되었음. 또한 플랫폼 관리 및 스토리지 문제를 포함하면 생물정보 관련 문제가 29%나 차지하고 있습니다.

그림 1. Most Significant Improvement to Your Next Generation Sequencing Workflow (출처: The Global Outlook for Next Generation Sequencing: Usage, Platform Drivers & Workflow, October 31, 2011. BioInformatics, LLC)

또한 이 보고서에 따르면, 현재의 연구 환경에서 연구자들이 느끼는 가장 큰 병목 지점으로 워크플로우(work flow)의 투명한 관리에 있다고 합니다. 즉, 워크플로우를 이용한 여러 단계의 심층 분석의 결과와 각 단계별 결과의 이력을 투명하게 관리하는 정보공학적 플랫폼이 가장 큰 병목 문제로 지목하였습니다. 이와 함께 워크플로우의 아웃소싱을 전략적으로 수행하기 위한 체계 구축도 중요한 이슈로 부각되었습니다. 오믹스 프로젝트가 생물학자, 시스템관리자, 개발자, 생물정보학자 등 다양한 유형의 인력이 협업을 해야하고 문제를 해결하기 위한 워크플로우 자체가 복잡하게 얽혀 있어서 어떤 부분을 나누어서 아웃소싱을 해야할 지 업무를 분장하기 어려운 것입니다.

이러한 경향은 NGS 고유의 문제라기보다는 소프트웨어의 일반적인 경향으로 파악할 수 있으며(그림 2, from www.geospiza.com), 우수한 품질의 소프트웨어 서비스를 개발하기 위해서는 하드웨어, 소프트웨어, 휴먼웨어의 상생적 역할이 중요하다고 생각합니다. 게다가 NGS는 새로운 분석 프로토콜이 다양하게 공개되고 있으므로 신속하게 수용하고 서비스할 수 있는 체계가 필요합니다. 

다음 글에서는 오픈소스와 상용 소프트웨어어 프레임워크에 대해서 다뤄 보겠습니다.

그림 2. 소프트웨어 빙산, 편리한 화면요소와 신뢰성 높은 알고리즘을 위해서는 내제된 많은 요소기술이 뒷받침되어야 함 (www.geospiza.com)

Collecting meta data from Entrez

It's often to show growth of sequence data of interest when one writes research proposal. For an example, you requires to collect number sequences from agricultural organisms and compare it to human if you want to explain how sequences regarding to agricultures grow faster than human data. Usually the gross statistics of GenBank, is posted on NCBI's Web page, might be not enough to describe details of the data growth.

By using show index, preview, and limit functions in Entrez, you can quickly collect meta information like number of entries.

dbEST	Total records	Records for last 3 years	Growth rate for last 3 years
human	8,315,231	177,492	2.1%
mouse	4,853,547	3,289	0.1%
cattle	1,559,494	45,232	2.9%
pig	1,620,570	144,207	8.9%
chicken	600,423	1,041	0.2%
insects	4,493,137	1,864,326	41.5%
bacteria	1,266	1,012	79.9%
fungi	2,893,583	1,508,814	52.1%
plant	22,633,681	7,290,397	32.2%

To complete the above table, we need to count total records for each species in dbEST at current date. It's quite simple as the following: 1) choose a database, EST, on NCBI's front page; 2) press search button; 3) click advanced search; 4) choose an "Organism" field; 5)  type "human" on query box; 6) click show index; 7) press "add to Search Box"; 8) press preview button (not search here),  and 9) repeat for every species from 3).

Select DB of interest

In advanced search

Now you can validate your query word whether correct term in Entrez and then choose "human (10064832)" among drop items. 

Here the number in parenthesis means number of entries for human organism. Please try yourself to choose "All fields" and same query, then compare the numbers. Maybe they are different counts, why?

One query was built and the results was stored as #4 symbol in Search History.

With repeats for every species, you can get similar results like the above.

Now we need to find whole numbers of EST for last 3 years. At first back to front page of dbEST and the click the limits link. Now you can see the following  screen and choose proper items in "Published in the last", and then press search button. When you see a list page of results, click "advanced search" again. If you see similar page like the below, it's correct.

In this page, #11 query means total numbers of ESTs published for the last 3 years

To find human ESTs published for the last 3 years, you better to use an "AND" Boolean operator with shortcut symbol. That means to type query as the following then push the "Preview" button.

By repeating for every organisms, you might get similar results like the below. The growth rates for each organisms should be calculated by your self.

Out-focusing of Lens

Photo 1. A toy car, a trial of out-focusing

Yesterday I went to a photo studio, LemonTerrace Dongtan-branch, to take pictures of my son for celebrating the first birthday of him, "Dol", a Korean word. The host of the studio, a specialist of camera, taught how to use my DSLR camera (EOS1000D) in details. Especially he demonstrated the out-focusing technique in case-by-case. A simple protocol for the beginner of DSLR user is the following; (1) Set 'Av' mode of your camera, (2) Set max number of ISO, my case is 1400, (3) adjust the distance between the camera and the subject as close as possible, (4) check shutter's speed by half pressing shutter button, (5) the speed should be lower than 1/60, (6) if too lower eg. 1/125, change ISO to be lower, (7) if higher speed, adjust the distance or set higher ISO, and (8) take a picture by full pressing. And he highly recommended for me to change my camera lens for portraits, 30 mm f/1.4. The studio facilities and peoples gave me very good impression.

Photo 2. Eunchae, my daughter, is cooking

Photo 3. Eunwoo, my son.

Protocol : Out-focusing using Cannon DSLR

1. Set 'Av' mode of your camera
2. Set max number of ISO, my case is 1400
3. Adjust the distance between the camera and the subject as close as possible
4. Check shutter's speed by half pressing shutter button, thenthe speed should be lower than 1/60
5. if the speed is too lower then 1/60, eg. 1/125, change ISO to be lower, eg. 400.
6. if higher speed, adjust the distance or set higher ISO
7. take a picture by full pressing.

END.

How to find a long indel from Nsp2 alignment

Motivation
Zhou et al. presented that an unique 30-amino-acid deletion in Nsp2-coding region is a key feature to classify whether a strain is a highly pathogenic porcine reproductive and respiratory syndrome virus (PRRSV). And the Nsp2, nonstructural protein 2, has been shown to undergo remarkable genetic variation, primarily in its middle region, while exhibiting high conservation in the N-terminal putative protease domain and the C-terminal predicted transmembrane region (Han et al. 2007). This post aims to show how to find a quite large deletion in a specific coding-region with positional tolerance.

Figure 1. The 30-Amino-Acid Deletion in the Nsp2 of Highly Pathogenic PRRSV. 
(Zhou et al. 2009)

Method and Implementation
Pairwise alignment between a sequence of interest and a reference sequence (ORF1a of VR-2332 strain) is an essential step for finding insertions and/or deletions, shortly indels. The two sequences were aligned with BLAST (Altschul et al. 1999), exactly saying bl2seq with blastx program because my query are assumed as DNA sequences from clinical samples. Remind that global alignment like ClustalW (Thompson, J.D. et al. 1994) is not suitable for this case. The reason is that query sequences were quite shorter than a reference sequence. After blasting two sequences, a Perl script using BioPerl parses a report file from BLAST and evaluates whether indels' position and size is tolerable for user's constraints.


A format of BLAST (bl2seq) result consists of several sections including a summary of query, alignment reports for each HSP (high scoring pair), parameters, and statistics. An example part of the result is shown as the below; Each subject begins with a symbol, ">", accession number and description are followed. In next part, summary statistics for each HSP are written, including bit score, E-value, identities, etc. Actual alignment between the query and subject are described in the last part. These format of data are repeated for all HSPs and subjects.

Figure 2. An example BLAST's result of the deletion in Nsp2.

Bio::SearchIO module supports to read and write search results in various file format of biological applications such as BLAST, FASTA, HMMER, Sim4, etc. Basic parsing a file from BLAST was done by Bio::SearchIO functions and I recommend to read a HOWTO document on the BioPerl. A next step is to extract a quite large indels from HSP alignment information.


At the first time I thought too pedantic to implement this problem --- this is why I write this post. An example of too much 'pedantic' idea is like this; calculating % gap for each columns, cut-off with a specific threshold, finding a series of deletions without break using graph, and evaluating each deletions. It's too much complex! I stop my work. After drinking a cup of coffee, I considered it again in mind of "Perl script should be Perlish".


The following code is to get a string containing a alignment of query sequence potentially including indels marked as "-".


my $query_alignment = $hsp->query_string();


Now I can use super strong features regular expression of Perl for find a series of indels; this codes extracts positions of start and end residues for all possible long indels, which are referred to a reference sequence. And also indels is constrained with minimal length of deletion.



my ($ref_start, $ref_end) = $hsp->hit->strand < 0 ?
                                    ( $hsp->hit->end,
                                      $hsp->hit->start ) :
                                    ( $hsp->hit->start,
                                      $hsp->hit->end );

while( $query_alignment =~ /(\-{$threshold_length,})/g ) {
      $indel_start = $-[1] + $ref_start;
      $indel_end = $-[1] + $ref_start + length($1) - 1;
      $indel_length = length($1);
      push  @indels, [$indel_start, $indel_end];
}


Further filtering constraints were implemented and passed from arguments, threshold % indentification (-p), expected start position of indel (-s), expected end position of indel (-e), and tolerant positional error(-d). Sample result is shown as the following:


[bckang@gxs ~]$ ./gvs_find_indel.pl nsp2.blsx -s 921 -e 949 -d 15
#Mutation dectector's parameters: p=30 s=921 e=949 d=15
>AAO13191(NA_prototype)
920     948     29
>AAO13191(NA_prototype)
>AAO13191(NA_prototype)
>AAO13191(NA_prototype)
>AAO13191(NA_prototype)
>AAO13191(NA_prototype)

It is true the simplest way is the best. Also it is better way that Perl is supposed to be Perlish and Python looks like Pythonic.

References

1. Zhou, L. et al. (2009) The 30-amino-acid deletion in the Nsp2 of highly pathogenic porcine reproductive and respiratory syndrome virus emerging in China is not related to its virulence. J. Virol.  83(10):5156-67.
2. Han, J. et al. (2007) Identification of Nonessential Regions of the Nsp2 Replicase Protein of PRRSV Strain VR-2332 for Replication in Cell Culture. J. Virol.  81(18):9878-90.
3. Wikipedia: PRRSV.
4. Altscul, S.F. et al. (1999) Basic local alignment search tool. J. Mol. Biol. 215:403-410.
5. Thompson, J.D. et al. (1994) CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22(22):4673-4680.
6. Perl, http://www.perl.org
7. BioPerl, http://www.bioperl.org

Running X-Window application without a screen

Batch processing of beautifying phylogenetic trees
I prefer to use treedyn when I decorate many number of phylogenetic trees at  one time. In use on Linux desktop, my Perl script works well to convert from raw trees to beautified trees. However,  it was stopped at running of treedyn because Tcl of treedyn essentially required a DISPLAY of X-Window when the script was integrated into an Web application.
Xvfb (X virtual framebuffer) presents a virtual screen for user or program, which works only on memory not real video device. And xvfb-run.sh, a script of Xvfb, is an utility for command line program.

$ xvfb-run.sh my_script.pl -infile tree.phy -outfile tree.ps

Please see details of xvfb-run.sh  in Dascalita's very concise post.


Additional comments
Chevenet, the author of treedyn, and Christen recently published ScripTree that is more flexible and functional program for script developers. This ScripTree is also required to use Xvfb in case of background execution without a screen even the program presents built-in CGI interfaces.

References

1. F. Chevenet, http://www.treedyn.org/
2. Wikipedia, X virtual framebuffer 
3. D. Dascalita, Building a Flex project with Maven using Hudson on CentOS
4. F. Chevenet and R. Christen, http://phylo.lirmm.fr/scriptree/

스마트폰을 바꾸다. 삼성 블랙잭에서 구글의 넥서스원으로

애플의 아이폰4를 기다리다가 지쳐 결국은 안드로이드 폰으로 교체를 했습니다. 그전에 쓰던 스마트폰은 삼성의 블랙잭으로 MITS 사용자 포럼에서 WM6.2 롬을 업그레이드후에는 그럭저럭 쓸많했지만, 너무 힘겨운 웹 브라우징과 안스러운 카메라 기능 때문에 기종 변경을 결심했습니다.

잠시, 블랙잭에서 주로 사용하던 어플리케이션을 살펴보면,

1. 주소록 관리 : 구글 Contact 연동으로 800건 가까이 되는 연락처를 관리했습니다. 기본 주소록에서 초성 검색이 잘되는 편이고, 롬 업그레이드 이후에는 속도도 만족할 만 했습니다.
2. 전자우편 : 구글 Exchgane 서버를 이용한 Gmail연동을 주 메일 관리자로 사용했습니다. 쿼티자판이 있어 꽤 편리하고 빠르게 글을 작성할 수 있었습니다.
3. 트윗 : MoTweet을 이용했는데, 글 작성은 편리했지만 작은 스크린 때문에 가끔 긴 글을 읽을 때는 스크롤 경계에서 읽기 어려운 경우도 많았고, 사진을 업로드하면 아주아주 오래 걸리거나 먹통이 되는 경우가 많았습니다.
4. 페이스북 : 마이크로소프트에서 제공하는 Facebook 어플리케이션을 이용했습니다. 주로 읽기 전용이였고, 글 작성은 대부분 PC에서만 했습니다.
5. 음성녹음 : 지금 세 돌, 7개월 된 아이들의 목소리를 녹음하는데 사용했답니다. 아이들의 목소리를 그들이 컸을 때 들려주고 싶더군요. 전화벨 소리로 지정하는 것도 쉬워서 꽤 재미있는 전화벨 소리를 세팅할 수 도 있었답니다.
6. 검색 : 구글 어플리케이션을 설치했습니다. 미국 사이트로 들어가야 어플리케이션을 다운 받을 수 있었답니다. 급한 키워드 검색과 지도 검색은 꽤 요긴하게 사용했습니다.
7. 카메라 : 카메라 화질이 나빠서 가끔 사용했습니다. 동영상도 찍기는 했지만 PC로 다운 받는 것도 번거럽더군요.

이렇게 정리하다보니 나름 블랙잭을 십분 활용했다는 생각이 드는군요.  100메가 데이터에 대해서 매월 5000원의 정액요금을 냈지만 한번도 용량을 초과한 적은 없습니다. 이렇게 제대로 활용했지만 최소 3일에 한번씩은 오동작을 해서 재부팅을 해야하는 것과 풀브라우징을 못해서 이메일이나 소셜 네트워크의 링크를 확인할 수 없는 것은 너무 불편하더군요.

아무튼 이제는 넥서스원으로 전향을 하고 첫 사용후기를 쓴다면,

• 편리한 점
• PC 싱크가 불필요 : 전 음악이나 동영상을 다운 받아서 플레이하지 않기 때문에 스마트폰 사용과 관련된 모든 것이 넥서스원 상에서 다 해결되더군요. 물론 PC에 있는 미디어를 일부 옮기고 싶을 경우에는 USB를 연결해야 겠지만 제가 활용하는 범위에서는 거의 그럴 이유가 없을 같습니다. 이 건 정말 편리합니다.
• 카메라 : 다른 스마트폰과는 비교할 수 없지만, 예전 것에 비해서 너무 좋아 졌습니다. 그리고 무엇보다 사진을 찍은 다음에 바로 웹으로 전송할 수 있어서 관리하기 쉬워졌습니다.
• 소셜네트워크 : 다양한 어플리케이션이 있어서 입맛에 맞춰서 쓸 수 있겠더군요. 그리고, 다른 분들이 올린 링크를 풀브라우징할 수 있으니 예전의 블랙잭은 반쪽 자리였다는 생각이 들었습니다.
• 불편한 점
• 전화걸기 : 이건 기본 기능이라고 생각드는데, 넥서스원의 전화걸기 UI는 너무 생경합니다. 적응하는데 시간이 많이 걸릴 것 같네요.
• 터치키보드 : 블랙잭의 쿼티자판을 쓰다가 이 걸 써보니 너무 불편하고 적응이 어렵습니다.
• 배터리 : 착탈식이지만 구매시에 추가 배터리가 기본으로 제공되지 않습니다. 새로 구매해야하는데...

앞으로는 차근 차근 어플리케이션을 최적해 시켜보면서 활용도를 높여가야할 것 같네요.

말을 배우기 시작한 은채의 낱말 사전

30개월을 막 넘기시작한 첫 째 딸, 은채는 두 음절 이상은 어려워하는 것 같다. 요즘 주로 하는 말들과 그 뜻을 정리해본다. (언젠가 아이 크고나서 이런 날들을 잊고 싶지않아서 남기는 포스트..)

• 아깜~ : 아, 깜짝이야!
• 씨꺼 : 시끄러워 (주로 TV 소리가 클 때)
• 고마~ㅂ: 고맙습니다.
• 은채민 : 은채 비타민(약국에서 파는 뿡뿡이 비타민류를 말함)
• 피커 또는 티커: 스티커
• 알탕 : 사탕
• 초쿄 : 작고 검은색 종류의 먹을 수 있는 모든 것 (심지어 홍삼, 청국장환 등)
• 이야 : 은우. 아내와 내가 "은유야~"라고 부르는 것이 이렇게 들리는 모양이다.
• 태~ : 책
• 라멘 : 라면. 꽤 정확한 일본식 발음을 한다. ㅎㅎ
• 아이스미 : 아이스크림 (응용사례: 초코아이스미)