Velvet의 무서움...

예전에 Velvet 사용시 input 대비 사용할 만한 메모리를
표시한 포스팅이 있었는데...
input 서열양이 많아지니...
일정 이상에서는.... 그 두배정도가 필요할듯 보이네요..

raw데이터 (멀티 라이브러리) 압축 푼 데이터가 10G정도인데
필터링 하면 10G이하로 줄어들었을텐데...
72G 서버에서 헐떡데는 꼴이라니.... 쩝..

예상보다 시간이 많이 걸려서 왜 그러나 보니깐..
중간에 swap으로 가고 있는 상황이라..

velvet으로 어셈블리할때에는 일단 메모리가 무한한 곳에서... ㅎㅎ


그리고 velvet은 매번 할때마다 데이터 결과가 다르다는... ㅋㅋ
동일 버전에서 K-mer 및 다른 옵션을 동일하게 두번 돌렸을때 결과값이
달라지니 혹시 같은 옵션으로 돌린결과를 확인 할때 통계값이 상이해도
놀라시지 않으셔도 됩니다.

Ubuntu에서 cummeRbund 설치시 주의 사항

Ubuntu에서 RNAseq 분석 후 비주얼라이제이션 관련해서
(저는 사용하고 있습니다. ㅋ) 사용하고 있으신
cummeRbund 패키지를 설치하실라 치면 XML 에러가 생기는 것을 확인 하실 수 있습니다.

XML관련 라이브러리가 Ubuntu에 없어서 그렇다고 하네요.. ㅎㅎ

Ubuntu 10대에서 사용하던 방법인데 13에서도 먹힘니다. :)

참고 사이트 R-help

또한 XML과 함께 RCurl설치시 에러가 생기는 경우도 비슷합니다. :)
다음과 같이 라이브러리를 설치하시면 아름다운 설치 결과를 보실  수 있으십니다. ㅎㅎ

> sudo apt-get install libxml2-dev

> sudo apt-get install libcurl4-openssl-dev

잠시 헤매고 있었는데...
구글에 검색하니.... 해결 방법이 뙁~!!!

GFF3에서 유전자 개수가 몇개인지 궁금할때?

요즘 de novo를 다루는 관계로
assembly 후 gene prediction 할 때 지난번에 포스팅 했던 maker를 사용하는 일이
빈번하다.

maker 결과 중 gff3 type (이 gff/gtf 파일의 형식이.. 버전마다 상이해서... 물론 본인은 차이점은 잘 모르겠다는게 문제.. 여하튼 다르다고 하니...)으로도 파일이 생성되는데
이 파일을 분석에 사용하시라고 분석자에게 보내드렸는데..
안타깝게도 gff 파일이 처음이셨던듯하다.
그런 분에게 gff파일을 보낸 내가 잘못했지만...
gff파일에서 유전자개수를 잘못 알고 계신 관계로.. ㅋㅋ
(지금까지 그렇게 알고 계시면 큰 낭패인데...)

여하튼..
gff파일에서 유전자 개수를 세시는데
$wc genome.gff
하신 듯.. (다르게 하면 그 숫자가 안나오고 wc하면 언급한 숫자가 나온다)

그래서 간단하나마 gff 파일에서 유전자 개수 세기를
언급하고자 한다.
대충 숫자만을 알고 싶다면 굳이 스크립트 필요없다.
$cut -f 3 genome.gff | grep gene | wc

자 이러면 유전자 개수를 알 수 있다.

다음부터는 wc만 하지 않길 바라는 간절한 마음뿐...

JGI Project List

DOE산하 JGI의 프로젝트 페이지를 가보면
헐... 이라고 보일수도 있습니다만...
2014년에 새 프로젝트 레폿 시스템을 구축중에 있다니깐.. ㅎㅎ

제가 오늘 포스팅할 내용은
JGI에서 진행하고 있는 프로젝트들은 어떤게 있는지..
걍 한번 R로 그래프 그려봤습니다.
좀더 이쁜 그래프들을 많이 그릴수 있을것 같은데..
실력이 미천한지라.. ㅎㅎ

