생명정보학을 공부하면서 궁금한게 있다면 이곳으로~!!

http://openbio.kr/questions/

혜식형님이 (시간이 남아돌아서..)만든 사이트

슬슬 글들이 올라오고 있습니다.

궁금하신것이나 알고 있는 정보들에 있어서 함께 공유하고
해결하면 좋을 것 같습니다. :)

RStudio IDE Server 설치기

RStudio IDE (Server)

잘 알려지고 널리 쓰여지고 있는 통계프로그램인 R을 
좀더 편리하게 사용할 수 있게 해주는 IDE 프로그램입니다.
오늘 말하려는것은 Desktop용이 아닌 Server용

처음에 이해한건 고성능 서버에 RStudio Server를 설치하고
데스크탑에는 RStudio Desktop라는 Client를 설치해서 사용하는 줄 알았는데..
하.... 그게 아니었더군요..;;;

걍 Desktop은 내 컴퓨터의 있는 R을 인식해서 좀더 사용하기 편하게
IDE로 보여주는 것으로 역할이 끝!!
-괜히 깔았어.... ㅎㅎ

여하튼...
R과 RStudio를 설치하면서...
좀 삽질을 하는 관계로..
R을 설치하는데 libR.so가 없어서 고생을..
그것도 개고생을... 별것도 아닌것이...

그래서 어쩔수 없이 현재 최신버전인 R-2.15.2는 사용못하고
rpm버전으로 제공되는 R-2.14.1을 사용해서 해결했다죠..
나중에 시간되면 R-2.15.2로 업그레이드 작업도 한번 할 예정입니다. ㅎㅎ

역시 설치할때 메뉴얼 숙지는 필수인데
저처럼 친절하지 않습니다. :)

R과 RStudio를 다운로드 받은 rpm으로  설치를 잘 마무리하고
실행시키는데 /usr/include 폴더 밑에 R 헤더파일이 없다고 투덜투덜투덜..

그래서 R 소스파일을 대충 컴파일해서 거기서 나오는 헤더파일을
해당 폴더에 샤샤샥 복사~ (컴파일한 파일은 제거~)

그리고 서버와 세션 설정 파일을 작성한후 구동시키면 끝~ :)
(example파일로 설정 파일을 제공하는 줄 알았는데;; 헐.... 걍 관리자가 파일 만들어서 하면 된다는... 이건 모.. ㅎㅎ)

아.. RStudio Server의 경우 웹서버를 구동하기 때문에
기존에 80포트를 사용하거나 추후에 웹서버를 운영할 게획이 있으신 분께서는
포트 변경해서 RStudio를 구동  시키시면 됩니다. :)

이렇게 위의 일련의 작업들을 정상적으로 잘 마무리하시고
URL란에 http://<server-ip>:<port>와 같이 입력하시고
엔터를 치시면 다음과 같은 화면을 웹브라우저를 통해서 확인하실 수 있습니다. :)

아.. user id와 passwd는 system account을 그대로 사용합니다. :)

RStudio (Server)의 장점은
고성능 서버의 자원을 사용해서
내 컴퓨터(모바일기기를 포함해서 웹브라우저가 되는 기기면)에서 작업하기 어려운
대규모(혹은 간단한) 통계작업이나 분석 작업을 진행할수 있다는 점~
그리고 plot, graph도 볼 수 있다는 점이라고 말씀드릴 수 있겠습니다. 고갱님~ :)

bam파일에서 fastq로 파일을 뽑을 수 있을까?

당연히 뽑을 수 있으니
글을 쓰는 것이겠지요? ㅎㅎㅎㅎ

그러나 원하는 서열이 bam파일에 있는 전체  서열이 아닌 한
약간의 작업을 해줘야 한 다는 것

현재 사용하고 있는 bam2fastq에 발등을 찍힌 관계로
align작업 후 얻어진 bam파일에서 곧바로 bam2fastq를 사용하여
 fastq를 뽑지 않고 있습니다.

약간 귀찮지만 다음 단계들을 거쳐서 뽑으면
본인이 원하는 서열들을 정확히 bam파일에서 뽑아 낼 수 있다는 것!!

bam2fastq나 그런 류의 프로그램만 사용하면 된다는 구글링 결과는
거짓부렁;;; 제길...

현재 다운로드 가능한 bam2fastq는 1.1.0 이다.
좋은 결과 있으시길~ :)


samtools view -H align.bam > align.mapped.sam
samtools view -F4 align.bam >> align.mapped.sam
samtools view -bS align.mapped.sam > align.mapped.bam
bam2fastq --aligned -o align#.mapped.fq align.mapped.bam

명령어 주석
-H는 헤더파일을 뽑는 옵션
-F4는 저도 정확히 모르겠지만 bam파일에서 -F4는 paired-end read가
모두 align되는 flag인듯 합니다.
-f4를 해서 저장한 파일들을 보면 align되지 않은 것들이 저장되는 것은 확인하였고,
-F4의 경우 align 정보가 표시되는 것으로 보아 맞는것으로 보입니다. :)
-F는 해당 flag를 제외한 결과를 return하는 옵션이고, 
-f는 해당 flag를 포함한 결과를 return하는 옵션입니다.
그러므로 -f4를하면 unmapped된 결과만 저장되고, -F4를 하면 unmapped되지 않은 결과가 저장되게 됩니다. :)
생성한 sam파일을 다시 bam파일로 변환하여 bam2fastq를 사용하여
fastq를 얻으면 됩니다. 다만, 구글링 결과에서 --no-aligned와 --aligned가 같다고
하는 글들이 있었는데..
--aligned를 해야 align된 paired read들만 fastq로 저장됩니다.
--no-aligned의 경우 결과가 상이한 것으로 나타나서 --aligned를 권장
--aligned와 동일한 결과를 보여주는 옵션은 --no-unaligned...
믿거나 말거나~ ㅎㅎ

MySQL에세 제공하는 스토리지 엔진

Posted at 2009/07/14


출처: afeleia

MySQL 스토리지 엔진에 대해서 좋은 글이 있어서 포스팅~ ^^;; ㅎㅎ
-생명공학으로 학부와 석사를 마쳤는데.. 왜.. MySQL을.. ㅋㅋ

현재 자신이 사용하는 스토리지 엔진을 확인하는 방법

mysql> show table status;
혹은
Information_schema 테이블 확인 (5.0이상부터 지원)

1. MyISAM
MySQL의 기본 스토리지 엔지으로 데이터 저장에 실제적인 제한이 없고 매우 효율적으로 저장한다. Full-Text 인덱스를 지원하며 특정 인덱스에 대해 메모리 캐쉬를 지원한다. 트랜잭션은 미지원/ 테이블 레벨의 락을 지원 잦은 변경및 삭제에는 좋은 성능이 나오지 못하나 데드락 발생은 예방

2. InnoDB
ACID 트랜잭션을 지원하며, MyISAM보다 데이터 저장비율이 낮고, 데이터 로드 속도가 느리다. 특정 데이터와 인덱스에 대해서 메모리 캐쉬를 지원하며 외부티를 지원한다. 데이터 압축이 불가능하고 자동 에러 복구 기능이 있다. 테이블 레벨이 아닌 ROW 레벨의 락을 지원한다.

