월요일, 8월 27, 2012

Tophat을 run할 때의 마음가짐

RNA-Seq 작업을 하면서 빈번하게 사용하는 Alignment tool로 TopHat을 꼽을 수 있다.
(나의 경우 그렇다. 아니면 말고.. 쳇~)

본인의 경우 대부분의 프로그램들의 default값을 사용하기 좋아라 하지만
최근 NGS관련 tool을 다루면서부터 default값은 신뢰하지 않기로 했다.
왜냐?

최근 각광받는 NGS 분석 tool들의 대부분의 default값들은 Human, Mouse같은 Model 종들에 대해서 적합한 것 들이지 내가 다루는 곰팡이나 식물은 전혀 Out of 안중이기 때문이다.

그래서 아주 죽을맛이다라는거다 ㅋㅋ
성능 짱 좋은 서버로 테스트 해보고 싶은 경우의 수를 모두 다 해보면 좋겠지만
논문내는건 시간싸움이다 보니 해보고 싶은 모든 경우에 대해서 테스트 못할 수 도 있다.

그래서 옵션 중에서 Key가 될만한 옵션들만 본인의 종에 맞게 조정해서 분석을 해야 그나마 시간 대비 분석 결과에 만족 할 수 있을 것으로 생각한다.

그 중 TopHat의 경우 intron-length를 분석하고자 하는 종에 맞춰서 값을 사용하기 바라는 바이다.
TopHat의 --max-intron-length의 경우 500,000bp인데 상식적으로 곰팡이 같은 종의 경우 한 유전자안에 500kbp짜리 intron이 있을리 만무하지 않겠는가?

그래서 이런 종 특이적인 정보를 사용하는 경우 본인이 분석하는 종을 대표할 수 있는 값을 사용하는 것이 보다 좋은 결과를 얻을 수 있을것이다.
(강릉 교육에서 들어서 요건 확인하고 한다는거.. ㅋㅋ)

사람이나 마우스 하는 분들은 걍 default 값 사용하면됩니다. (요건 좀 부럽습네다. ㅎㅎ)

아... intron길이 구하는건 스스로, 그걸 누가 매번 알려줄수는 없잖아~
구글링하면 어느정도 커버 할수 있을 자료 찾을 수 있습니다.
요즘 NGS때문에 denovo도 꽤나 하는듯 하니..
-대신 없으면 추가로 denovo하시면 될듯... 전략만 잘 짜면... 괜찮을듯한데.. ㅎㅎ


그래서 NGS 작업을 위해선..
스크립트언어라도 배우는게 좋다는 점~
간단한 코드는 짤 수 있어야 한다는 점~
텍스트 파싱은 할 줄 알아야 한다는 점~




화요일, 8월 21, 2012

Tophat2에서 libz.so.1 에러에 대처하는 우리들의 자세

RNA-Seq 작업을 하시는 분들의 경우
많은 분들께서 TopHat과 Cufflinks 조합으로 분석을 진행하리라 생각합니다.

본 글은 좀 old한 리눅스 시스템에서
TopHat 그것도 TopHat2의 바이너리를 사용하여 작업을 하실 때
libz.so.1 관련 에러가 나는 문제가 발생했을 때 대응 할 수 있게 해줍니다.
(경험치 +1)

기존 시스템에서 사용하고 있는 libz.so.1의 버전이 옛날것이라
이미 컴파일 되어 있는 Tophat의 바이너리파일에 저장되어 있는 정보랑 맞지 않아
발생 하는 것으로 보입니다.
fc12에서 TopHat-1.4.0에서는 전혀 문제가 없었는데..
fc12에서 TopHat2에서는 문제가 발생해버렸네요.
(그리고 fc14에서는 문제가 발생하지 않습니다.)

그러므로 다른 에러는 저도 모르겠습니다. ㅋ

/lib64/libz.so.1: no version information available

위의 에러를 만나시게 된다면
다음 링크에 있는 파일(fc14의 파일입니다.)을
리눅스의 /lib64/폴더 밑에 다운로드 받아 저장하시고,
링크를 새로 만들어 주시면 됩니다. :)


파일 다운로드 libz.so.1.2.5

원래 시스템에 있는 libz.so.1 링크는 삭제

>ln -s /lib64/libz.1.2.5 /lib64/libz.so.1

이렇게 하면 다음부터는 위의 libz.so.1 에러는 발생하지 않을 것입니다. :)

Good luck.



목요일, 7월 19, 2012

Local에서 BLAST+ 작업하기 x64

