You can write your love letter in……DNA.

1. Save the above scripts like ‘translate.py’
2. change permission

chmod +x translate.py

3. compose your love letter or something in ‘plain english’ and save it as text. Let’s assume you wrote something like this

I love you so much, babe..

4. Then execute the script.

./translate.py love.txt
>Iloveyousomuchbabe
ATGATTTTATAGGTTGAATATTAGTGATCTTAGATGTGATGTCATGTGGCTGTGGAATAA

Now you have a functional ‘gene’ for your ‘love letter’!. Even it will be translated in the cell cause it has start and stop codon.

GIST in wordpress.com

이전까지 소스코드를 블로그에 적을떄 <pre>를 쓰거나 아니면

태그를 썼었는데, 알고보니 wordpress.com에서는 gist 를 지원하고 있었다. ;;;

gist에 코드쪼가리를 만든다음, 워드프레스에서 글 작성하면서 그냥 gist URL 쳐넣으면
Screen Shot 2012-12-29 at 2.57.43 AM

gist 코드가 요기 임베딩되어있넹? 뭔코드인지는 묻지마셈 걍 흔한 awk 파싱해서 MySQL 쳐넣처넣

Multifasta에서 빠르게 일부 레코드만 추출해오기

이런 작업은 사실 일상적으로 하는 작업인 관계로 Perl이나 Python 과 같은 스크립트 언어를 조금이라도 아는 사람이라면 그닥 어렵지 않게 할 수 있는 일이긴 하지만, 이런 작업도 수GB 이상의 파일이라면 그리고 불러올 파일의 갯수가 늘어나면 적지 않은 시간이 소요된다. 좀 더 빠르게 이런 작업을 마칠 수 없을까?

Profile HMM 에 의해 도메인 및 단백질 검색을 하는 유명한 소프트웨어인 HMMER 3.0 에 딸려있는 유틸리티인 esl-sfetch를 이용하면 아주 빠르게 이런 작업을 할 수 있다.

HMMER3 을 설치하고 hmmer3 의 실행파일에 PATH를 걸어둔 다음

esl-sfetch -h

하면

$ esl-sfetch -h
# esl-sfetch :: retrieve sequence(s) from a file
# Easel h3.0 (March 2010)
# Copyright (C) 2010 Howard Hughes Medical Institute.
# Freely distributed under the Janelia Farm Software License.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: esl-sfetch [options]          (one seq named )
Usage: esl-sfetch [options] -f   (all seqs in )
Usage: esl-sfetch [options] --index        (index )

 where general options are:
  -h     : help; show brief info on version and usage
  ....중략.

가령 지금 나는 vibrio.fasta라는 274MB짜리 fasta 파일이 있고, 이 파일은

$tail vibrio.fasta
SDGATWLDRRLIHGETADLAPTGFGQQVREAMDQRREHHIEQGDATRSRDSRVFYRRNLLAILREREVAG
VGSDMALSKGLPFRAATDGESVSGKFTGTVHLSSGKFAVVEKSHEFTLVPWRPIIDRQLGREVMGIVQGG
SVSWQLGRQRGLER
>gi|424025505|ref|ZP_17765144.1| drug resistance transporter, Bcr/CflA subfamily protein [Vibrio cholerae HC-62B1]
MTTTRPAWAYTLPAALLLMAPFDILASLAMDIYLPVVPAMPGILNTTPAMIQLTLSLYMVMLGVGQVIFG
PLSDRIGRRPILLAGATAFVIASLGAAWSSTAPAFVAFRLLQAVGASAMLVATFATVRDVYANRPEGVVI
YGLFSSMLAFVPALGPIAGALIGEFLGWQAIFITLAILAMLALLNAGFRWHETRPLDQVKTRRSVLPIFA
SPAFWVYTVGFSAGMGTFFVFFSTAPRVLIGQAEYSEIGFSFAFATVALVMIVTTRFAKSFVARWGIAGC
VARGMALLVCGAVLLGIGELYGSPSFLTFILPMWVVAVGIVFTVSVTANGALAEFDDIAGSAVAFYFCVQ
SLIVSIVGTLAVALLNGDTAWPVICYATAMAVLVSLGLVLLRLRGAATEKSPVV

요런 식으로 되어 있는 multifasta 파일이다. 여기서 gi|424025505|ref|ZP_17765144.1| 와 같은 accession 으로 파일을 찾고싶다고 치자.

gi|260362785|ref|ZP_05775654.1|
gi|260880072|ref|ZP_05892427.1|
gi|28901226|ref|NP_800881.1|
gi|153835952|ref|ZP_01988619.1|
gi|308095182|ref|ZP_07663192.1|
gi|254224936|ref|ZP_04918551.1|
gi|229529188|ref|ZP_04418578.1|
gi|254285269|ref|ZP_04960234.1|
gi|421351525|ref|ZP_15801890.1|
gi|229515173|ref|ZP_04404633.1|