3. Cluster (NDB)
트랜잭션을 지원하고 모든 데이터와 인덱스가 메모리에 존재하여 매우 빠른 데이터 로드 속도를 자랑하며 PK 사용시 최상의 속도를 나타낸다.

4. Archive
MySQL 5.0부터 새롭게 도입된 엔진으로 자동적으로 데이터 압축을 지원하며 다른 엔진에 비해 80% 저장공간 절약 효과를 자랑한다. 그리고 가장 빠른 데이터 로드 속도 또한 자랑하지만, INSERT와 SELECT만이 가능하다.

5. Federated
MySQL 5.0부터 새롭게 도입된 엔진으로 물리적 데이터베이스에 대한 논리적 데이터베이스를 생성하여 원격 데이터를 컨트롤 할 수 있다. 실행속도는 네트워크 요소에 따라 좌우되면 테이블 정의를 통한 SSL 보안 처리를 한다. 분산 데이터베이스 환경에 사용한다.

Convert a InnoDB table to MyISAM

Posted at 2008/08/25



기본 셋팅으로 테이블을 만들경우

InnoDB로 테이블 생성

InnoDB는 성능보다는 무결점(?)을 우선으로 하기때문에

성능이 MyISAM 타입보다 현저히 떨어진다고 합니다. ㅎㅎ


기존에 InnoDB를 MyISAM으로 바꿔주는 쿼리문.


ALTER TABLE tablename ENGINE = MyISAM;

Delete Command

Posted at 2008/08/25


DELETE

delete from table_name
해당 table의 내용을 삭제

ALTER
alter table table_name {change|add|drop} col_name col_name contents
table의 해당 컬럼을 변환하거나, 추가하거나, 삭제

MySQL UTF-8 셋팅

Posted at 2009/04/23

테스트 환경 Fedora 10 기본으로 깔리는 MySQL 5.0.77

참고 사이트
http://blog.artworker.biz/335
http://towis.net/2689923

vi /etc/my.cnf  [기본으로 설치되는 mysql의 경우]
-존재하는 섹션이 있고 존재하지 않는 섹션이 있다
 섹션이 없으면 추가하면된다.

[client]
default-character-set = utf8

[mysqld]
..생략..
default-character-set = utf8
default-collation=utf8_general_ci
init_connect=set collation_connection = utf8_general_ci
init_connect=set names utf8
character-set-server=utf8
collation-server=utf8_general_ci
character-set-client-handshake=TRUE
..생략..

[mysql]
default-character-set=utf8


파란색으로 Bold처리한 글씨는 새롭게 추가
나머지는 기존에 있는 것이었슴(Fedora 10, mysql 5.0.77의 경우)
-각 값들이 무엇을 하는 것들인지는 저도 잘..

환경 변수 수정 후 mysql 재시작
> /etc/init.d/mysql restart

mysql> show variables like 'c%';
+--------------------------+----------------------------+
| Variable_name            | Value                      |
+--------------------------+----------------------------+
| character_set_client     | utf8                       |
| character_set_connection | utf8                       |
| character_set_database   | utf8                       |
| character_set_filesystem | binary                     |
| character_set_results    | utf8                       |
| character_set_server     | utf8                       |
| character_set_system     | utf8                       |
| character_sets_dir       | /usr/share/mysql/charsets/ |
| collation_connection     | utf8_general_ci            |
| collation_database       | utf8_general_ci            |
| collation_server         | utf8_general_ci            |
| completion_type          | 0                          |
| concurrent_insert        | 1                          |
| connect_timeout          | 10                         |
+--------------------------+----------------------------+
14 rows in set (0.00 sec)

* MySql에서 데이터베이스 생성 
mysql>CREATE DATABASE {database_name} DEFAULT CHARACTER SET utf8 COLLATE utf8_general_ci;

대용량 HDD, mount하기

근래 본의아니게
컴퓨터 셋팅을 하게되서
간단한 작업 로그를 정리합니다.

대용량 HDD 포맷 후 mount하기.
Windows가 아니라 Linux 용입니다. ㅋ

요즘 단일 하드로 4TB가 나오는데..;;;
아직 그걸 사용하기에는 가격 매리트가 전혀 없는 관계로
(작업상 당연히 필요하지만 가난한 관계로 필요하다고 걍 살수없습니다.
 많은 랩들이 그러하듯이 ㅎㅎ :) )

anyway,
그래서  4T하드보다 가격 매리트가 아름다운 2T를 선호하고 있습니다.
2T보다 큰 파일들은 어떻게 할것이냐? 그래서 RAID를 사용하고 있지요 :)
그런데 공식적으로 windows나 linux에서 전통적으로 사용하는 파티션 프로그램인
fdisk는 2T 이상의 용량을 하나의 파티션으로 설정 할 경우 비추를 하고 있습니다.
왠지 저한테 묻지 마세요 저도 몰라요 ㅎㅎ

대신 parted라는 프로그램을 제공하고 있습니다.
물론 리눅스에서 GUI용 프로그램으로도 제공하고 있다고 합니다.
저는 CUI만 쓰니깐 몰라요~ :)

일단 fdisk와는 다르게 parted는 설정 후 'w: write'단계가 없습니다.
설정하면 설정되는 겁니다. 실시간으로 중간에 취소하기 없는겁니다.

RAID카드에서 2T이상의 용량으로 RAID를 잡았거나 4T 하드를 달았거나 동일합니다. :)
  #fdisk -l   
명령어를 사용하여 현재 시스템에 어떤 HDD들이 인식되어 있는지 확인합니다.
parted로 파티션을 나누고 포맷할 장치를 확인합니다. :)
예) /dev/sdb 장치를 설정해야 한다고 한다면....

#=======================================
#parted /dev/sdb
(parted) mklabel gpt //대용량 하드 형식인 gpt를 사용해서 라벨링을 하겠다라는 의미
(parted) mkpart //파티션을 만들겠다는 명령어

     partition name [primary]? [Enter] //파티션 이름으로  걍 엔터
     File system type? [ext2]? [Enter] // 어차피 ext2로는 안할것이니 엔터
     Start? //파티션할 장치의 시작을 처음으로 잡을 경우 0을 기입하고 엔터
     End? //파티션의 마지막을 장치의 어느부분으로 할지 설정하는 단계, 전체를 잡을 경우 100%라고 하면 되고, 아닐 경우 원하는 숫자를 GB단위로 기입하고 엔터
(parted)q

#mkfs.ext4 /dev/sdb1

#mount /dev/sdb1 /data
#=======================================
이렇게 하면 /data 위치에 /dev/sdb1 마운트가 되서 사용가능해집니다. :)

Warning Warning BAM header...

Warning: BAM header too large File

TopHat - cufflinks 조합으로 RNAseq을 분석하는 분들 중에서
과연 얼마나 접해보셨을까 하는데요..
혹 수백,수천개의 chr을 가진 genome을 분석하시는 분께서는 보셨을지도..

그렇습니다. 아직 genome project가 완벽하게 완료되지않아
chr이 완벽하게 정리되지 않아 그렇게 길지 않은 scaffold로 존재하고 있는 경우
마주할 수 있는 문제입니다.
BAM header 파일이 너무 커서 즐!! 이라고 내뱉는것입죠

