Local BLAST 설치 및 사용 (in OS X)

한참전에 윈도우 기준으로 로컬 PC에 BLAST를 설치하고 돌리는 글을 이전 블로그에 쓴 적이 있다. 그런데 그것도 벌써 한 9년 전의 일이고, BLAST 자체도 기본적으로 바뀌지는 않았지만 BLAST+ 라는 형태로 바뀌었다. 그래서 업데이트하는 느낌으로 OS X 에서 Local BLAST를 설치하고 돌리는 과정을 적어보고자 한다. (윈도우는 어떻게 해염 하려면 옛날 글 을 참조하셈)

1. 다운로드

요기에 가서 dmg 포맷의 설치파일을 받는다. 세상 많이 좋아졌다. dmg 포맷의 설치파일이 나오니..

제일 아랫넘이다.

2. 설치

패키지 선택하고

오오 GUI Installer! (그러면 뭐해 어차피 나중엔 다 커맨드 라인인데)

여튼 이렇게 하면 쨘~ 하고 설치를 끝낸다. 그런데 대체 어디에 설치를 해놓은 거냐.

3. 설정

이렇게 설치를 하면 기본적으로 /usr/local/ncbi/blast 내에 파일을 복사해 놓는다. 

옛날에 커맨드 라인 BLAST  (지금 현재 Legacy BLAST라고 이야기하는) 를 설치해 본 사람이라면 명령어가 좀 틀려졌다는 느낌을 받을 것이다. 즉 이전에는 blastall 이라는 커맨드에서 blastp, blastn 등을 지정했지만 이제는 별도의 명령어가 등장했다.

여튼 여기서 추가로 설정을 해주어야 한다. blast 실행파일이 있는 /usr/local/ncbi/blast/bin에 대한 PATH는 설치과정에서 되니 (/etc/path.d 에 해놓는다) 할 필요가 없지만, blast가 db 파일을 찾을 디렉토리와 환경변수 ($BLASTDB)는 따로 만들어 두는 것이 좋다.

/usr/local/ncbi/blast에서
su mkdir db


로 db 디렉토리를 만들고, ~/.profile 을 편집한다.

vim ~/.profile

난 vi 못씀! 하면 다른 걸 써라. TextEdit이든 Sublime Text이건..
로 .profile에 들어가서
BLASTDB=/usr/local/ncbi/blast/db
export BLASTDB

를 추가.

source ~/.profile

을 하거나 쉘을 재시동.

그 다음에 db를 설치한다. 만들어진 db를 아까 만든 db 디렉토리에 넣어도 되지만 FASTA 형식으로 된 시퀀스를 가지고 db를 설치할 수도 있다. 이걸 해보자.

cd $BLASTDB
pwd
/usr/local/ncbi/blast/db


만들어놓은 DB를 다운받는것도 좋지만 걍 FASTA 형식의 시퀀스를 받아서 db를 만들어보자.

sudo wget ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/pdbaa.gz

해서 다운받고

sudo gunzip pdbaa.gz

압축풀고, 이제 blast용 db를 만들어보자. 이전에는 formatdb라는 것을 썼지만 요즘은 makeblastdb 라는 것을 써야한다고 한다.


makeblastdb -help
USAGE
makeblastdb [-h] [-help] [-in input_file] [-input_type type]
-dbtype molecule_type [-title database_title] [-parse_seqids]
[-hash_index] [-mask_data mask_data_files] [-mask_id mask_algo_ids]
[-mask_desc mask_algo_descriptions] [-gi_mask]
[-gi_mask_name gi_based_mask_names] [-out database_name]
[-max_file_sz number_of_bytes] [-logfile File_Name] [-taxid TaxID]
[-taxid_map TaxIDMapFile] [-version]

….

긴데, 알아야 될 옵션은 두가지뿐.
-in 에서는 blastdb를 만들 입력파일 (여기서는 pdbaa)
-dbtype 에서는 시퀀스의 종류 (핵산의 경우에는 ‘nucl’, 단백질의 경우에는 ‘prot’)


sudo makeblastdb -in pdbaa -dbtype prot

명령을 내리면..

Building a new DB, current time: 02/22/2015 01:26:10
New DB name: pdbaa
New DB title: pdbaa
Sequence type: Protein
Deleted existing BLAST database with identical name.
Keep Linkouts: T
Keep MBits: T
Maximum file size: 1000000000B
Adding sequences from FASTA; added 75499 sequences in 3.41043 seconds.

pdbaa.pin, pdbaa.psq 라는 파일이 생기면 완료.

4. BLAST

이제 BLAST search를 날려봅시다. 검색할 대상이 있는 시퀀스가 저장된 디렉토리로 가서, 가령 Arp2.fasta 식으로 뭐 이런 시퀀스가 텍스트로 저장되어 있는 디렉토리로 이동했다고 치고, 여기서


blastp -query Arp2.fasta -db pdbaa

와 같이 query 뒤에 search하려는 파일, -db 뒤에 검색하려는 db (위의 예에서는 아까 만든 db인 pdbaa) 를 부르면


BLASTP 2.2.30+

Reference: Stephen F. Altschul, Thomas L. Madden, Alejandro A.
Schaffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J.
Lipman (1997), “Gapped BLAST and PSI-BLAST: a new generation of
protein database search programs”, Nucleic Acids Res. 25:3389-3402
….

식으로 화면에 출력.
저장하려면


blastp -query Arp2.fasta -db pdbaa -out outfile
or
blastp -query Arp2.fasta -db pdbaa > out


포맷을 바꾸어서 hit 의 목록만 표시하려면

blastp -query Arp2.fasta -db pdbaa -out outfile -outfmt 7

하면