음.. 잘 보일지는 모르겠지만...
(JGI에서 export된 파일을 받아보면 아시겠지만 공란이 은근 많습니다.)
DOE산하 연구소답게 동물보다는 확실히 미생물이나 곰팡이/식물 프로젝트를
많이 진행하고 있는듯 합니다.
(archaea의 오타는... ㅋ 데이터가 원래 저러다 보니.. ㅋ)
곰팡이가 좀 증가한것처럼 보이는데..
저건 아마도 F1000의 영향인것 같네요..

위의 그래프는 JGI 프로젝트들 중 곰팡이 관련 프로젝트가  현재 어떤 상태로
있는지 보여주는 그래프 입니다. 예전에 시작했던 프로젝트들은 역시나 완료되어있습니다.(중간에 살짝 미확인 결과도 있긴 합니다만....)

다음에는 좀더 이쁘게 그래프를 그릴수 있기를 기대하며... :)

Regular Expression of Python

아...  간만에 re모듈 사용;;;
ㅋㅋㅋㅋ

역시.. 정규표현식은 녹녹치 않다는

ENSEMBL에서 제공하는 newick포맷을
MEGA에서도 import하게 해주는 변환시켜주는 스크립트를
re 모듈을 사용해서 구현..
이라고 해봤자 지저분한건 마찬가지...

re.split(r"\[\&\&NHX:[\w+=\w+:,\w+=\w+.\w+]+T=\d{4,10}\]",tree_string)

다른 re모듈대신 re.split를 사용한 이유는
tree정보가 한줄에 저장되어 있는 관계로... :)
그리고 match되는 pattern을 취하는것이 아니라 버리는 것이므로 split하면
pattern을 구분자로 리스트로 만들어 주므로 굳이 따로 작업을 할 필요가 없다는.. :)

다음을 해석해 보자면..
re.split(r"\[\&\&NHX:[\w+=\w+:,\w+=\w+.\w+:]+T=\d{4,10}\]",tree_string)


녹색: pattern이 대괄호([])로 쌓여져 있다는 것을 확인
빨강: 대괄호 이후 &&NHX: 문자열이 있다는 것을 확인
오렌지색: 한개이상의 문자열 "=" 한개이상의 문자열 그리고 ":" 혹은 한개 이상의 문자열 "=" 한개 이상의 문자열 "." 한개이상의 문자열 그리고 ":"으로  이루어진 문자열이 한번이상  있다는것 을 확인
노란색: 그 이후 T=로 시작하는 4자리에서 10자리의 숫자가 있다는 것을 확인

import re 

f_open = (open(file_name,"r")).readline()
ns = re.split(r"\[\&\&NHX:[\w+=\w+:,\w+=\w+.\w+:]+T=\d{4,10}\]",f_open)
print "".join(ns)

genewise 설치 관련 Tip

출처: 9 by 6


최신 버전의 genewise인 2.4버전대인 경우는 source를 컴파일해야 하는데
이게 잘 안될때가 있다는 점...

src폴더에 들어가서 make all하면
conflicting type for 'getline'이라는 에러가 계속 떠서 찾아봤더니
아주 좋은 해결방법이... ㅋㅋ
getline이라는 함수가 getline_new로 바뀐듯.. ㅋㅋ

위의 블로그에 나와있듯이..

sed -i.old 's/getline/getline_new/' HMMer2/sqio.c
sed -i.old 's/isnumber/isdigit/' models/phasemodel.c

phasemodel은 상관없는데 같이 묶어놔서.. getline처럼 함수가 바뀌어서
에러가 나는 경우인듯...

모 이렇게 해주면...
착하게도 에러없이 컴파일이 잘 되고
genewise 2.4.x를 사용하실 수 있습니다. :)

Maker란

Maker는 Gene annotation 작업을 하는 pipeline으로 EVM과 함께 많이 사용된다고 합니다.

요즘같이 자고일어나면 DNA sequencing 가격이 계속 떨어지는 세상에서는 많은 연구자들이 de novo sequencing을 하여 생명체의 genome을 확보하기가 몇년전과 비교해보더라도 확연하게 쉬워진것을 알 수 있습니다.