samtools view -H input.bam #bam파일 Header만 print하는 명령어

구글링을 통해 얻은 결론
소스코드 수정후 새로 컴파일을 해야한다는 것!
(아놔... precompile된것만 편하게 사용하고 있었는데...
 소스 컴파일한다고 더 제대로 작동한다는 보장도 없는데 말이죠 아놔;;; )

여하튼 cufflink 소스 파일과 필요한 패키지들을 (DNS가 문제인지 외부로 직접은 안붙고
내부만 붙어서 다른 서버에서 다운받아서 복사한 후) 어찌어찌해서 설치 ㅋㅋ
현재까지는 잘 작동있다는 점~
컴파일하는동안 내내 warning이 화면을 도배했다는 점~
이거마저 안되면 난 모르겠다는 점~


cufflinks 패키지 설치시 많은 난관들이 있었지만 구글링을 통해 해결
그 경험을 정리해서 필요한 것과 수정해야 하는것들을 순서대로 다시 정리하자면

1. cmake 설치
 cmake 다운로드
 생각하지도 말고 root로 접속하여
 압축 풀고 폴더 안으로 침투하여
 >./bootstrap
 >make && make install
 을 나비처럼 날아서 타이핑과 엔터를 치면 나도 모르는 사이 설치가 되고 있다는 사실!!


2. boost 설치
 boost 다운로드
 이것 역시 다운로드 후 압축 풀고 root권한으로 접속 하여 설치작업을 진행하는것이 여러모로 건강에 도움이 될 듯 하다.
>cd boost/tools/build/v2/
>./boostrap.sh
>./b2 install
>./bjam #심심하면 이것도 실행을... 아... 기억이.. ㅎㅎ
boost 설치에 대해서는 cufflinks 홈페이지를 참조하는 것도 나쁘지 않는듯
cufflinks 튜토리얼
 root 권한으로 걍 설치하면 BOOST_ROOT path 지정하는게 필요가 있을까 하는 생각도..


3. samtools 설치
 samtools 다운로드
이건 모 설치하는데 크게 어렵지 않는 관계로 걍 본인 계정으로 설치, 그냥 압축 풀고
make 실행시키면 설치될듯합니다.
다만 이후 head파일이나 library파일을 위에서 언급한 tutorial페이지에 나와있는대로
올바른 위치에 복사를 해주어야 정신건강에 좋을듯 하다는 말을 남기면서 다음 단계로 고고씽~!!!


4. eigen 설치
 eigen 다운로드
이것 역시 압축을 해제한 후 계정을 root권한으로 갈아타고
압축 해제한 폴더로 들어가서 서브 디렉토리중 하나인 Eigen 폴더를 통채로
시스템 헤더 파일이 있는 곳으로 복사하면 OK
복사할 디렉토리는 cufflinks 튜토리얼을 참고하시길..


5. cufflinks 설치
 cufflinks 다운로드
마지막으로 대미를 장식할 오늘의 주인공 cufflinks
이것은 꼭 root 권한으로 설치 안해줘도 상관없다.
각자 개인 계정에 압축을 풀고 설치 과정을 시작하면 된다.
단 이 글에서 문제가 되었던 위험 요소를 제거를 하기 위해서 파일하나를 수정할 필요가 있다. :)
>cd cuffinks-2.x.x
>vi src/hits.cpp
MAX_HEADER_LEN = 4 * 1024 * 1024 로 되어 있는 것을
각자 genome 사정에 맞게 수정하면 된다. 
본인의 경우 MAX_HEADER_LEN = 128 * 1024 * 1024
>./configure --prefix=/path/to/cufflinks/install
(위에서 root권한을 사용해서 기본설정으로 설치를 안해주었다면
 --with-boost/--with-eigen/--with-bam 경로를 설정해 주어야 한다.
>make && make install


5단계를 거치고 나서 cufflinks가 설치가 완료되면
이제 새로 컴파일한 실행 파일로 실행하면 일단 warning을 뱉어내지 않으면서
일을 시작할 것이다.

중간에 세그먼트 폴트 에러를 뱉어내면서 죽지 않기를 바랄뿐이다. :)

SFF 파일은 어떻게 작업을 해야하나...

아....

딱히 인연이 없던 454 파일을 작업할 기회가 생겨서..
다음과 같이 스크립트를 좀 작성했습니다.

454에서 제공하는 Data analysis를 이용하지 않아서 좀 거시기합니다.
(Homopolymer trimming은 제공하지 못하고 있습니다. ㅎㅎ)

convertsff.py 

SFF파일을 fastq 혹은 fasta, qual 파일로 변환하는 스크립트
fastq로 변환하는 경우 illumina 1.3/1.5+ score로 변환 됩니다.
biopython이 설치되어 있어야 합니다.


SFF_Filter.py

illumina 데이터와는 일단 길이가 차이가 나니.... ㅎㅎ
filtering 스크립트를 간단히 만들었습니다.
cutoff base quality와 cutoff read length는 사용자가 설정 할 수 있게 하였습니다. :)
그리고 N의 포함 정도와 cutoff base quality 포함 정도는 고갱님의 의견을 반영하여
박하지 않게 설정해서 fixed시켜 놨습니다.
맘대로 수정하셔도 무방합니다. :)
다만 더 좋은 옵션이나 방법으로 업데이트 하셨다면 공유를 해주시면 더더욱 감사드리겠습니다.

데이터 변환 및 필터링이 끝나고 나면
QC를 해봐야 겠죠?

좀 간단히 결과를 이쁘게 그려주는게 어디 없을까 하고 있었는데
prinseq라는 프로그램이 있어 잠시 사용해봤습니다.
사용방법은 어렵지 않아요~ :)


Base Quality가 Phred +33인 경우

perl prinseq-lite.pl -verbose -fastq <input.fq> -graph_data <output.gd> -out_good null -out_bad null

Base Quality가 Phred +64인 경우

perl prinseq-lite.pl -verbose -fastq <input.fq> -graph_data <output.gd> -phred64 -out_good null -out_bad null

Quality Check 결과물을 Html로 확인하는 경우
perl prinseq-graphs.pl -i <output.gd> -html_all -o <output_name>

Quality Check 결과물을 png로 확인하는 경우
perl prinseq-graphs.pl -i <output.gd> -png_all -o <output_name>





추후에 시간이되면
Data analysis 프로그램을 설치해서 작업하는 단계나 방법에 대해서 설명하도록 하겠습니다. 그리고 추가적으로 NGS QC toolkit을 사용해서 QC하는 것도..
S대 L모군이 찾아논건데 괜찮아보여서 테스트 해볼까 합니다.
Illumina 외에 454도 지원하고 제일 매력적인건 multi-thread를 지원한다는것!!

모 여하튼...

다음기회에~ :)

TopHat을 바라볼때 중요한것

Read Manual!!! 

TopHat manual


사실 알고리듬 모르니...
라고 생각한다면.. 모 어쩔수 없고?? ㅎㅎ :)

하지만 무엇인가 알고 돌리는것과 모르고 자연에 출판된 protocol만 따라 돌리는것에는
많은 차이가 있으니..