예전에 포스팅한 Local에서 BLAST 돌리기는 32bit 버전이었는데
요즘 64bit OS를 사용하고 있는 관계로 (본인 또한 다 64bit ㅎㅎ) 업데이트를 해보기로 한다.

일반적으로 Local이라 함은 데스크탑 즉, 윈도우 환경이 대다수일 것이라고 생각된다.
(물론 리눅스나 맥을 데스크탑으로 사용하시는 능력자분들도 있으시겠다.)

Blast를 윈도우 환경에서 작업하고 싶은 경우
NCBI 사이트에 가서 BLAST 프로그램을 다운받으면 된다.

Blast+ 64bit, Blast+ 32Bit, Blast 64Bit, Blast 32Bit

Blast와 Blast+의 차이는 엄청나다 Blast는 기존에 간단한 옵션과 사용방법을 그대로 유지하고 있지만 Blast+의 경우 드라마틱한 속도 개선과 성능이 향상 되었다(는 모르겠고 옵션과 사용방법은 확실하게 드라미틱하게 복잡해졌다)고 한다.

일단 위에서 본인의 OS Bit수에 맞는 Blast를 다운로드 받고, exe파일을 더블 클릭하여 압축을 해제한다. 단, 알수없는 많은 파일들이 눈앞에 펼쳐지는 것을 보기 싫다면 별도의 폴더를 만든 후 더블 클릭하시길..

더블클릭하여 압축을 해제하게 되면 bin, data, doc 폴더가 생성되게 된다.

윈도우에서 Blast작업은 "명령 프롬프트" 창에서 하던가 아니면 별도의 스크립트를 작성해서 실행 시킬 수 있다.

기존의 Blast와는 다르게 BLAST 프로그램들인 blastp, blastn등과 같은 프로그램들이 각각 분리되어졌다.




Step 1. BLAST DB 생성

사용예
makeblastdb -in <input_file.fasta> -input_type {asn1_bin|asn1_txt|blastdb|fasta|xml} -dbtype {nucl|prot} -parse_seqids -hash_index

추가적으로 masking 작업을 위한 masker 프로그램이 동봉되어 있다.
-mask_data {dustmasker|segmasker|windowmasker}
사용하실려면 사용하시길... :)



그리고 output되는 파일의 파일 용량이 1G보다 큰 경우 blastdb 파일이 여러개로 쪼개 질 수 있다. 이런 상황을 방지하기 위해서는 다음 옵션을 사용하여 output 파일의 파일당 최대 용량을 늘려주면 된다.

-max_file_sz <size>, size는 "1GB", "2GB" 이렇게 작성하면 된다.

예제
~/blast+/ncbi-blast-2.2.25+/bin/makeblastdb -in input.fa -input_type fasta -dbtype nucl -parse_seqids -hash_index -max-file-sz 4GB


--------------------------------------------------------------------------------
Building a new DB, current time: 07/18/2012 13:00:27
New DB name: input.fa
New DB title: input.fa
Sequence type: Nucleotide
Keep Linkouts: T
Keep MBits: T
Maximum file size: 4294967296B
Adding sequences from FASTA; added 211174 sequences in 151.743 seconds.
--------------------------------------------------------------------------------


위의 작업이 끝나면 다음과 같은 파일들이 생성된다.
input.fa.{nhd|nhi|nhr|nin|nog|nsd|nsi|nsq}


Step 2. BLAST 실행

예전의 BLAST와 달리 BLAST+에서는 blast 프로그램들이 모두 각각의 수행 파일로 존재한다. 사용되는 옵션은 대부분 유사하니 크게 걱정할 필요는 없다. 다만 parameter 이름이 약간 달라진것 제외하고는 :)

사용예

blastn -query <input_file> -db <blastdb_file> -out <output_file> -evalue <e-value> -outfmt {0..11} -num_threads <number thread>

예제

blastn -query query.fa -db database.fa -out output.txt -evalue 1 -outfmt 6 -num_threads 8

"database.fa 파일에 query.fa파일을 8개의 thread를 사용해서 blastn을 하는 작업으로 e-value가 1이하인 것만 저장하고 결과 파일은 output.txt파일에 저장한다. 그리고 결과 형식은 tabular형식으로 저장한다" 라는 의미를 담고 있는 명령어임. :)


모 그럼... 이정도로.. BLAST+ 간단 사용법에 대한건 마무리하는걸로..

금요일, 7월 13, 2012

BLAT and BLAST Output Format Fields