그래서 이런 gene annotation tool들이 필요해졌죠
genome sequence만 있어서는 알수 있는게 별로 없으니깐요
생명체 안에서 일을하는 것은 단백질이고 그것을 만들 설계도는 gene이니
내가 sequencing해서 genome을 가지고 있다고 해서 연구 끝이 아니라는 얘기.. :)

근데 왜 EVM이 아니라 Maker를 언급하는걸까요?
걍 제가 써봤으니깐 언급한 겁니다. 다른 이유는 딱히 없습니다. ㅎㅎ :)

Maker의 경우 장점이라고 할 수 있는게
genome의 repeat masking을 pipeline에서 해준다는거 정도? 꼽을 수 있겠습니다. :)

그거 말고는 EVM이랑 비슷한듯 합니다.
Annotation 결과 품질이나 알고리즘면으로는 모...
알수가없으니..
단점은 홈페이지가 심심하면 다운된다는 정도?? ㅎㅎ

그럼 Maker를 믿을 수 있겠느냐?
그래서 한번 확인해 봤습니다.

중고등시절 들어봤을 플라나리아
그리고 애국가에도 나오는 소나무(종이 좀 다를듯합니다. ㅋ)
최재천 교수님께서 좋아하시는 개미 몇종.. 등등
GMOD 사이트를 방문하시면 확인 하실 수 있습니다.

다음에 기회가 된다면
좀더 경험을 해 본 다음에..
더 좋은글로 찾아뵙겠습니다. :)



ps. GMOD에서 NESCent라는 곳에서 매년 Gene Annotation 관련된 school이 열리는 듯 합니다. 2013년 써머스쿨은 지나갔고 관심있으시고 여력이 되신다면 한번 참석해보시는 것도 나쁘지 않을 듯 합니다 :)

SOAPaligner 결과에 Header넣고 SAM/BAM 변환

NGS alignment tool 중 BGI의 SOAPaligner가 있는데, 그 중 SAOP2 (SOAP3는 GPU로 인해 제한적으로 사용가능하기 때문에 SOAP2를 주력으로 사용하고 있다.) 작업 후 sam->bam으로 만드는 작업 중 문제에 대해서 얘기하고자 한다.

SOAPaligner 작업을 하면 sam/bam이 아닌 text파일이 생성된다.

그러나 이 text파일은 우리가 routine하게 사용하는 프로그램의 input으로는
사용할수가 없다는 것이 함정!!!

그래서 일련의 변환 작업을 필요로 한다.

오늘은 SOAP의 sam/bam 변환 작업에 대해서 알아보도록 하겠습니다.

이제와서 갑자기 SOAP을 얘기하냐고 물으신다면
요즘 SOAP을 가까이 하는 관계로 ㅎㅎㅎㅎ
팔은 안으로 굽지만 포스팅은 많이 사용하고 있는 tools 중심이니깐~ :)


첫번째 SOAPaligner는 알아서 잘 alignment 하도록 합니다.

이제 본격적인 변환 작업입니다.
다행히 BGI에서 soap output 형식을 sam 형식으로 변환시켜주는 perl 스크립트를 제공하고 있습니다.
그이름도 아름답도다 soap2sam.pl

근데 함정은 soap으로 alignment된 결과만 sam 형식으로 변환해준다는 점.
sam에 있어야할 header같은 정보는 무시한다는 것!!!


soap에는 없는 sam의 header파일 만드는 방법은
그냥 간단히 reference파일을 사용하여 header로 변경하는 script코드 작성하여도 되고

import os,sys

from Bio import SeqIO

try:

    input_fa = sys.argv[1]

except:

    print "python convert_ref_sam_header.py <ref> > <outpu.header>"

    exit(1)

print "@HD\tVN:1.0\tSO:unsorted"

for rec in SeqIO.parse(open(input_fa), format='fasta'):

    name = rec.description

    seq = rec.seq.tostring()

    name = (name.split()[0]).strip()

    print "@SQ\tSN:%s\tLN:%s" % (name,len(seq.strip()))

print "@PG\tID:soap2\tPN:soap2\tVN:2.2.1"

아니면

samtools faidx reference.fa

작업 후