T사의 K박사님의 정보로 TopHat 2.0.5를 허벌나게 사용중에 있습니다.
-한달 전만해도 TopHat 2.0.4를 사용중에 있었습니다.
-그 석달 전?? 반년 전 만해도 TopHat 1.0.3?을 사용하고 있었다는...


여하튼...
이번에 TopHat 2.0.5를 사용하면서 기존과 다르게 사용한 옵션이 있으니

--read-realign-edit-dist


그리고 사용안한 옵션도 있으니

-G / --GTF

옵션 이름 만으로도 대충 감들 잡으셨을 테니 옵션에 대한 설명은 패스하고,
왜 -G/--GTF를 사용안하냐?
(엄밀히 말하자면 known gene과 prediction gene의 문제..)
이 옵션을 사용하게 되면 --read-realign-edit-dist를 active시킨 의미가 없어지기 때문입니다.

이번에 --read-realign-edit-dist를 사용하면서 running 시간이 dramatically하게 증가하는 것을 경험했는데, S대 L군의 말로는 자기는 running 시간이 차이가 많이 나지 않는 다는 것!!
둘의 차이가 모였냐하니.. -G옵션을 사용하고 안하고 차이였습니다.

-G 옵션 설명에 gtf 정보를 사용하여 transcript sequence를 뽑아내서 거기에다가만 mapping을 한다는 것;; (역시 지도 교수님은 위대하다는 ㅎㅎ, 본인의 경우 해당 페이지를 몇번을 보고도 그냥 지나쳤었는데.. ㅎㅎ)

여하튼... -G를 사용하고 --read-realign-edit-dist 옵션을 사용하는것도 의미가 있겠지만 -G를 사용하지 않는게 더 좋은 결과를 낼 수 있지 않을까하는 단상을 끄적여 봅니다.

각자 실험하는 개체에 따하 gtf 사용여부를 판단하시면 되고 어떤 결과를 보느냐에 따라
--read-realign-edit-dist를 사용 여부를 결정하시면 됩니다.

제 경우 이게 그냥 자연에 출판된 protocol에 나온 방법보다 좋을것 같다는 생각이 듭니다.
이제 조만간 결과가 나오니 확인해보고 다시 글을 쓰도록 하겠습니다.


그리고 아시다시피 TopHat을 돌렸으면 cufflink도 돌리셔야죠.. ㅎㅎ :)
(아님 말고 ㅎㅎㅎㅎ )


ps. 누누이 말하지만 Human/Mouse는 default와 자연에 출판된 protocol이 甲이 맞는듯 합니다. ㅎㅎ 

그렇게 좋은 PacBio에 손이 안가는 이유...

"진정 우리꺼는 여러분들에게 좋으면 좋지
해를 안끼친다는.... "

- PacBio 본사 시니어 연구원느님의 발표


그렇게해도 PacBio는 정이 안간다는 ㅎㅎ

Illumima/ Life Tech.는 "우리거 좋아, 한번 써봐" (라는 우리꺼 안쓰면 니네 좀 후회할껄?)라는 느낌이라면,

PacBio는 "이번 논문에도 나왔듯이 우리꺼쓰면 울트라 캡숑 짱 따봉 좋아요 한번 써보세요" (라는 느낌?)

점심먹으면서 K군과 담소를 나누면서
Microorganism/ Meta genome 분야에서는 454에 비해 확실히 경쟁력이 있는데
(미국에서 1K Fungal genome project에서 PacBio를 사용하고 있다고 합니다.)
그외에는 과연 얼마나 경쟁력이 있는지... 잘 모르겠다는.... ㅎㅎㅎㅎ

그리고 제일 중요한건,

개인적으로 PacBio를 선듯 사용하지 못하는 이유는
비용문제에 대해서 확실한 해결책을 제시하지 못하고 있다는것도 큰 문제인듯..

PacBio를 가장 괴롭히는 것이 Error ratio문제인데
어차피 random error니깐 depth가 많으면 된다는 점~

다만, 다른 시퀀서의 QV를 맞추기위해 그 depth만큼
시퀀싱을 하면 비용 증가로 이어진다는것.

지구상에 재료비에 제한을 두지 않고 풍족하게 사용가능한 랩을 제외하고
사용 가능한 QV에 맞는 depth만큼 시퀀싱할 랩 아니면 ㅎㄷㄷㄷ

모 어차피 시퀀싱 업체에 맡기면 되니깐~  :)

ps. 약간의 글 수정이 있었습니다.
ㄴㅈㅊ에 다니는 지인의 염려가 있어 약간 수정을 하였습니다.
기술적인 부분이 아닌 현실적인 문제인 비용문제에 대해서 언급했으니
모 문제가 있겠냐마는.. ㅎㅎㅎㅎ

왜 샘플수가 400이어야 했는가.

왜 300을 하지?? ㅋㅋ

사실 KOGO학회 이전에
창범형님께서 facebook 담벼락에 공지하지 않았다면
그냥 지나갔을 법한 일이었습니다.

사업 취지와 목표, 방법에 대해서 주저리주저리 작성해놓은 문서 보면서
그냥 들었을만하 의구심..
물론 학회장에서 피뽑는다고 해서 신기한 마음에 걍 했지만 서도.. ㅎㅎ :)

샘플수집이 걍 지원자 400명?
2배수 3배수 뽑고 선별이나 무작위 추첨해서 분석하는것도 괜찮을듯한데..
그리고 지역별로도 차이가 날텐데.. 흠..
그리고 400명은 어떻게 나온건데?

궁금한점이 몇가지 들었죠 ㅎㅎㅎㅎ

근데 생각하다보니 이거 구상한 분들도 나랑 똑같이 귀차니즘이구만?
하고 생각을 마무리하게 됐지요 ㅎㅎㅎㅎ

그 이유..

1. 샘플 수집
이런 류의 사업을 진행하때 샘플수집을 어떻게 하는지 잘 모르겠지만..
안내문을 보면 서울에서 400명의 지원자를 받아서 사업을 진행하는 것 같이 보였는데,
이런 사업을 걍 지원자만 가지고 진행하는건지.. 좀 궁금하고..
서울에서만 진행하는건지 좀.... 그랬는데..

서울과 수도권을 포함하면 대한민국 절반이니.. 쿨럭;;;
-이거 완전 서울지상주위자로 몰리겠는데;;; 제가 이번에 시민이되서 그런게 아니라구요 ㅋㅋ

2. 대표성
위에서 언급한 수집 방법과 연결된건데...
Koread reference genome사업인데
대한민국을 대표할수있는 유전체가 될려면...??
이북도 해야하는건거 아닌가? 라는 말이 나올수 있습니다. ㅎㅎ
이북을 포함하지 않은 대표서열을 만든다면
지역적안배를 고려해볼만했을텐데
이북을 포함해서 대표서열을 구축하고자 한다면 굳이 지역적으로 나눠서 할 필요가??
ㅎㅎㅎ 어차피 이북에서 샘플을 못구하는데..
이미 이남에 편향되어 있다고 볼수 있지 않나하는 생각이...
그래서 서울에서만 해도 크게 문제가 되지 않고 대표성에도 문제가 될까하는
생각이 들었습니다.
-물론 서울에서 하니 새터민 분들도 참여를 했냐라고 물으신다면...
 아니요 라고 말씀드릴수있습니다. 확률상 서울이 높다 이거죠;;
 그리고 그분들은 이런 과제 있다는거 모르고 사실겁니다.
 저도 창범형님아니었으면 그분들과 똑같이 모르고 지나갔을테죠 ㅎㅎ