BLAST (-m 8)과 BLAT 결과를 보면 table 형식으로 되어 있는데 head들이 친절하게 설명되어 있는 것도있지만 대량분석할 때에는 귀찮아서 Header 정보를 제거하고 결과를 뽑아서 가끔씩 헷갈릴때가 있다. 본인은 그렇다.. :)

그래서 다시 정리를... 쿨럭.. ㅎㅎ


NCBI Blast Tabular output format fields

(Blast Head의 경우 부연 설명이 필요 없을 정도로 simple하다. :))
QueryIdSubjectId
Identity percent
AlnLength
mismatchCount
gapOpenCount
QueryStart
QueryEnd
SubjectStart
SubjectEnd
Evalue
bitScore

위의 링크된 사용자가 심플하게 parsing하는 예제를 함께 보여주고 있는데 
참고하면 좋을듯 :)


-Python 

for line in open(“myfile.blast”):
(queryId, subjectId, percIdentity, alnLength, mismatchCount, gapOpenCount, queryStart, queryEnd, subjectStart, subjectEnd, eVal, bitScore) = line.split(“\t”)


-Perl 

while (<>) {
($queryId, $subjectId, $percIdentity, $alnLength, $mismatchCount, $gapOpenCount, $queryStart, $queryEnd, $subjectStart, $subjectEnd, $eVal, $bitScore) = split(/\t/)
}





BLAT Spec

matches int unsigned , # Number of bases that match that aren't repeats
misMatches int unsigned ,  # Number of bases that don't match
repMatches int unsigned ,  # Number of bases that match but are part of repeats
nCount int unsigned ,      # Number of 'N' bases
qNumInsert int unsigned ,  # Number of inserts in query
qBaseInsert int unsigned , # Number of bases inserted in query
tNumInsert int unsigned ,  # Number of inserts in target
tBaseInsert int unsigned , # Number of bases inserted in target
strand char(2) ,           # + or - for query strand, optionally followed by + or – for target strand
qName varchar(255) ,       # Query sequence name
qSize int unsigned ,       # Query sequence size
qStart int unsigned ,      # Alignment start position in query
qEnd int unsigned ,        # Alignment end position in query
tName varchar(255) ,       # Target sequence name
tSize int unsigned ,       # Target sequence size
tStart int unsigned ,      # Alignment start position in target
tEnd int unsigned ,        # Alignment end position in target
blockCount int unsigned ,  # Number of blocks in alignment. A block contains no gaps.
blockSizes longblob ,      # Size of each block in a comma separated list
qStarts longblob ,         # Start of each block in query in a comma separated list
tStarts longblob ,         # Start of each block in target in a comma separated list


-Python

for line in open(“myfile.blat”): 
(matches, misMatches, repMatches, nCount, qNumInsert, qBaseInsert, tNumInsert, tBaseInsert, strand, qName, qSize, qStart, qEnd, tName, tSize, tStart, tEnd, blockCount, blockSizes, qStarts, tStarts) = lines.split("\t")



수요일, 5월 16, 2012

Python 설치 순서

1여년만에 업그레이드 기념 OS 재설치 후 필요한 프로그램 재설치 중에 있습니다.
오늘은 제가 주력으로 사용하는 스크립트 언어인 Python 설치 순서에 대해서 정리하는 포스팅을 하겠습니다.
(이 순서는 지극히 개인적인 생각으로 하는 것이오니 순서바뀐다고
 경찰차 출동안합니다. 쇠고랑 안찹니다. :) )

지금 사용하는 OS가 64비트이기 때문에 기존에 호환되는 라이브러리들이 잘 작동하지 않습니다.
특히나 Rpython의 경우 64bit 윈도우에서는 정상적으로 작동 잘 안되는 것같습니다.
-Rpython이 지원하는 R과 pythpn버전이 정확히 일치해야 작동하는 듯 보입니다.

제가 주로사용하는 라이브러리들은 모 정해져 있으니.. :)

- Python 2.6.6 (Python 2.6버전에서 msi파일은 2.6.6이 마지막입니다.)

- Base-x.x.x (Base-12.5.7)

- mxBase-x.x.x (egenix-mx-base-3.2.4)

- numpy-x.x.x (numpy-MKL-1.6.2rc1)

- scipy-x.x.x (scipy-0.10.1)

- Biopython-x.x.x (biopython-1.59)

- distribute-x.x.x (distribute-0.6.26)

- PIL-x.x.x (PIL-1.1.7)

- reportlab-x.x.x (reportlab-2.5)

모.. 요정도??