여튼, 구체적인 옵션은 blastp -help 하면 나옴. BLAST Parsing은 여기를 참조. 지금 21세기도 15년이 다 되었는데 BLAST Parser 따위는 새로 짜지마라. 제발…

Radial Profile Plugin을 이용한 원형 물체의 Profile 측정 (ImageJ)

잊을까봐 적어놓는 ImageJ PlugIn 사용법.

요즘은 요런 그림을 자주 들여다보는데 (쥐 난자의 형광염색사진임. 적색은 actin, 청색은 DNA, 초록색은 ‘어떤’ 단백질) 이런 것을 보다보면 intensity의 profile이 궁금할때가 있음. 그럴때는 대충 이렇게 선을 선택해서

하면

선을 따라서 채널별로 플롯을 그릴수도 있고, 이것을 데이터로 익스포트할수도 있다.

그러나 위의 녹색 채널과 같이 잡시그널이 많은 경우에는

요렇게 노이즈가 많이 나옴. 이걸 해결할 수 있는 방법은 없을까?

가령 요런 식으로 여러개의 선을 그려서 측정한 다음 평균을 내면 어떨까? 그러나 손으로 이걸 하려면 오차도 많이 나고 좀 무리데스.

이런 경우에 쓰라고 어떤 양반이 ImageJ Plugin을 만들어 놓은 게 있다.

Radial Profile Extended 

위의 링크에서 Radial_Profile_Angle_Ext.jar 를 받고 ImageJ의 PlugIn 폴더에 넣는다. (맥인데요, 아이콘만 있고 폴더는 어디있어염? 하는 분이라면 이대로 한다)

여튼 해당 플러그인을 설치한 후 PlugIn메뉴에서 “Radial Profile Angle”을 선택

요런 메뉴가 뜨고, 데이터 화면에는 아래와 같이 원이 그려지는데, 저 원을 주변으로 profile이 그려진다. 화살표를 이용하여 원호의 중심을 조절하거나 Up/down으로 원의 반경을 조절할 수도 있고, 메뉴에서 적절히 조절도 가능.

그 다음에 “Calcurate ROI Profile”을 선택하면

전 영역에서 평균된 intensity profile이 표시됨.

원이 완벽하지 않은 경우에는 아래와 같이 일부 각도에 대해서만 profile을 재는 것이 좀 더 샤프한 profile을 얻을 수 있는 방법.

이전보다 샤프하지?

녹색 채널

두개의 데이터를 export하여 플로팅을 하면

ImageJ에서 stack을 합치는 방법

XDS.INP for Pohang Accelator Laboratory 7A

최근에 포항방사광가속기 beamline 7A에서 데이터 컬렉션을 해왔음. 포항에서 데이터 컬렉션을 한것은 처음인데, 비록 규모는 APS나 Spring-8과 같은 공룡빔라인에 비해서 작지만 비교적 잘 관리되고 업데이트된 곳이라고 생각함.

근데 대개의 빔라인이 그렇듯이 Data Reduction은 HKL2000를 사용하고 있고 HKL2000용 Site.def는 제공되고 있음. 그렇다면 열라비싼 상용프로그램인 HKL2000을 사용하지 않는 사람은 사이트에서 Data Reduction 해가야 되겠네? 이런 경우에는 무료로 제공되는 XDS를 이용하면 편함. 그러나 이건 GUI가 그닥 제대로 되어있지 않지? ㅋㅋㅋ 이런 것을 위해서 한번 XDS.INP를 만들어 보았음.

참고로 여기서 사용하는 Detector는 ADSC의 Q270. 기존의 ADSC용 XDS.INP를 조금 고치면 됨.

DETECTOR=ADSC  MINIMUM_VALID_PIXEL_VALUE=1  OVERLOAD= 65000
DIRECTION_OF_DETECTOR_X-AXIS= 1.0 0.0 0.0
DIRECTION_OF_DETECTOR_Y-AXIS= 0.0 1.0 0.0
TRUSTED_REGION=0.0 1.05 !Relative radii limiting trusted detector region

MAXIMUM_NUMBER_OF_JOBS=2  !Speeds-up COLSPOT & INTEGRATE on a Linux-cluster
MAXIMUM_NUMBER_OF_PROCESSORS=2!0)

X-RAY_WAVELENGTH=0.979          !Angstroem
INCIDENT_BEAM_DIRECTION=0.0 0.0 1.0
FRACTION_OF_POLARIZATION=0.90 !default=0.5 for unpolarized beam;0.90 at DESY;
POLARIZATION_PLANE_NORMAL= 0.0 1.0 0.0
AIR=0.001                     !Air absorption coefficient of x-rays
SPACE_GROUP_NUMBER=0  !0 for unknown crystals; cell constants are ignored.
UNIT_CELL_CONSTANTS= 59.57  69.69  56.97        90  90 90

!FRIEDEL'S_LAW=FALSE !Default is TRUE.
!STARTING_ANGLE=  0.0      STARTING_FRAME=1
!used to define the angular origin about the rotation axis.
!Default:  STARTING_ANGLE=  0 at STARTING_FRAME=first data image
!RESOLUTION_SHELLS=10 6 5 4 3 2 1.5 1.3 1.2
!STARTING_ANGLES_OF_SPINDLE_ROTATION= 0 180 10

!TOTAL_SPINDLE_ROTATION_RANGES=30.0 120 15
!REFERENCE_DATA_SET= CK.HKL   !Name of a reference data set (optional)

INCLUDE_RESOLUTION_RANGE= 40.0 2.0
!==================== SELECTION OF DATA IMAGES ==============================
!Generic file name, access, and format of data images
 NAME_TEMPLATE_OF_DATA_FRAMES=../crystal1_1_1_???.img  SMV DIRECT
 DATA_RANGE=1  360      !Numbers of first and last data image collected
 BACKGROUND_RANGE=1 6  !Numbers of first and last data image for background