ps. 일단 본 과제에 최소한 하나는 이북 서열입니다. (저요 ㅎㅎ)


3. 샘플 숫자
아... 이건 좀...
제 생각이 맞다면 진짜 안습인데;;;
과제란게 한정된 자원으로 진행하다보니.. 쿨럭..
아시겠죠? ㅋㅋㅋㅋ
이것저것 빼고나니 시퀀싱할 비용으로 할수있는 샘플수가 400명
언저리;;;
이렇게 생각하고 싶진 않지만...ㅋㅋㅋㅋ
그리고 지원자로만 채워질 숫자라면 400도 그리 쉬운숫자가 아니었을겁니다.
이바닥에서 놀고있는 저도 창범형님아니었으면 그냥 지나가고
나중에 저런것도 했었구나 라고 할판이었으니깐요;;;; ㅎㅎㅎㅎ


여하튼... 이번 사업으로
M/D/T사 하나는 NGS로 영업이익은 꽤나 날듯(순이익은 아님..  ㅋ)
시퀀서 규모로는 M사가 유력하긴한데..
그리고  추가적으로 K사도???
아니면 걍 국책과제에 클라우드 리소스 무상 대여;;;;
그냥 그렇다고요 쿨럭;;; ㅎㅎㅎㅎ

파일의 포맷을 변환하는데 필요한 것들

내가 아니란 말이닷!!! ㅋㅋ

python에서 Biopython을 이용하여
간단하게 convert하는 샘플 코드를 제공하고 있으니
여러분들도 쉽게 만들수 있어요~ :)
Biopython에서 제공하는 Tutorial 


오늘 문의가 들어온 파일은 sff파일
Roche의 454 GS FLX? sequencing 결과파일로....
ABI와 함께 illumina한테 밀려서 뒷방으로 들어앉은 파일 포맷입니다.
그러나 아직도 쓰는 이유는 read 길이가 길기때문 :)

그렇습니다. PacBio도 Nanopore다 디립다 길게 sequencing해준다는
애들이 있습니다. 그런데 왜 옛날꺼 쓰냐?? PacBio는 base quality가 안습이고,
Nanopore는.... 언제 출시일지 전 잘 모르겠습니다. 업자가 아닌관계로 ㅎㅎ

그래서 위의 길게 sequencing 해준다는 시퀀서를 제외하고는 Roche의 454가 read 길이가 가장 길다고 할 수 있겠습니다. NGS중에선 말이죠

그런데 sff파일을 보려고 하면 문제가 생깁니다.
권모씨께서 문의를 한것이 그것때문인지는 모르겠지만 걍 일반인이
sff파일을 걍 직접 볼수가 없습니다. 왜냐 binary파일이니깐요(sff파일이 binary라고
알고 있는데  직접 다뤄본적이 없어서... ㅎㅎ )

그래서 사람이 볼수 있게 파일을 변환시켜줘야 한다는 겁니다.

convertSff.py

#!/usr/bin/python

import os, sys
from Bio import SeqIO

try:
inputSFF = sys.argv[1]
outputPREFIX = sys.argv[2]

except:
print "Usage: python convertSFF <input.sff> <output_name>"
print ""
exit(1)


SeqIO.convert(inputSFF,"sff","%s.fasta"%(outputPREFIX), "fasta")
SeqIO.convert(inputSFF,"sff","%s.quality"%(outputPREFIX), "qual")
SeqIO.convert(inputSFF,"sff","%s.fastq"%(outputPREFIX), "fastq")

권모씨의 요청으로 급조한 날림 convert python 코드 ㅋㅋ
이 스크립트를 수행하면 세개의 파일이 나오게 될것으로 예상됩니다. ㅎㅎ
안나오면 어쩔수없고... ㅎㅎ


아.. 그리고 사족으로 LT사의 SOLiD의 경우 우리가 알고 있는 서열과 달리
첫 염기 서열만 서열이고 그 다음부터는 A/G/T/C 알파벳이 아닌 숫자로 되어있는데..
이걸 굳이 변환해서 reference geneome에 mapped 작업하지 말라고 합니다.
Re-sequencing하는 경우라면 변환해서 mapping하지 말고 원래 원본 파일 그대로를
input으로 하는 align 프로그램을 사용해서 mapped한 다음에 그 다음 작업을
일반적으로 사용하는 samtools나 GATK같은 프로그램을 사용하라고 합니다.
(다들 알고있는거 한번더 상기 시켜드렸습니다. 혹시 아나요 SOLiD 포맷을 분석하게 될지.. ㅎㅎ)

분석시 raw 파일을 사용해야 하는 이유는 SOLiD만의 월등한 quality 효과를 볼수 있어서
그러지 않겠나하는....  믿거나 말거나 저 혼자만의 생각입니다.. ㅎㅎ
다만, 타사 제품과 다르게 복잡하게 숫자로 표현한건 아니겠죠...
나름의 숨은 뜻이.... 쿨럭.. (설마... 간지용;;;;; )

Re-sequencing이 아닌 denovo일 경우 모 어쩔수 없이 fastq파일로 변환을 해야 하지 않을까 합니다. assembly 프로그램을 작동시키려면 아무래도 SOLiD format보다는 fastq 포맷이
수월하니깐요.. :)

그럼....

Illumina Adapter Sequence

Illumina Sequencing에서 사용되는
Adapter중 TruSeq (분석할때 받는 데이터들이 다 요녀석으로 되어 있어서...) DNA/RNA Adapter Sequence를 확인해서 확인해봤습니다. ㅎㅎ

Type	Sequence
TruSeq Universal Adapter	AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT
TruSeq Adapter, Index 1	GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 2	GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGATGTATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 3	GATCGGAAGAGCACACGTCTGAACTCCAGTCACTTAGGCATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 4	GATCGGAAGAGCACACGTCTGAACTCCAGTCACTGACCAATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 5	GATCGGAAGAGCACACGTCTGAACTCCAGTCACACAGTGATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 6	GATCGGAAGAGCACACGTCTGAACTCCAGTCACGCCAATATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 7	GATCGGAAGAGCACACGTCTGAACTCCAGTCACCAGATCATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 8	GATCGGAAGAGCACACGTCTGAACTCCAGTCACACTTGAATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 9	GATCGGAAGAGCACACGTCTGAACTCCAGTCACGATCAGATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 10	GATCGGAAGAGCACACGTCTGAACTCCAGTCACTAGCTTATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 11	GATCGGAAGAGCACACGTCTGAACTCCAGTCACGGCTACATCTCGTATGCCGTCTTCTGCTTG
TruSeq Adapter, Index 12	GATCGGAAGAGCACACGTCTGAACTCCAGTCACCTTGTAATCTCGTATGCCGTCTTCTGCTTG