요런 내용을 list.txt로 저장했다고 치자
제일 먼저 할 일은 fasta 파일을 인덱싱하여 빨리 검색을 할 수 있게 하는거.

$esl-sfetch --index vibrio.fasta
Creating SSI index for vibrio.fasta...    done.
Indexed 665311 sequences (665311 names).
SSI index written to file vibrio.fasta.ssi

이러면 파일이름.ssi 식의 인덱스 파일이 생기고 이제 검색을 할 수 있다.

$esl-sfetch -f vibrio.fasta list.txt

하면 바로 쫙 stdout 으로 파일이 추출된다. ㅋ

사실 이런 작업은 hmmer 등을 수행해서 나온 결과파일을 가지고 해당하는 fasta 파일을 데이터베이스에서 추출할때 유용한데, hmmer를 수행해서 테이블 형식으로 다음과 같은 출력파일이 나왔다고 하자.
Screen Shot 2012-12-21 at 12.59.44 AM

이 출력물을 awk 등으로 파싱해서 원하는 gi 값만 esl-sfetch에 전달해 주면 곧바로 hmmer hit 에 상응하는 시퀀스 파일을 뽑을 수 있다.

awk '{if (substr($0,0,1)!="#") {print $1}}' VopL.fasta.hmmer | esl-sfetch -f vibrio.fasta - > hmmerresults.fasta

참고로 samtools 에도 faidx 라는 동일한 기능 (fasta 를 indexing한후 추출) 이 있으나 테스트해 본 결과 esl-sftech가 더 빠르다.
약 247MB의 fasta에서 10개 정도의 fasta를 추출할때 시간을 재보면

time awk '{if (substr($0,0,1)!="#") {print $1}}' VopL.fasta.hmmer | xargs samtools faidx vibrio.fasta
real	0m1.033s
user	0m0.921s
sys	0m0.110s

time awk '{if (substr($0,0,1)!="#") {print $1}}' VopL.fasta.hmmer | esl-sfetch -f vibrio.fasta
real	0m0.007s
user	0m0.005s
sys	0m0.006s

수십 GB 정도의 파일이라면 좀 더 차이날 것 같다.;;