!====================== INDEXING PARAMETERS =================================
!Never forget to check this, since the default 0 0 0 is almost always correct!
!INDEX_ORIGIN= 0 0 0          ! used by "IDXREF" to add an index offset

!Additional parameters for fine tuning that rarely need to be changed
!INDEX_ERROR=0.05 INDEX_MAGNITUDE=8 INDEX_QUALITY=0.8
!SEPMIN=6.0 CLUSTER_RADIUS=3
!MAXIMUM_ERROR_OF_SPOT_POSITION=3.0


!================== CRITERIA FOR ACCEPTING REFLECTIONS ======================
 VALUE_RANGE_FOR_TRUSTED_DETECTOR_PIXELS= 6000 30000 !Used by DEFPIX
                   !for excluding shaded parts of the detector.
!INCLUDE_RESOLUTION_RANGE=20.0 0.0 !Angstroem; used by DEFPIX,INTEGRATE,CORRECT

!used by CORRECT to exclude ice-reflections
!EXCLUDE_RESOLUTION_RANGE= 3.93 3.87 !ice-ring at 3.897 Angstrom
!EXCLUDE_RESOLUTION_RANGE= 3.70 3.64 !ice-ring at 3.669 Angstrom
!EXCLUDE_RESOLUTION_RANGE= 3.47 3.41 !ice-ring at 3.441 Angstrom
!EXCLUDE_RESOLUTION_RANGE= 2.70 2.64 !ice-ring at 2.671 Angstrom
!EXCLUDE_RESOLUTION_RANGE= 2.28 2.22 !ice-ring at 2.249 Angstrom

!WFAC1=1.0  !This controls the number of rejected MISFITS in CORRECT;
        !a larger value leads to fewer rejections.


!============== INTEGRATION AND PEAK PROFILE PARAMETERS =====================
!Specification of the peak profile parameters below overrides the automatic
!determination from the data images
!Suggested values are listed near the end of INTEGRATE.LP
!BEAM_DIVERGENCE=  0.473  !arctan(spot diameter/DETECTOR_DISTANCE)
!BEAM_DIVERGENCE_E.S.D.=   0.047 !half-width (Sigma) of BEAM_DIVERGENCE
!REFLECTING_RANGE=  1.100 !for crossing the Ewald sphere on shortest route
!REFLECTING_RANGE_E.S.D.=  0.169 !half-width (mosaicity) of REFLECTING_RANGE

!NUMBER_OF_PROFILE_GRID_POINTS_ALONG_ALPHA/BETA=9 !used by: INTEGRATE
!NUMBER_OF_PROFILE_GRID_POINTS_ALONG_GAMMA= 9     !used by: INTEGRATE

!CUT=2.0    !defines the integration region for profile fitting
!MINPK=75.0 !minimum required percentage of observed reflection intensity
!DELPHI= 5.0!controls the number of reference profiles and scaling factors

!PATCH_SHUTTER_PROBLEM=TRUE         !FALSE is default
!STRICT_ABSORPTION_CORRECTION=FALSE !TRUE  is default


!=========== PARAMETERS DEFINING BACKGROUND AND PEAK PIXELS =================
!STRONG_PIXEL=3.0                              !used by: COLSPOT
!A 'strong' pixel to be included in a spot must exceed the background
!by more than the given multiple of standard deviations.

!MAXIMUM_NUMBER_OF_STRONG_PIXELS=1500000       !used by: COLSPOT

!SPOT_MAXIMUM-CENTROID=3.0                     !used by: COLSPOT

!MINIMUM_NUMBER_OF_PIXELS_IN_A_SPOT=6          !used by: COLSPOT
!This allows to suppress spurious isolated pixels from entering the
!spot list generated by "COLSPOT".

!NBX=3  NBY=3  !Define a rectangle of size (2*NBX+1)*(2*NBY+1)
!The variation of counts within the rectangle centered at each image pixel
!is used for distinguishing between background and spot pixels.

!BACKGROUND_PIXEL=6.0                          !used by: COLSPOT,INTEGRATE
!An image pixel does not belong to the background region if the local
!pixel variation exceeds the expected variation by the given number of
!standard deviations.

!SIGNAL_PIXEL=3.0                              !used by: INTEGRATE
!A pixel above the threshold contributes to the spot centroid


!================= PARAMETERS CONTROLLING REFINEMENTS =======================
!REFINE(IDXREF)=BEAM AXIS ORIENTATION CELL !DISTANCE
!REFINE(INTEGRATE)=!DISTANCE BEAM ORIENTATION CELL !AXIS
!REFINE(CORRECT)=DISTANCE BEAM ORIENTATION CELL AXIS

사실 여기서 제일 중요한 게 다른것보다는 ORGX=2133.0 ORGY=2080.0 !Detector origin (pixels). ORGX=NX/2; ORGY=NY/2
이걸 정하는건데 대충 NX, NY (이미지의 X Width, Y Height) 의 절반이라고 나와있지만 사실 빔센터가 항상 정확한게 아님. ㅋ

제일 좋은 방법은 빔센터를 직접 확인하는 것인데, adxv와 같은 뷰어로 img 파일을 불러보면

요런식으로 빔센터의 X,Y Pixel 좌표를 찍어서 이것을 ORGX, ORGY로 잡으면 적절함.

Simple GUI using Bottle.py and twitter bootstrap

https://github.com/arq5x/toy_bottle_bootstrap_app