나중에 급할때 찾기 좀 애매해서...
요기다가 급 정리 ㅎㅎ

데이터 받았는데 TruSeq DNA/RNA Adapter인데 Index 12번보다 큰 경우
TruSeq Small RNA Index를 사용하는 것이라고 하네요
기본적인 Adapter sequence는 TruSeq DNA/RNA 인데 Index만 TruSeq Small RNA..
모 그렇다고 합니다. :)

좀더 자세한 Illumina Adapter Sequence에 대해서 알고 싶다면
다음 링크 참조 LINK

Velvet 사용 Tip

작년부터 NGS, NGS 해서 한두번쯤은
많은 NGS 프로그램들을 들어보셨을 것이고
조금이라도 관련되어 있는 일을 하시는 분이라면
이미 많이사용해 보셨을 것이라고 생각됩니다.

오늘은 그중에서 Velvet, sequence assembly 프로그램에 대해서..

Velvet은 de Bruijn graph 알고리듬을 사용했다고 합니다(라고 메뉴얼에 적혀 있습니다.
물어보지 마십시요. 저도 몰라요 ㅠ.ㅜ ).

여하튼 assembly 프로그램은 제가 알고 있는 바로는 크게 두가지
velvet과 같이 graph 알고리듬을 사용하는 프로그램과
아니면 sequence overlap 방법을 사용하는 프로그램이 있습니다.

그중에서 전 velvet이 좋은데 왜냐?
사용하기가 간편해서 입니다. ABySS나 ALLPATH-LG와 같이 복잡미묘하지 않기때문이죠
(간단하기로는 BGI의 SOAPdenovo가 최고봉인듯;;; )
단, velvet의 경우 big size의 genome을 assembly 할때는 고성능 CPU도 중요하지만
대용량 RAM도 함께 준비하시기 바랍니다.
물론 Big Size genome이라서 RAM이 많이 필요 한게 아니라 assembly read 개수가 많기 때문에
대용량 RAM이 필요한것이니.. ㅎㅎ assembly에 사용할 read 양이 적으면 상관없습니다.
다만 assembly에 사용할  read양이 많으면 대용량 RAM은 필수라는 것 잊지마시기 바랍니다.
-RAM이 작아도 작동은 하지만 swap을 사용하므로 결과 보시기 힘드실겁니다. ㅋ

곰팡이를 assembly 해본 경험으로 미루어 Velvet 에서의 RAM 요구량을
대략 적어봤습니다.

Total Input (Giga)	Memory (Giga)
40	256
20	128
10	64
5	32
2	16
1	8

Velvet으로 assembly할때는 일단 RAM은 많으면 많을 수록 좋습니다. :) ㅎㅎ

그럼 오늘 포스팅은 여기까지 ㅋ

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

awk 사용법

Blast를 수행 시 -m 8을 하면 자료 뽑아내기가 쉬운 것을 알고 있다.
근데 리눅스에서는 awk라는 명령어를 사용하여 별도의 코딩을 하지 않고
일정 값 이상/이하의 결과들만 골라서 볼 수 있다.
-biodb의 wiki에 정리해둔 것이 있었는데...;;;

구글 뒤져보니 본인이 원하는 기능의 awk 기능을 잘 설명해 놓은 것이 있어서
그대로 옮겨보도록 한다. ㅋ

출처: :+:하늘을 닮은 호수:+:

-m 옵셥에서 8값을 선택한 결과 파일 (blastout.file) 에서 score가 100 이상의 결과들만 뽑길 원하는 경우
> awk '$12 > 100 { print $2 }' blastout.file

본 문장을 응용하면 결과값에 특정 문자열만 들어가 있는 것을 포함/제외 하고
출력하기, 결과값이 중복된 것이 있으면 sort를 이용하여 제거 할 수도 있으니
참 간편하지요? ;;;

그걸 몰랐던 학부 시절 때 하나하나 python으로 삽질하던 기억을 하면;; 아놔;;;
그러나 요즘도 걍 python으로 작업한다는 훈훈한 이야기가 전해내려 온다능..;;

CPAN이용해서 Perl 모듈 설치하기

Circos를 설치하기 위해서 perl 모듈 설치에 좀더
편하게 하는 방법이 있어서.... ㅎㅎ

CPAN (Comprehensive Perl Archive Network)는 각종 perl 모듈이
모아져 있는 사이트인데 perl 5.8이후부터 cpan이라는 명령어를 지원하여
보다 편리하게 perl 모듈들을 설치할 수 있게 되었다.
(윈도우의 cgwin은 역시 잘 안된다능;;;; )

5.8 이전의 경우
$perl -MCPAN -eshell
또는
5.8 이후
$cpan

이라고 명령어를 실행시키면 처음에 실행시키면 이것저거 무엇을 한다고 하고
 [yes]를 입력하면 알아서 작업을 한후 cpan이라는 프롬프트를 보여줍니다.
cpan>install 모듈명
으로 쉽게 설치 우후훗...

HTTP Status 상태 코드별로 페이지 작성하기

출처:http://linux.tini4u.net/stories.php

일반적으로 웹호스팅을 받아보면 HTTP Status code page가
윈도우의 기본적인 그것이 아니라 깔끔하게 제작된 페이지를 본적이 있을 것입니다.
그것은 아파치에서(혹은 사용자가 [Override enabled]) 서버에러에 대한
응답을 지정해주어서 그렇습니다.

아파치 문서에 보면 ErrorDocument 부분이 정의되어 있으니
한번쯤 읽어보시면 이해하시기가 쉬우실 겁니다.
대략적으로 정리하면 총 3가지 방법으로 응답을 지정할 수가 있는데,
그것은 아래와 같습니다.


1. 보통의 텍스트

ErrorDocument 500 "The server made a boo boo."

문자열인 경우엔 " " 안에 문자열을 넣으면 됩니다.
추가=> [" "] 표시는 텍스트임을 알려주는 것으로서 그 자체는 출력되지 않습니다.


2. 내부 전환

ErrorDocument 404 /missing.html

서버 내부의 파일로 전환하는 방법인데 주의할 점으로는
절대 경로로 지정했을 경우엔 틀린 방법이라는 것입니다.
문법의 최상위 경로(/)는 DocumentRoot를 의미하기 때문입니다.
즉, 내 도메인이 /home/foobar/www/로 DocumentRoot가 설정되어 있고
전환할 페이지가 /home/foobar/www/missing.html 에 있다면

ErrorDocument 404 /missing.html

위와같이 설정을 해주셔야 작동을 합니다.

다만 이것이 사용자단이 아닌 아파치단에서(httpd.conf) 설정되었을 경우
모든 도메인에 대해서 404 코드는 /missing.html 파일을 찾게 됩니다.
이런 경우 사용자 입장에서는 /missing.html 파일을 사용하지 못합니다.

만약 서버 관리자로써 모든 Status Code 페이지를 전환하려면

Alias /Error "/usr/local/apache/htdocs/Error/"
ErrorDocument 404 /Error/missing.html

이런식으로 Alias를 만들어 주면 해결이 됩니다.
왜냐하면 모든 도메인은 /Error/missing.html를 찾을 것이고
/Error은 /usr/local/apache/htdocs/Error/ 으로 보내지기 때문입니다.
/Error/missing.html를 풀이해보면
/usr/local/apache/htdocs/Error/missing.html 이 되겠죠.

