FASTA-Reader Warning

BLAST+의 makeblastdb를 사용하여
blast용 db를 만들고자 할 때
항상 에러없이 잘 만들어지지 않을 수 있다.
거기다 다음과 같은

Error: (1431.1) FASTA-Reader: Warning: FASTA-Reader: Ignoring FASTA modifier(s) found because the input was not expected to have any.

내맘에 안드니깐 무시하고 만든다
라는  에러를 내놓을때는.....

db를 만들고자 하는 fasta파일 각 서열의 header를 간단하게 바꿔보는것도 나쁘지 않은 방법일듯... :)

왜 header를 수정해줘야 하는가?
본인이 header 수정해줘서 정상적으로 blastdb를 만들었으니깐요? :)

NGS 물량공세 플랫폼 오픈

오늘 페북에 내 스팸 메일 제외하고
가장 핫한 이슈는 일루미나가 새로운 NGS기계를 내놓았고
그걸 마크로젠이 냉큼 구매했다는 이야기...

>Illumina Sequencing System Spec<

물량으로는 일루미나는 때려잡아도 못잡을듯..
결과가 TB단위가 나올줄이야... ㅋㅋ

그리고 HiSeq X Ten 이름에 있는 것과 같이
HiSeq X를 사려면 10대라는 최소 주문 수량을 만족해야 한다는...
"우리 일루미나 고객님들은 쪼잔하지 않아요"가 2014년도 일루미나 슬로건인가? ㅎㅎ

일년에 한두번씩 휴양지에서 세미나 개최해서 초정도 해주고 그런것 같다만...
이게 무슨 청첩장도 아니고 최소 주문 수량을...;;;; ㅎㅎ

여하튼...
일루미나 덕분에
Open the real hell gate... ㅋㅋ

CentOS에서 genewise 설치하기

지난번에 genewise 설치 관련 포스팅은
사용하지 않는 함수를 바꿔주는것으로 해결을 봤는데..

이번에는 잘 설치 되어있는 라이브러리 헤더파일을 못잡아서
컴파일이 제대로 되지 않는 문제가 발생하였는데.. ㅋㅋ

이것 역시 솔루션이 ㅎㅎ

해결책은 9 x 6 에..

makefile에 있는 옵션 중 glib 옵션을 명확하게 바꿔줘서 해결한다는... ㅎㅎ

find ./ -type f -name "makefile" -exec sed -i.old 's/glib-config --libs/pkg-config --libs glib-2.0/g' "{}" +;
find ./ -type f -name "makefile" -exec sed -i.old 's/glib-config --cflags/pkg-config --cflags glib-2.0/g' "{}" +;

여하튼 해결되서 다행 ~ ㅎㅎ :)

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년 써머스쿨은 지나갔고 관심있으시고 여력이 되신다면 한번 참석해보시는 것도 나쁘지 않을 듯 합니다 :)