스크립트/명령행 기반 커맨드에 간단한 웹 기반 인터페이스를 다는 예제. 사실 비슷한 짓거리를 하고 있었는데 BedTool 을 만든 아저씨가 아주 간단한 예를 만들어 올려놓았다. ㅋ

즉 파일 하나로 해결되는 웹 마이크로프레임웍인 Bottle.py으로 간단한 서버와 웹프레임웍을 구축하고 Twitter Bootstrap 으로 프론트엔드를 만듬. (출력은 Jinja2 를 Template Engine으로 이용)

스크립트 실행하고 웹브라우저에서 localhost:8088에 접속하면 인터페이스 실행 (사실 이것도 스크립트 실행후에 웹브라우저 띄워줄수도 있지만)

자세한 것은 링크 참조.

GIST in wordpress.com

이전까지 소스코드를 블로그에 적을떄 <pre>를 쓰거나 아니면

태그를 썼었는데, 알고보니 wordpress.com에서는 gist 를 지원하고 있었다. ;;;

gist에 코드쪼가리를 만든다음, 워드프레스에서 글 작성하면서 그냥 gist URL 쳐넣으면

	#!/bin/bash
	mysql -u root allbacteria < gff.sql
	mysql -u root allbacteria < accession.sql
	awk '{ if ($3=="CDS") print $1"\t"$4"\t"$5"\t"$7"\t"$9}' $1.gff | sed -e 's/ID=.*Name=//' -e 's/;Parent.*$//' > $1.txt
	cat $1.txt | mysql -u root allbacteria –local-infile=1 -e "LOAD DATA LOCAL INFILE '/dev/stdin' INTO TABLE gff;"
	awk '/^>/ {split($0,a,"|");print a[4],"\t",$1}' $1.fasta | sed 's/>//' > $1.txt2
	cat $1.txt2 | mysql -u root allbacteria –local-infile=1 -e "LOAD DATA LOCAL INFILE '/dev/stdin' INTO TABLE accession;"

view raw

gistfile1.sh

hosted with ❤ by GitHub

gist 코드가 요기 임베딩되어있넹? 뭔코드인지는 묻지마셈 걍 흔한 awk 파싱해서 MySQL 쳐넣처넣

Multifasta에서 빠르게 일부 레코드만 추출해오기

이런 작업은 사실 일상적으로 하는 작업인 관계로 Perl이나 Python 과 같은 스크립트 언어를 조금이라도 아는 사람이라면 그닥 어렵지 않게 할 수 있는 일이긴 하지만, 이런 작업도 수GB 이상의 파일이라면 그리고 불러올 파일의 갯수가 늘어나면 적지 않은 시간이 소요된다. 좀 더 빠르게 이런 작업을 마칠 수 없을까?

Profile HMM 에 의해 도메인 및 단백질 검색을 하는 유명한 소프트웨어인 HMMER 3.0 에 딸려있는 유틸리티인 esl-sfetch를 이용하면 아주 빠르게 이런 작업을 할 수 있다.

HMMER3 을 설치하고 hmmer3 의 실행파일에 PATH를 걸어둔 다음

esl-sfetch -h

하면

$ esl-sfetch -h
# esl-sfetch :: retrieve sequence(s) from a file
# Easel h3.0 (March 2010)
# Copyright (C) 2010 Howard Hughes Medical Institute.
# Freely distributed under the Janelia Farm Software License.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: esl-sfetch [options]          (one seq named )
Usage: esl-sfetch [options] -f   (all seqs in )
Usage: esl-sfetch [options] --index        (index )

 where general options are:
  -h     : help; show brief info on version and usage
  ....중략.

가령 지금 나는 vibrio.fasta라는 274MB짜리 fasta 파일이 있고, 이 파일은

$tail vibrio.fasta
SDGATWLDRRLIHGETADLAPTGFGQQVREAMDQRREHHIEQGDATRSRDSRVFYRRNLLAILREREVAG
VGSDMALSKGLPFRAATDGESVSGKFTGTVHLSSGKFAVVEKSHEFTLVPWRPIIDRQLGREVMGIVQGG
SVSWQLGRQRGLER
>gi|424025505|ref|ZP_17765144.1| drug resistance transporter, Bcr/CflA subfamily protein [Vibrio cholerae HC-62B1]
MTTTRPAWAYTLPAALLLMAPFDILASLAMDIYLPVVPAMPGILNTTPAMIQLTLSLYMVMLGVGQVIFG
PLSDRIGRRPILLAGATAFVIASLGAAWSSTAPAFVAFRLLQAVGASAMLVATFATVRDVYANRPEGVVI
YGLFSSMLAFVPALGPIAGALIGEFLGWQAIFITLAILAMLALLNAGFRWHETRPLDQVKTRRSVLPIFA
SPAFWVYTVGFSAGMGTFFVFFSTAPRVLIGQAEYSEIGFSFAFATVALVMIVTTRFAKSFVARWGIAGC
VARGMALLVCGAVLLGIGELYGSPSFLTFILPMWVVAVGIVFTVSVTANGALAEFDDIAGSAVAFYFCVQ
SLIVSIVGTLAVALLNGDTAWPVICYATAMAVLVSLGLVLLRLRGAATEKSPVV

요런 식으로 되어 있는 multifasta 파일이다. 여기서 gi|424025505|ref|ZP_17765144.1| 와 같은 accession 으로 파일을 찾고싶다고 치자.

gi|260362785|ref|ZP_05775654.1|
gi|260880072|ref|ZP_05892427.1|
gi|28901226|ref|NP_800881.1|
gi|153835952|ref|ZP_01988619.1|
gi|308095182|ref|ZP_07663192.1|
gi|254224936|ref|ZP_04918551.1|
gi|229529188|ref|ZP_04418578.1|
gi|254285269|ref|ZP_04960234.1|
gi|421351525|ref|ZP_15801890.1|
gi|229515173|ref|ZP_04404633.1|