추가=> 스크립트나 SSI로도 내부 전환이 가능합니다.


3. 외부 전환

ErrorDocument 402 http://www.example.com/subscription_info.html

서버 외부의 파일로 전환하는 방법인데 김정균님 경험상으로
외부 전환시 CGI와 htaccess 인증시에 505 Status Code가 발생했다고 합니다.
원래 요청과 관련있는 환경 변수의 상당수가 스크립트에 전달되지 못한다는 점을 알고 있어야 합니다.


※ 설정하면서 주의해야 할 사항
위의 설정대로 설정했고 아무런 문제도 없는데
IE 자체의 에러페이지가 보이는 경우가 있습니다.
이런 경우는 대략 몇가지 문제가 있는데 대표적인것은 아래와 같습니다.

1. 완전한 HTML 페이지가 아닐 경우
<BODY>로 시작해서 </BODY>로 확실하게 끝나지 않았거나
혹은 <HTML>로 시작해서 </HTML>로 확실하게 끝나지 않았을 경우 입니다.
또는 PHP의 exit 명령으로 종료할때 </BODY></HTML>가 출력되지 않는
경우에도 IE 기본 페이지가 나옵니다.

2. Custom error page의 사이즈가 기준치보다 작은 경우
Custom error page 라는게 크기가 정해져 있기 때문에
이 크기가 넘지 않는 경우에도 IE 기본 페이지가 나옵니다.
Custom error page 의 각 크기 기준치는 아래와 같습니다.

Code Description File Size


400 Bad Request > 512 bytes
403 Forbidden > 256 bytes
404 Not Found > 512 bytes
405 Method Not Allowed > 256 bytes
406 Not Acceptable > 512 bytes
408 Request Time-out > 512 bytes
409 Conflict > 512 bytes
410 Gone > 256 bytes
500 Internal Server Error > 512 bytes
501 Not Implemented > 512 bytes
505 HTTP Version Not Supported > 512 bytes

Fasta 형식의 파일에서 빈서열 제거

가끔씩 Blast를 수행하고자 formatdb를 수행 할 때,
다음과 같은 에러를 접한 적이 있으리라 본다.

[formatdb] WARNING: Cannot add sequence number XXXXX XX.XXX.XXX.
 because it has zero-length.
[formatdb] FATAL ERROR: Fatal error when adding sequence to BLAST database.

formatdb를 수행하려는 fasta 서열에
빈 서열을 가지고 있기 때문에 나오는 에러로
빈 서열을 제거하면 OK!

vi check.py

import glob,sys
from Bio import SeqIO


file = sys.argv[1]


ow = open(file.split('.')[0]+'.check','w')
seqs = SeqIO.parse(open(file), format='fasta')
 for rec in seqs:
        name = rec.description
        seq = rec.seq.tostring()


        if len(seq.strip()) != 0:
                ow.write('>'+name+'\n')
                ow.write(seq+'\n')


ow.close()

[test]$python check.py test.fasta

Finding Homo Polymeric

출처: Finding homopolymer stretches in contigs

필요에 의해서 perl로 있는거 이용해서 약간 수정

$ vi findingHomopolymer.pl

$file = $ARGV[0];
open(data,$file);
$s = <data>;
$min = 4;


while ( $s =~ /(A{$min,}|T{$min,}|G{$min,}|C{$min,})/g) {
   $end = pos($s);
   $start = $end - length($1) + 1;
   print "$start, $end, $1 \n";
}

$ perl findingHomopolymer.pl seq.fa

※주위: seq.fa에는 서열만 한줄로 있어야 작동합니다.
fasta 폼 이런거 인식 못합니다. ㅋ

MSSQL Data Type

숫자형
bit : 1bit, 0,1
tinyint : 1byte, 0~255
smallint : 2byte, -2^15 ~ 2^15-1
int : 4byte, -2^31~2^31-1
bigint : 8byte, -2^63~2^63-1
smallmoney : 4byte, -2^31~2^31-1
money : 8byte, -2^63~2^63-1
numeric /decimal : -10^38~10^38-1
float : -1.79E+308 ~ 1.79E+308
real : -3.40E+38 ~ 3.40E+38

문자형
char/varchar : 8000자 이하의 문자
nchar/nvarchar : 4000자 이하의 유니코드 문자
text : 8000자 넘는 문자
ntext : 8000자 넘는 유니코드 문자

날짜형 
smalldatetime : 4byte, 1900-1-1 ~ 2079-6-6, 1분 단위
datetime : 8byte, 1753-1-1 ~ 9999-12-31, 1 ms 단위

이진/특수형
binary/varbinary : 8000 byte 이하의 이진값
image : 8000 byte 이상의 그림등의 바이너리 값
cutsor : 커서 이름을 변수로 처리
uniqueidentifer : 전역 유일 구분자
teimstamp : 데이터베이스 안에서 유일한 수

Eclipse Plugin 설치

(지금은 또 어떻게 얼마나 많이 변했을지 모르는 Eclipse.... ;; )


Eclipse 사용시 사용하고자 하는 plugin을 설치할 일이 있을 것이다

웬만하면 [HELP] - [Software Updates]에서 해결 가능하다.
그러나 가끔씩 맘에 안들게 이 메뉴로 해결이 불가능 할 때가 있다.

그럴때에는 manual하게 설치를 해줘야 한다.
-모 그냥 심심풀이로 사용하고자 하는 경우에는 굳이
 스트레스 받으면서 할 필요 없다.

3.4버전 Eclipse인 Ganymede는 Europa와 달리
Plugin폴더에 파일만 복사하면 plugin이 설치가 안된다.
-그래서 조낸 힘들었다.. ㅠ.ㅜ
-Europa도 안써봤었는데.. 제길..
앞으로 만날 Maven, Ant가 무섭다.. xml 설정 같은거 지랄같이 못하는데..

여하튼... Ganymede에서 수동으로 plugin을 설치하려면
두개의 설정 파일과 두개의 폴더에 관련 파일들을 복사해 주어야
Ganymede가 기분좋게 인식해준다.
-[HELP]-[Software Updates]에서 의존성 문제로 설치안되던 녀석들도
너무 깔끔하게 설치된다는 사실.. 제길...

일단 수정되어야 할 파일
eclipse/artifacts.xml
eclipse/configuration/org.eclipse.equinox.simpleconfigurator/bundles.info

그리고 수동으로 설치 할 plugin관련 파일들을 저장할 폴더 두곳
eclipse/features/
eclipse/plugins/

그런데 문제는 각 파일과 폴더를 들여다 보면 막막할것이다.
파일안에 어떻게 내용을 넣어줘야하며, 폴더에는 어떤 파일들을 넣어줘야 하는지..

그래서 본좌는 개발용으로 사용하는 eclipse외에 버전별로(왜 버전별인지 플러그인 설치하다가 당해보면 알것이다.) 다운로드 받아놨다. ^^