-C 옵션을 이용하면 multifasta 내의 subsequence을 잘라오는 것도 가능하다. (BED file 등을 파싱해서 빠르게 시퀀스를 얻는것도 가능..아 그건 bedtools쓰면 되나..여튼.

어쩄든 NGS 시대가 되면서 처리해야 할 데이터의 양이 부쩍 늘어난 상황에서 파싱과 같은 간단한 작업의 경우에도 가능하면 unix 기반의 텍스트 툴및 파이프, 리다이렉션을 이용하는 것이 수행속도 면에서도 훨씬 이점이 있는 것 같다는.

1+1은..

요즘 파이썬으로 작성한 스크립트에 최소한의 UI를 씌워보고자 몇가지를 살펴보고 있는데, 요즘 유행하는 Microframework 을 이용하여 간단히 웹 기반의 인터페이스를 씌우는 게 가장 편한 일인듯.

Microframework 의 경우에도 여러가지가 있겠지만, 단 한개의 파일로 해결되는 bottle 을 써볼 예정. 테스트용으로 만든 코드..;;

#!/usr/bin/python
# -*- coding: utf-8 -*-
from bottle import route, run, template
import webbrowser

@route('/:name')
def index(name):
    if int(name) < 6:

        list = u"{0}+{0}은".format(name)
        i = int(name)
        for k in range(i):
            list = list + u"귀요";
        list = list + u"미";
    else:
        list = u"그만해;;;;지겹다";
    return template('<b>{{name}}</b>!', name=list)
webbrowser.open('http://localhost:8080/1')
run(host='localhost', port=8080)

이 파일을 run.py 뭐 이런식으로 저장해서 디렉토리에 넣고, bottle.py 를 다운로드받아 같은 디렉토리에 넣고, 터미널에서

./run.py

로 실행하면…
Screen Shot 2012-12-19 at 12.18.58 AM
…..;;;;;;;

http://localhost:8080/2 로 주소를 바꿔보면..

Screen Shot 2012-12-19 at 12.20.38 AM

6+6은 뭘까? 직접 알아보도록 하고 (쿨럭) 여튼 이렇게 하면 간단히 파일 하나로 웹어플리케이션 구축가능. ㅋ 스크립트 형태로 있는 것에 웹 UI를 달아야 할때와 같을때 매우 유용할 듯 싶다. 특히 twitter bootstrap와 같은 템플릿과 같이 이용하면 최소한의 노력으로 그럴싸한 외관의 웹 UI를 씌울 수 있게 될듯.

Bioawk

Heng Li (author of bwa, samTools..) made a awk version supplemented with bioinoformatic centric functions.
https://github.com/lh3/bioawk
Some examples
Reverse Complement of FASTA

awk -c fastx '{print ">"$name;print revcomp($seq)}' sequence.fasta
Create FASTA from SAM
   samtools view aln.bam | \
awk -c sam '{s=$seq; if(and($flag, 16)) {s=revcomp($seq)} print ">"$qname"\n"s}'

Get the mean Phred quality score from FASTQ:

   awk -c fastx '{print ">"$name;print meanqual($qual)}' seq.fq.gz

Detailed documentation is here.

PacBio로 Gap Filling 하기

오늘 소개할 논문은 아래 논문

 

English et al., Mind the Gap: Upgrading Genomes with Pacific Biosciences RS Long-Read Sequencing Technology, Plos One 2012

 

NGS 시대가 왔다고 하지만, 결국 모델생물이 아닌 잡생물의 지놈 어셈블리 퀄리티는 피니싱까지 마친 모델생물의 지놈 퀄리티에 결코 미치지 못하는 것이 사실. 게다가 short read에 의한 de novo assembly 의 경우 나오는 퀄리티가 기존의 생거 시절의 드래프트 퀄리티보다도 훨씬 못함.

그래서 이것을 좀 개선해 보려고 여러가지 시도가 되고 있는데, 역시 가장 유망한 것은 PacBio나 Oxford Nanopore와 같은 3세대 시퀀싱 기술이 일반화되어 몇kb씩 되는 read를 쑥쑥 뽑아내는 것일것임. 그러나 PacBio는 기껏해야 85% 정도의 base accuracy 를 가지고 있다는 문제, Oxford Nanopore는 뭐 기계가 나와야지 원…ㅋㅋㅋ

그래서 그나마 현재 데이터가 나오는 PacBio 데이터를 좀 보완하는 방법들이 몇가지 발표되고 있는데 이전에 발표된 것은 PacBio 의 long read에 illumina의 short read 데이터를 align 해서 에러교정을 하고, 이 long read 로 assembly를 하는 방법.

Image

그러나 저 위 링크에 나온 방법은 기존의 드래프트 어셈블리를 그대로 두고 여기에 PacBio read 를 얼라인한후 gap 을 메꾸는 넘을 찾아서 Gap 부분만 로컬 어셈블리해서 (한마디로 PacBio 데이터를 gap에 잘라붙이는;;;) 어셈블리를 개선하는 방법. 새롭게 어셈블리를 하지 않고 기존의 어셈블리에 근거해서 어셈블리를 개선하는 것의 장점이라면 어노테이션을 새로 홀라당 다 할 필요가 없다라는 이야기. (어셈블리가 틀려지면 contig 번호, base 넘버가 다 틀려지므로 기존의 유전자 위치정보가 다 도루묵이 됨 ㅋ)

Image

그냥 간단하게 말해서 gap 사이 통과하는 롱다리 PacBio read 골라내고 잘라붙이기 한다고 이해하면 되겠음.

그래서 이런 파이프라인을 거치면 얼마나 어셈블리가 개선되나?

Image

샘플에 따라서 틀리지만 어쨌든 6X 정도의 PacBio  데이터를 이용해서 잘라붙이기 신공을 발휘하면 90%에서 60% 정도의 gap을 메꿀 수 있다고..

물론 몇가지 개선되어야 할 부분은 있는데,

1. 위 파이프라인의 경우 scaffold 내의 gap (sequence gap) 을 메울 수는 있지만 scaffold 간의  gap (physical gap) 을 메우는 것은 현재로써는 불가능한듯. 사실 sequence gap 보다 scaffold 간의 link 를 확립하는 것이 보더 더 유용한 정보를 제공하는데 이 부분에 대해서는 아직 ㅋ 뭐 다음 페이퍼라구염? ㅋㅋㅋ

2. 아무런 에러교정 없이 PacBio 데이터를 잘라붙여서 Gap을 메꾸는 셈이므로 다른 데이터가 없이 PacBio 데이터만 있는 부분이라면 PacBio 의 에러가 같이 스리슬쩍 들어오겠졈? 아놔 Gap 으로 NNNN 하는 것보다는 좀 틀려도 시퀀스가 있는게 더 낫지? 라고 생각할지 모르겠지만 어셈블리가 만들어진 다음 genbank 에 들어갈때 ‘이부분은 PacBio 데이터 밖에 읍음. 에러 있어도 몰라 ㅋ’ 하고 써있는것도 아니고..뭐 하긴 모든 genome sequence는 완벽한게 아니니깐 글타 치자..

3. 결국 Long Read Sequencing 기술의 발전이 더 있어야 이런 꾸질한 짓거리를 안할 수 있다는 말. Oxford Nanopore가 빨리 나오든지 해야지..

4. PacBio는 이제서야 이런 식으로 ‘쓸모가 전혀 없는 것은 아니다’ 라는 내용이 나오긴 하지만..주가를 보면 과연 이 회사 얼마나 버틸까에 대한 의구심이 드는 게 사실ㅋ

Image

 

쥬커버그 “저런 애들도 있는데 나보고 뭐라하는거임? ㅋㅋㅋ”