요런 내용을 list.txt로 저장했다고 치자
제일 먼저 할 일은 fasta 파일을 인덱싱하여 빨리 검색을 할 수 있게 하는거.

$esl-sfetch --index vibrio.fasta
Creating SSI index for vibrio.fasta...    done.
Indexed 665311 sequences (665311 names).
SSI index written to file vibrio.fasta.ssi

이러면 파일이름.ssi 식의 인덱스 파일이 생기고 이제 검색을 할 수 있다.

$esl-sfetch -f vibrio.fasta list.txt

하면 바로 쫙 stdout 으로 파일이 추출된다. ㅋ

사실 이런 작업은 hmmer 등을 수행해서 나온 결과파일을 가지고 해당하는 fasta 파일을 데이터베이스에서 추출할때 유용한데, hmmer를 수행해서 테이블 형식으로 다음과 같은 출력파일이 나왔다고 하자.

이 출력물을 awk 등으로 파싱해서 원하는 gi 값만 esl-sfetch에 전달해 주면 곧바로 hmmer hit 에 상응하는 시퀀스 파일을 뽑을 수 있다.

awk '{if (substr($0,0,1)!="#") {print $1}}' VopL.fasta.hmmer | esl-sfetch -f vibrio.fasta - > hmmerresults.fasta

참고로 samtools 에도 faidx 라는 동일한 기능 (fasta 를 indexing한후 추출) 이 있으나 테스트해 본 결과 esl-sftech가 더 빠르다.
약 247MB의 fasta에서 10개 정도의 fasta를 추출할때 시간을 재보면

time awk '{if (substr($0,0,1)!="#") {print $1}}' VopL.fasta.hmmer | xargs samtools faidx vibrio.fasta
real	0m1.033s
user	0m0.921s
sys	0m0.110s

time awk '{if (substr($0,0,1)!="#") {print $1}}' VopL.fasta.hmmer | esl-sfetch -f vibrio.fasta
real	0m0.007s
user	0m0.005s
sys	0m0.006s

수십 GB 정도의 파일이라면 좀 더 차이날 것 같다.;;