그래서 원하는 plugin을 설치가 되는 eclipse에 설치 된 후,
그 eclipse에 저장된 폴더들과 파일의 관련 부분만을 긇어서 원래 개발용
eclipse에 첨가시켜주면 OK!!
젠장.. 이거 깨닫는데 한달 걸렸다..;;;

오픈소스 라이선스 가이드

ㅋ 이것마저 어렵군..
GPL
CPL
Apache License
MIT License 도 있고.. ㅋㅋ

어떤놈은 코드공개,
어떤놈은 공개안하고 돈받고 팔아도 되고...

그러나 일단 배포를 안하면 코드공개와는 상관없다는...

오라클 데이터 타입

문자형 데이터
CHAR : 고정길이 문자형 데이터 타입 (~2k Byte)
VARCHAR2 : 가변길이 문자형 데이터 타입 (~4k Byte)
NCHAR : 고정길이 유니코드 문자형 데이터 타입 (~2k Byte)
NVARCHAR2 : 가변길이 유니코드 문자형 데이터 타입 (~4k Byte)

숫자형 데이터
BINARY_FLAOT : 32비트 부동 소수 (4 Byte)
BINARY_DOUBLE : 64비트 부동 소수 (8 Byte)
NUMBER : 가변 숫자 타입 (~21 Byte)


날짜형 데이터
DATE : 고정 길이 날짜,시간 데이터 (~7 Byte, DD-MM-YY)
INTERVAL YEAR TO MONTH
INTERVAL DAY TO SECOND
TIMESTAMP : DATE 보다 정밀도가 높음
TIMESTAMP WITH TIME ZONE : 지역 시간대 적용크리에이티브 커먼즈 라이센스

hmmer Manual

Blast와 함께 보편적으로 사용되는 Hmmer에 대한 설명서
hmmbuild/ hmmcalibrate/ hmmsearch에 대해서 설명
-물론 제가 사용하는 옵션에 대해서만 blast만큼 많지 않음. default로 사용해도 문제가 없으니깐~ 문제를 모르는것일 수도.. ㅎㅎ


hmmbuild: hmm matrix 만들어 줌
hmmbuild [-options] <hmmfile output> <alignment file>

-F 기존에 동일 이름의 hmm파일이 있으면 삭제하고 새로 만듬. 이 옵션 설정 안해주면 hmmbuild 아예 실행안됨.

ex) hmmbuild -F your_file.hmm your_file.aln

-f/ -g/ -s algorithm styles을 설정하는 옵션 이번에 사용하면서 이런 옵션을 처음 봤습니다. 왠지 hmm 멋져보이는 이유는.. ㅋ

ex) hmmbuild -f your_file.hmm your_file.aln

--amino/ --nucleic 강제로 alignment file이 어떤 서열인지 알려주는 것입니다.

ex) hmmbuild --amino your_file.hmm your_file.aln

-sequence weighting strategies
- model construction strategies
위의 무엇인가 고급스러운 것을 최대한 안건드리면 사용하는게
제 생활신조입니다. default인 이유는 그런 이유가 있을 것이다 라는.. ㅋ
개인적으로 잘 아시는 분만 선택해서 사용하시면 됩니다.
사용방법은 옵션을 그냥 적어주시면 됩니다.

ex) Alternative model construction strategies중 --fast 옵션 사용
      hmmbuild --fast your_file.hmm your_file.aln

hmmcalibrate: 만들어진 hmm matrix를 보정 시켜줌
hmmcalibrate [-options] <hmmfile>
--cpu: 프로그램 수행에 사용할 cpu 갯수 설정, 멀티 코어의 경우 가능. 단, 컴파일 및 바이너리 파일을 받을때 cpu옵션이 on 되어 있는 것을 받아야 사용 가능

--seed: hmmcalibrate를 몇번 수행할것인지 설정 하는 옵션 인듯.

본인의 hmmcalibrate 사용 예

ex) hmmcalibrate your_file.hmm

hmmsearch: 만들어진 hmm 파일을 이용해서 유사한 서열을 찾음.
hmmsearch [-options] <hmmfile> <sequence file or database>

-A <n>: 상위 n개 까지만 출력
-E <x>: blast의 e-value cutoff와 같은 것
-T/ -Z옵션도 안좋은 값을 짤라내기 위한 옵션

--cpu : 프로그램 수행에 사용할 cpu 갯수 설정, 멀티 코어의 경우 가능. 단, 컴파일 및 바이너리 파일을 받을때 cpu옵션이 on 되어 있는 것을 받아야 사용 가능

 --domE <x> / --domT <x>
위의 -T/ -Z의 옵션과 같이 도메인에서 필터링 하는 옵션인듯. 사용 안해봤음. ^^
<sequence file or database>는 fasta format 파일이면 사용 가능함.

본인이 으레 쓰는 방법임. hmmsearch 결과는 '>'로 빼주면 됨.

ex) hmmsearch -E 0.001 your_file.hmm your_database.fasta > result.output

Blastall Manual

NCBI에서 제공되는 Blastall에 대한 메뉴얼
blast-2.2.18을 기준으로 작성합니다. 현재 2.2.20이 나와있죠??
아마 옵션은 거의 동일할것입니다.
제가 많이 사용하는 것을 중심으로 설명합니다.
-지금은 더 업되어 있을 겁니다.
 그리고 이제는 슬슬 BLAST+로 옮겨타보려고 계획중입니다. :)

-p 5개의 기본 blast 프로그램중 하나를 선택하는 옵션

ex) -p {blastn|blastp|blastx|tblastn|tblastx}

-d blast를 돌리기 위한 데이터베이스 선택하는 옵션

ex) -d {nr|nt|your_database_file}

blast에서 데이터베이스로 사용하기 위해서는 fasta파일을 formatdb로 blast에 사용할 수 있는 데이터베이스로 변환시켜주어야 사용 가능. formatdb 수행후 붙는 확장자 명은 적어주지 않아도 됨. 파일이름 적음.

-i 검색해보고 싶은 서열(들) 입니다. Query 파일은 fasta format으로 되어있어야 함.

ex) -i your_query_file.seq 현재폴더에 있는 서열 파일      
      -i /your/home/path/query.fasta 다른 폴더에 있는 서열 파일

-e Expectation value를 정해줘서 설정된 값보다 크면 결과에 포함시키지 않는 옵션. 일반적으로 blastn의 경우 1e-06/1e-12, blastp의 경우 1e-03/1e-06으로 설정하고 상황마다 조정하면서 사용.

ex) -e 1e-06

-m 결과 파일을 저장할때의 format 결정 옵션. 일반적으로 로컬에서 blast를 돌리시려는 분들은 대량의 서열을 분석하기 위함이니, -m 8이 결과 파일을 분석하기 용이함,

ex) -m 8

-o Blast 결과 파일 설정하는 옵션

ex) -o your_output_file

-M blast를 실행시킬때 Matrix를 사용하게 하는 옵션. 서열과 서열을 비교하면서 weight를 주어서 peptide 서열을 검색할때 사용됨. 기본값은 BLOSUM62.
Matrix는 /your_blast_folder/data/ 밑에 있음.

ex) -M {BLOSUM62|PAM250|your_matrix}

-a CPU가 1개 이상일때 blast 수행시 하나 이상의 cpu를 사용하게 하는 옵션

ex) -a 2