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 작업을 위해선..
스크립트언어라도 배우는게 좋다는 점~
간단한 코드는 짤 수 있어야 한다는 점~
텍스트 파싱은 할 줄 알아야 한다는 점~

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.

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+ 간단 사용법에 대한건 마무리하는걸로..

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")

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)

모.. 요정도??

MEGA5 Usage

오늘은 간단히 MEGA (Molecular Evolutionary Genetics Analysis) 사용법 중

Multi Fasta Sequence가지고 Alignment하고 Phylogenetic tree를 그리는 것에 대해서 

간단히 알아보도록 하자. :)

MEGA 사이트에서 알아서 개인정보를 팔아서 다운로드를 받던지 주위에 이미 받아논

지인에게 달라고 해서 얻기를 바란다. 

일단 설치 후 실행 시키면 다음과 같은 화면을 접할 수 있다.

Alignment하고 싶은 파일을 Open하도록 하자.

위의 [Open A File/Session ...]을 클릭하게 되면 다음과 창이 뜨게되며 알맞은 파일을 선택하면 된다.

알맞은 파일을 선택 한 후 [열기]를 선택하면 다음과 같은 창이 뜨게 되는데 당황하지 말고 걍 [Align] 버튼을 클릭하면된다. :)

[Align]버튼을 클릭하면 다음과 같이 보여지게 된다.

이제 Align을 해보도록 하자. Alignment 프로그램 중에 Clustalw를 사용하도록 하자.

다음과 같은 [Alignment] -> [Align by ClustalW]를 클릭하면 된다.

만약 서열들을 선택해주지 않았다면 다음과 같은 경고창이 보이게 된다. 이때에는 [OK] 버튼을 클릭하면 된다.

이제 ClustalW를 실행시키기 위한 Parameter를 설정하게 되는데 사실 건드릴 일이 그렇게 있을까? 필요하면 기호에 맞게 수정해서 사용하시길...

Protein Weight Matrix만 잘 사용하면 크게 문제가 없을 것으로 보인다.

Parameter설정 후 [OK] 버튼을 클릭하면 다음과 같이 Alignment를 수행하게 된다.

단, 서열이 많을 수록 소요 시간은 기하급수적으로 늘어난다는 점만 주의하시길..

Alignment가 완료되면 다음과 같이 Align된 결과가 보여지게 된다.

MEGA에서는 MEGA만의 저장 format를 지원하는데 Alignment된 결과를 포함하여 저장하기 때문에 naming rule을 잘 정의해서 사용하면 동일한작업을 두번 할 필요가 없게 된다.
아래 그림은 Alignment 결과를 저장하는 MEGA Aln 포맷은 mas로 저장하는 그림이다.

Alignment를 하였다면 대게 Phylogenetic tree까지 그리고자 할 것이다.
그리기 위해서는 mas파일이 아닌 MEGA Format 파일이 필요하다. 이것은 Alignment를 수행 한 창에서 [Data] -> [Export Alignment] -> [MEGA Format]를 클릭하면  MEGA 포맷으로 저장할 수 있다.

드디어 Phylogenetic tree그리는 시간이다. :)
MEGA Main 화면에서 [Analysis] -> [Phylogeny] 선택 후 기호에 맞는, 상황에 맞는 방법을 선택하여 Tree를 생성하면된다. :)

위의 순서대로 Method를 선택 한 후 MEGA format의 파일을 열면 다음과 같은 Parameter 설정 창이 보이게 되된다. 여기서도 본인 분석에 맞게 설정하여 작업하면 트리가 나오게 된다.

그럼,

Good Luck. :)

Python에서 통계 작업하기

통계작업할 때 지금까지는 대부분 엑셀,
조금 데이터양이 많아지거나 고급 통계기법이 필요하면 R을 사용하였다.

근데 아쉽게도 R에서의 문자열 핸들링관련해서 지식이 민망하여
통계처리 중간에 문자열 조작이 필요한 경우 중간에 값을 텍스트로 빼서
파이썬에서 작업하고 다시 R로 가져가는 상당히 귀찮은 작업을 거쳤다.

그래서..
결국 미루다 미루다..

통계작업 중간에 문자열작업이나 다른 수정이 필요업는 단계까지
파이썬에서 작업을 마치는 방법을 터득하기로 했다.

근데.. 이게 파이썬에서 통계작업을 한번도 안해본지라
몇날 몇일 헤매고 있었는데

Scipy(파이썬 설치할때 기본적으로 설치하는 모듈인데... 쓰지는 않고 있었다)에서
통계 함수들을 제공한다니.. ;;;

Scipy Statistical Function

여하튼...
이것때문에.. 어느정도 고민 해결~ ㅎㅎ

Velvet과 Oases 설치

Velvet 다운로드

$tar zxf velvet_1.2.03.tgz
$cd velvet_1.2.03
$make

default 설정으로 설치하는 경우 그냥 make하면 되지만..
K-mer값을 큰 값으로 사용하고자 하는 경우나 multi-thread를 사용하고자 하는 경우
옵션값을 조절해줘야 한다.
아니면 make할 때 옵션을 설정해 주어도 가능하다. :)

$vi Makefile
...
MAXKMERLENGTH=99
CATEGORIES=16
...
OPENMP=1
LONGSEQUENCES=1
...


$make 

Oases 다운로드

$tar zxf oases_0.2.06.tgz
$cd oases_0.2.06
$make


Oases 또한 Velvet과 마찬가지로 default 설정으로 설치하고자 하는 경우에는 make를 해주면 되지만, velvet경로도 내부 서버 상황에 맞게 설정 해야 하는 경우 가 있으므로, 다음과 같이 Makefile을 수정 후 make를 하면된다. :)

$vi Makefile
...
MAXKMERLENGTH=99
CATEGORIES=16
...
OPENMP=1
LONGSEQUENCES=1
VELVET_DIR=../velvet_1.2.03
...


$make 

tabix vs intersectBed

tab으로 구분되어 있는 txt 파일 처리하는 tabix에서 대해서
혜식 형님께 조언을 구하는데 대량으로 처리하려면
tabix보단 intersectBed가 더 좋다는 정보 획득!! ㅎㅎ

intersectBed 홈페이지

COUNTIF 함수

아... 어제 작업하면서 이런 기능 엑셀에 이런기능이 설마 없을까?
했는데... 역시 있었다. 착한 놈들 ㅎㅎㅎㅎ

이름도 기특하다 ! COUNTIF !
cell을 세는데 특정 조건에 맞는 셀 개수를 세는 함수

모 사용방법은 그렇듯이 항상 간단하다

=COUNTIF(범위, 조건)
범위는 우리가 자주 쓰는 A2:A30
조건은 개수를 세려는 셀의 조건이나 문자열 등등등..
내 경우 0.01보다 작은것이 몇개나 있는지 확인하기 위해서 이녀석을 찾았다.

그러므로 내가 사용한것은 =COUNTIF(A1:A10000, "<0.01")