-C 옵션을 이용하면 multifasta 내의 subsequence을 잘라오는 것도 가능하다. (BED file 등을 파싱해서 빠르게 시퀀스를 얻는것도 가능..아 그건 bedtools쓰면 되나..여튼.

어쩄든 NGS 시대가 되면서 처리해야 할 데이터의 양이 부쩍 늘어난 상황에서 파싱과 같은 간단한 작업의 경우에도 가능하면 unix 기반의 텍스트 툴및 파이프, 리다이렉션을 이용하는 것이 수행속도 면에서도 훨씬 이점이 있는 것 같다는.

1+1은..

요즘 파이썬으로 작성한 스크립트에 최소한의 UI를 씌워보고자 몇가지를 살펴보고 있는데, 요즘 유행하는 Microframework 을 이용하여 간단히 웹 기반의 인터페이스를 씌우는 게 가장 편한 일인듯.

Microframework 의 경우에도 여러가지가 있겠지만, 단 한개의 파일로 해결되는 bottle 을 써볼 예정. 테스트용으로 만든 코드..;;

#!/usr/bin/python
# -*- coding: utf-8 -*-
from bottle import route, run, template
import webbrowser

@route('/:name')
def index(name):
    if int(name) < 6:

        list = u"{0}+{0}은".format(name)
        i = int(name)
        for k in range(i):
            list = list + u"귀요";
        list = list + u"미";
    else:
        list = u"그만해;;;;지겹다";
    return template('<b>{{name}}</b>!', name=list)
webbrowser.open('http://localhost:8080/1')
run(host='localhost', port=8080)

이 파일을 run.py 뭐 이런식으로 저장해서 디렉토리에 넣고, bottle.py 를 다운로드받아 같은 디렉토리에 넣고, 터미널에서

./run.py

로 실행하면…

…..;;;;;;;

http://localhost:8080/2 로 주소를 바꿔보면..

6+6은 뭘까? 직접 알아보도록 하고 (쿨럭) 여튼 이렇게 하면 간단히 파일 하나로 웹어플리케이션 구축가능. ㅋ 스크립트 형태로 있는 것에 웹 UI를 달아야 할때와 같을때 매우 유용할 듯 싶다. 특히 twitter bootstrap와 같은 템플릿과 같이 이용하면 최소한의 노력으로 그럴싸한 외관의 웹 UI를 씌울 수 있게 될듯.

여러개의 PDB 파일을 한꺼번에 렌더링하기 : PyMOL Batch Renderer

단백질 구조를 가지고 놀다보면 디렉토리에 PDB 파일이 쌓이게 되고..가끔은 이 파일이 무슨 구조의 파일인지조차 헷갈리는 경우가 있다.

물론 일일히 PyMOL과 같은 프로그램으로 열어보든지 아니면 PDB 사이트에서 검색을 해보셈…할수도 있겠으나 한두개면 모르겠는데 갯수가 늘어나면 난감하기 마련.

“디렉토리에 들어있는 PDB 파일을 한번에 몽땅 그림으로 보여주는 방법이 없을까?”

그런 방법 별로 생각이 안나네염.ㅋㅋ 그러나 이제 생겼음.

#!/usr/bin/python
#
# Batch Pymol Renderer
# Require Python 2.7 (argparse) and Pymol
#
#

import sys, subprocess, re, os, glob
import argparse

def pdbparsing(filename, regex):
#
# Parse pdb content and extract using regex
#
	if os.path.exists(filename):
		pdb = open(filename)
		pdbcontent = pdb.readlines()
		matched = ''
		matchedcount=0
		reg = re.compile(regex)
		for line in pdbcontent:
			match = reg.match(line)
			if match:
				if matchedcount==0:
					gap = " "
				else:
					gap = ""
				matched = matched+gap+match.group(1).strip()
		return (matched)
	else:
		print "{0} is not found!".format(filename)
		sys.exit()


def htmlout(description):
#
# Generate render.html which display rendered png files
#
#
	htmlheader = """
<!DOCTYPE html>
<html>
	<head>
    	<title>List</title>
 		<meta charset='utf-8'>
	</head>
	<style>
	</style>
	<body>
"""
	imgtag = """
		<img src="./{0}">
		<p><a href="http://www.rcsb.org/pdb/explore/explore.do?structureId={1}">{2}</a></p>
"""

	htmlfooter = """
	</body>
</html>
"""
	htmloutput = open('render.html','w')
	htmloutput.write(htmlheader)
	for pdb,desc in description.items():
		htmloutput.write(imgtag.format(pdb[:-4]+".png", pdb[:len(pdb)-4], pdb[:len(pdb)-4]+":"+desc))
	htmloutput.write(htmlfooter)
	htmloutput.close()


def main(arglist):
#
# PyMol executable path
# In cases of MacPymol, It assumed that pymol is installed in 'Application' folder.
# Other cases, plase set up the exact location of PyMOL executable
#
	PyMOLPath = '/Applications/MacPyMOL.app/Contents/MacOS/MacPyMol'
	if len(arglist.file)>0:
		subjectList = arglist.file
	else:
		subjectList = glob.glob('*.pdb')
	if len(subjectList)>0 :
		pymolLoadingFile = open('render.pml','w')
		pymolLoadingFile.write('bg_color white\n')
		description = {}

		style = "cartoon" #Defaults
		if arglist.ribbon :
			style = "ribbon"
		if arglist.line :
			style = "line"
		if arglist.stick :
			style = "stick"

		if arglist.view:
			if os.path.exists(arglist.view):
				f = open(arglist.view)
				view = f.readlines()
				if re.match("^set_view",view[0]):
					for line in view:
						pymolLoadingFile.write(line)

		for pdb in subjectList:	
			if os.path.exists(pdb):			
				print pdb
				pymolLoadingFile.write("load {0}/{1}\n".format(os.getcwd(),pdb))
				pymolLoadingFile.write("hide all\n")
				pymolLoadingFile.write("show {0}, {1}\n".format(style, pdb[:-4]))
				if arglist.cbc:
					pymolLoadingFile.write("util.cbc\n")
				else:
					pymolLoadingFile.write("spectrum count,rainbow,{0}\n".format(pdb[:-4]))
				if arglist.surface:
					pymolLoadingFile.write("show surface,{0}\n".format(pdb[:-4]))
				if not arglist.view:	
					pymolLoadingFile.write("orient {0}\n".format(pdb[:-4]))
				if arglist.ray:
					pymolLoadingFile.write("ray\n")
				pymolLoadingFile.write("png {0}.png\n".format(pdb[:-4]))
				description[pdb]=pdbparsing(pdb, '^TITLE    [ \d](.*)$')
		pymolLoadingFile.close()
		print "Running PyMol..."
		p = subprocess.Popen([PyMOLPath,'-c','render.pml'],
							  stdout=subprocess.PIPE)
		p_stdout = p.stdout.read()
		htmlout(description)
	else:
		print "Usage : render.py"
		sys.exit()		


if __name__ == "__main__":
	parser = argparse.ArgumentParser()
	parser.add_argument('-r', '--ribbon', action='store_true', dest='ribbon',default=False,
	                    help='Draw as ribbon')
	parser.add_argument('-l', '--line', action='store_true', dest='line',default=False,
	                    help='Draw as Line')
	parser.add_argument('-t', '--stick', action='store_true', dest='stick',default=False,
	                    help='Draw as Stick')
	parser.add_argument('-s', '--surface', action='store_true', dest='surface',default=False,
	                    help='Draw as Surface')
	parser.add_argument('-c', '--color_by_chain', action='store_true', dest='cbc',default=False,
	                    help='Color by Chain')
	parser.add_argument('-rt', '--ray_trace', action='store_true', dest='ray',default=False,
	                    help='Ray tracing')
	parser.add_argument('-f', '--files', nargs='*', dest='file',default=[],
	                    help='File to process')
	parser.add_argument('-v', '--fixed_view', action='store', dest='view')
	results = parser.parse_args()
	main(results)

위 스크립트를 render.py 등으로 저장하고 실행권한을 준다음 PATH가 설정되어 있는 디렉토리에 넣든 pdb 가 디글거리는 파일에 넣든지 하시고,

Python 2.7에서 지원되는 argparse 모듈을 사용하고 있으므로 Python 2.7 이상을 사용하든지 argparse를 설치. 
그리고 PyMol 이 설치되어 있어야 하고, MacPyMol 이 아닌 경우에는 75라인의 
PyMOLPath = '/Applications/MacPyMOL.app/Contents/MacOS/MacPyMol' 
을 자신의 PyMol 설치경로로 바꾸어 주어야 한다 

실행옵션을 보면

./render.py -h
usage: render.py [-h] [-r] [-l] [-t] [-s] [-c] [-rt]

optional arguments:
  -h, --help            show this help message and exit
  -r, --ribbon          Draw as ribbon
  -l, --line            Draw as Line
  -t, --stick           Draw as Stick
  -s, --surface         Draw as Surface
  -c, --color_by_chain  Color by Chain
  -rt, --ray_trace      Ray tracing

일단 아무 옵션 없이 render.py 를 실행시켜보면 render.pml 이라는 파일이 생성되고 자동적으로 스크립트내에서 PyMOL 을 실행시킨후 디렉토리내에 있는 모든 PDB를 하나씩 부른후 이를 png 로 저장한다.

그게 다가 아니라 render.html 이라는 html 형식의 카탈로그를 만들어줌. 설명을 클릭하면 rcsb pdb의 해당 pdb 페이지로 이동.

링크

여러가지 옵션을 제공하는데 (디폴트는 cartoon 방식의 렌더링이고, 색 지정은 N말단 – 청색 C말단 – 적색의 스펙트럼)

./render.py -r

Ribbon 으로 렌더링

./render.py -t
stick 으로 렌더링

./render.py -t
line으로 렌더링

./render.py -s
surface로 렌더링 (다른 표현방식에 비해 장시간이 걸리는 것에 유의)

./render.py -c
color by chain (chain별로 다른 color지정하여 렌더링)방식으로 색지정 변환.

./render.py -c -s
surface로 rendering하면서 color by chain

./render.py -rt
Raytracing (고품질의 렌더링이 가능하나 시간이 매우 오래 걸림에 유의)

추가로 부연설명하면

1. PyMol을 Python Script에서 제어하는 방법으로 가장 속편한 방법이라면 여기서 쓰는 방법처럼 PyMol script (위의 예에서는 Render.pml) 를 만들고 PyMol을 수행하여 작업을 수행하는 것임. 이렇게 생성된 스크립트를 대충 살펴보면

bg_color white
#바탕색 하얀색
load /Users/suknamgoong/Dropbox/pdbsearch/S_1G9O.pdb
#로딩. 
hide all
#몽땅지우고
show cartoon, S_1G9O
#S_1G9O이라는 pdb를 cartoon으로 그림
util.cbc
#color by chain
orient S_1G9O
#최적의 위치로 카메라위치 옮김
png S_1G9O.png
#S_1G90.png라는 이름으로 그림 png로 저장
.....반복

2. PyMol 을 커맨드 라인에서 수행하려면

PyMol -c script.pml

와 같은 방식으로 수행하면 편함.

BLAST를 이용한 PDB 자동 다운로드 및 처리

PDB에서 유사한 구조, 혹은 도메인을 모두 검색해서 구조를 서로 비교해 본다면 어떻게 해야 할까? 1. PDB에 가서 2. 키워드로 검색하고 3. 나오는 PDB 파일을 하나하나 클릭해서 저장하고 4. PyMol 등에 불러와서 5. Structure alignment 를 하고..
물론 몇 개 정도의 단백질이라면 이렇게 해도 된다. 그렇지만 비슷한 구조를 가진 단백질을 다 비교해 보고 싶다면? 가령 Kinase 로 구조가 풀려있는 PDB를 모두 검색해서 이들의 구조를 다 비교해 보고 싶다든지 한다면…이거 수작업으로는 쉽지 않다.

특정 아미노산 서열을 가지고 PDB에서 BLAST 검색해서 PDB 다운로드하고, 이걸 PyMol에 로딩한후 제일 처음 파일 기준으로 Structural Alignment를 해주는 스크립트. 자세한 설명은 나중에 시간나면…안나면말구
단, Pymol 버전이 1.3이상이거나 cealign (http://pymolwiki.org/index.php/Cealign) 이 설치되어 있어야 Structural Alignment 가 됨.

#
#
# 1. Blast search of PDB using RESTful web service form www.rcsb.org
# 2. Retrieve PDB ids and PDB files
# 3. Generate PyMOL script which load all of PDBs and structural alignment using cealign (http://pymolwiki.org/index.php/Cealign)
# Written by MadScientist (https://madscientist.wordpress.com
#
import urllib2, urllib, os
import xml.etree.ElementTree as ET
#
# The first Query is BLAST search using sequence and eValue Cutoff
# and second filter is high resolution Cutoff
# the last filter is for excluding title contains "arp" (To exclude Arp 2/3 complex)
# 
url = 'http://www.rcsb.org/pdb/rest/search'
queryXML = """
<?xml version="1.0" encoding="UTF-8"?>
<orgPdbCompositeQuery version=\"1.0\">
	<queryRefinement>
		<queryRefinementLevel>0</queryRefinementLevel>
			<orgPdbQuery>
				<queryType>org.pdb.query.simple.SequenceQuery</queryType>
				<sequence>{0}</sequence>
				<eCutOff>{1}</eCutOff>
				<searchTool>blast</searchTool>
			</orgPdbQuery>
</queryRefinement>
<queryRefinement>
	<queryRefinementLevel>1</queryRefinementLevel>
		<orgPdbQuery>
			<queryType>org.pdb.query.simple.ResolutionQuery</queryType>
			<description>ResolutionQuery </description>
			<refine.ls_d_res_high.comparator>between</refine.ls_d_res_high.comparator>
			<refine.ls_d_res_high.max>{2}</refine.ls_d_res_high.max>
		</orgPdbQuery>
</queryRefinement>
<queryRefinement>
	<queryRefinementLevel>2</queryRefinementLevel>
		<orgPdbQuery>
			<queryType>org.pdb.query.simple.StructTitleQuery</queryType>
			<struct.title.comparator>!contains</struct.title.comparator>
			<struct.title.value>arp</struct.title.value>
		</orgPdbQuery>
</queryRefinement>
</orgPdbCompositeQuery>
"""
#
# Guess what is this sequence 🙂 
#
#
sequence = """
MDDDIAALVVDNGSGMCKAGFAGDDAPRAVFPSIVGRPRHQGVMVGMGQKDSYVGDEAQS
KRGILTLKYPIEHGIVTNWDDMEKIWHHTFYNELRVAPEEHPVLLTEAPLNPKANREKMT
QIMFETFNTPAMYVAIQAVLSLYASGRTTGIVMDSGDGVTHTVPIYEGYALPHAILRLDL
AGRDLTDYLMKILTERGYSFTTTAEREIVRDIKEKLCYVALDFEQEMATAASSSSLEKSY
ELPDGQVITIGNERFRCPEALFQPSFLGMESCGIHETTFNSIMKCDVDIRKDLYANTVLS
GGTTMYPGIADRMQKEITALAPSTMKIKIIAPPERKYSVWIGGSILASLSTFQQMWISKQ
EYDESGPSIVHRKCF
"""
resolution="2.8"
#
# if you want to search remote homolog, increase E value cutoff. In this case, this protein is very highly conserved protein, so I used very stringent threshold.
#
eCutOff="1e-30"
biological_unit = True
reprot = """
http://www.rcsb.org/pdb/rest/customReport?pdbids=0&customReportColumns=structureId,structureTitle,resolution,experimentalTechnique,depositionDate,structureAuthor,classification,structureMolecularWeight&service=wsdisplay&format=xml
"""

queryText = queryXML.format(sequence, eCutOff, resolution)
print "querying PDB...\n"

#
# Fetch pdb list using XML query in queryText
# and pdb list was in pdblist
#

req = urllib2.Request(url, data=queryText)
f = urllib2.urlopen(req)

result = f.read()

if result:	
	pdblist = []
	pdblist = result.split()


	if biological_unit:
		pdburl = "ftp://ftp.wwpdb.org/pub/pdb/data/biounit/coordinates/all/{0}.pdb1.gz"
	else:
		pdburl = "http://www.rcsb.org/pdb/files/{0}.pdb.gz "
	
	f=open('load.pml','w')
	
	for pdb in pdblist:
		print "download {0}".format(pdb)
		urllib.urlretrieve(pdburl.format(pdb.lower()), pdb+".pdb.gz")
		os.system("gunzip {0}.pdb.gz".format(pdb))
		f.write("load {0}/{1}.pdb\n".format(os.getcwd(),pdb))
		
	f.write("hide all\n".format(pdb))
	f.write("show cartoon\n".format(pdb))

	for i in range(1,len(pdblist)):
		f.write("cealign {0}, {1}\n".format(pdblist[0],pdblist[i]))
	f.close()
	
	req = urllib2.Request(reprot.format(",".join(pdblist)))
	f = urllib2.urlopen(req)
	result = f.read()

	if result:
		f=open('list.txt','w')	
		root = ET.fromstring(result)
		for record in root.iter('record'):
			recordline = []
			for field in record:
				recordline.append(field.text)
			f.write("\t".join(recordline)+"\n")
		f.close()


	else:
	    print "Failed to retrieve results" 
			
else:
    print "Failed to retrieve results" 

어쨌든 실행하면

querying PDB...
download 1ATN
download 1C0F
download 1C0G
download 1D4X
download 1DEJ
download 1EQY
download 1ESV
download 1HLU
download 1J6Z
download 1KXP
download 1LOT
download 1MA9
download 1MDU
download 1NLV
download 1NM1
download 1NMD
download 1NWK
download 1P8Z
download 1QZ5
download 1QZ6
download 1RDW
download 1S22
download 1SQK
download 1WUA
download 1YAG
download 1YVN
download 1YXQ
download 2A3Z
download 2A40
download 2A41
download 2A42
download 2A5X
download 2ASM
download 2ASO
download 2ASP
download 2BTF
download 2D1K
download 2FF3
download 2FF6
download 2FXU
download 2GWJ
download 2GWK
download 2HF3
download 2HF4
download 2HMP
download 2OAN
download 2PAV
download 2PBD
download 2Q0R
download 2Q0U
download 2Q1N
download 2Q31
download 2Q36
download 2Q97
download 2V51
download 2V52
download 2VYP
download 3A5L
download 3A5M
download 3A5N
download 3A5O
download 3CHW
download 3CI5
download 3CIP
download 3DAW
download 3EKS
download 3EKU
download 3EL2
download 3HBT
download 3M6G
download 3MMV
download 3MN5
download 3MN6
download 3MN7
download 3MN9
download 3SJH
download 3U4L
download 3U8X
download 3U9D
download 3U9Z
download 3UB5
download 3UE5
download 4EFH

BLAST에서 hit 된 pdb 파일을 쳐묵쳐묵 다운로드하고

이렇게 생성된 load.pml 파일을 더블클릭하거나 PyMol에서 @load.pml 식으로 불러오면

(파일이 수십개가 넘어가면 cealign 단계에서 시간 꽤 걸림. 인내심을 가지고 기다려보삼)

결과물. 헉 단백질이 디글디글. 너무 정신없으니까 좀 엇나간 넘들 좀 끄면..

대충 특정한 단백질 (Actin)에 붙는 바인딩 파트너들이 어떠한 방식으로 붙는지를 개괄적인 양상을 볼 수 있음.

list.txt를 엑셀 등으로 열어보면 다음과 같이 다운로드된 파일의 목록을 볼 수 있음.

위 스크립트를 응용하면 좀 더 다양한 조건을 지정하여 PDB의 목록을 뽑고, 이를 다운로드하고, PDB에 대한 관련정보들을 저장하여 추후분석에 사용할 수 있을 것임.

시간이 나면 CCP4의 ‘superpose’ http://www.ccp4.ac.uk/html/superpose.html 를 이용하여 PyMol 의 cealign 을 거치지 않고 바로 Structure superimpose를 하도록 스크립트를 고칠 예정. 시간나면. ;;;;