cut -f 1,2 reference.fa.fai | awk '{print "@SQ\tSN:"$1"\tLN:"$2}' > sam.header

필요한 컬럼과 문자열만 조합하여 추출하여도 된다.
(이 경우 상하단에 @HD와 @PG 헤더를 추가해주어야 한다는 불편함이 있다. )


이런 작업을 거치면
soap을 sam을 거쳐 bam파일을 만들어서 다음 작업을 이어나갈수 있게 된다. :)


혹 fai파일을 만들기 위해 samtools faidx 작업을 수행하다가 에러가 발생하는 경우
fai파일을 만드려는 reference파일을 확인해보기 바란다.
fasta형식의 서열간 빈 line이 있는 경우 에러를 발생시킨다.
이런 경우 간단히 vi로 들어가서":g/^$/d"로 공란을 제거해주면 정상적으로 fai파일을 만들게 된다.

그리고 가장 쉬운 방법은 나중에 알려주는 법
(역시 메뉴얼을 잘 읽어야 한다는... 아놔....)

samtools view -bT reference.fa input.sam > output.bam

위에와 같이 뻘짓 필요없다는... ㅋㅋㅋㅋㅋㅋ ㅠ.ㅜ 아...;;; 제길..

Velvet으로 Assembly하기

결론은 de novo는 흠좀무라는... ㅋㅋ

species가 무엇이 되던간에 말입니다. ㅎㅎ


요즘 몇몇 곰팡이를 조립중에 있는데
SOAPdenovo, Velvet, SPAdes, SSPACE를 중점적으로 사용하고 있습니다.
ALLPATHS-LG는 read 준비하는 것까지는 작동을 하는데
그다음에 본 작업에서 error를 뱉고 죽어버리는 어처구니 없는 상황을 연출해주는 바람에
일단 pass했습니다.

여러 strategy으로 접근 중에 있는데...

아... 본 글은 그게 중요한게 아니라... (다음 기회에 언급 하도록 하겠습니다.)

Velvet assembly할 때 long insert library (Mate-Paired Library)를 한개만 지원하는 것 같아그 문제를 해결하가 위한 tip입니다.

해결 방법은 다음 링크에 나와 있습니다.


단 EMBOSS의 "revseq"를 사용하면 안습이 된다는 것 빼고는 말이죠~ :)

그래서 준비 했습니다.

#!/usr/bin/python
import os,sys
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC

def readread(s):
        return [s.readline(),s.readline(),s.readline(),s.readline()]


try:
        fq = sys.argv[1]

except:
        print "python command.py <input_fastq>"
        print ""
        exit(1)

fq_open = open(fq,"r")

fq_write = open("%s.rev" % (fq),"w")

s1read = readread(fq_open)

while(s1read[0]):
        h1, h2 = s1read[0].strip(), s1read[2].strip()
        s1, q1 = s1read[1].strip(), s1read[3].strip()

        rev = Seq(s1,IUPAC.unambiguous_dna).reverse_complement()
        qual = q1[::-1]

        fq_write.write("%s\n%s\n%s\n%s\n" % (h1,rev,h2,qual))

        s1read = readread(fq_open)

fq_write.close()

-위 소스는 python 2.6 이상에서 작동합니다. :)

velvet으로 조립 시 mate-paired Library가 여러개라면 서열을 변경 한 이후 사용해 보시기 바랍니다. :)

The Great Good Place는 폐쇄중;;

다름이 아니라 무료로 호스팅하고 있는 Hostinger 측에서
제가 사용하고 있는 계정에서 다수의 이메일이 발송되어
자신들의 IP가 블랙리스트에 등록될까 염려되어 제 계정을 막았다는
소식을 접했습니다.

그래서 문의 결과 어떻게 할 수 없다는 답변을 받고 일단
제 도메인들을 모두 이 gwlee.blogspot.kr 포워딩 시켜 놓은 상황입니다.

어차피 아직 공식 개장 하지 않아서 적어놓은 글은 많지 않았지만..
WP셋팅하고 블로깅을 하려고 했는데 이런 시련이 닥치네요 ㅎㅎ

조만간 다른 안식처를 마련하여 공개하도록 하겠습니다. :)