ncbi 횽아들은 어디까지 만들어 낼 것인가

이것저것 작업하면서

ncbi tool들을 다시 사용하고 있는데..

훗.. 역시 우리 ncbi 훃아들의 위대함을 다시 한번 느꼈다는...

NCBI BLAST에서 taxonomy로 제한 거는 기능을 당연히 stand-alone에서도 사용할 수 있는데 NCBI의 -taxid의 숨은 함정이 종 수준의 taxid만 제한 걸 수 있다는..

(근데 써보면 종 수준의 taxid만 제한이 걸리는지 갸우뚱 거리긴 함.. )

여튼 종 수준의 taxid만 제한할 수 있다는 것이 무엇이냐면..

종보다 상위 class의 taxid인 Enterobacterales의 taxid를 사용하면 정상적으로 작동을 안하게 된다는 말씀.

그러므로 NCBI BLAST 프로그램을 다운 받았을 때 함께 있는 get_species_taxids.sh를 활용하면 이 문제를 피해갈 수 있다고 합니다.

사실 최근까지 get_species_taxid.sh가 왜 있는지 관심은 없을 뿐더러
왜 쓰잘떼기 없는 shell script는 왜 넣어놨는지 했다는 ㅎㅎ 

여튼 언제나 NCBI 훃님들께 감사인사를... :)

참고 URL: https://www.ncbi.nlm.nih.gov/books/NBK569846/

출처: @ye._.vely618

mummer4, 미처 알아보지 못했다

mummer는 서열 정열프로그램으로 꽤나 오래부터 사용되었던...
NGS시대에 접어들면서 일반적으로 사용되지는 않고...
Reference 제작할때 종종 쓰이는...

MUMmer이라고 보통 쓰는데... MUM의 의미가 "Maximal Unique Matches"라고.. 여기에 나와있었네요..

그리고 이번에 알았는데.. MUMmer는 당연히 target와 query는 fasta형식만 입력될줄 알았는데 input으로 fastq도 사용할 수 있다는..
당연히 fasta 형식으로 변경해서 사용하려고 했는데... 작동해서 잠깐 화장실좀 다녀왔다능..

input으로 fastq형식을 받을 수 있는 버전의 tools이 2018년도에 논문으로 출판되어 한번 들추어 보았습니다.

이름하여 MUMmer4: A fast and versatile genome alignment system

이전 MUMmer3 이후 데이터구조를 32bit에서 48bit로 증가시켜 이론적으로 비교가능한 크기가 141Tbp (다시 한번 화장실을 다녀오게 만드는.....)로 늘렸다는.. 근데 이거 141Tbp이 입력 파일이면 입력 파일 로딩하는데만 한세월 아닌가...

여튼..

MUMmer4 논문 작성 할때 언급해준 정렬 프로그램으로
그냥 생명정보학 분석한다고 할때 기본 옵션인 BWA/Bowtie,
PacBio 데이터 할때 사용하는 BLASR 들을 언급해 주셨는데
이제 fastq형식의 파일을 지원해서
reference vs reference 비교 프로그램이 NGS용 정렬 프로그램들에게 어깨를 나란히 할 수 있는 기회를...

mummer라고 하면 갱장한 legacy 프로그램이라고 생각할 수 있지만...
다방면으로 개선 시킨 mummer4를 하나하나 뜯어보면 갱장히 힙해졌다고 볼 수 있습니다.

NGS 시대에 다양하게 쏟아져 나오는 reference크기에 대응할 수 있도록 비교 사이즈의 증가(141Tb)되었고, 또한 긴 서열 작업시 메모리 문제가 발생할 것을 대비해서 긴 서열들을 자동으로 분할해서 작업해주는 --batch라는 옵션기능까지...

mummer4의 새로운 기능들이 갱장히 많아졌습니다.

reference 서열을 제작하는 de-novo assembly하시는 분들 외에도 다양한 작업시 사용할 수 있을것으로 생각되니 한번 기능들 구경해보시고 활용해보시기 바랍니다.

그런데... 사실 NGS 서열 정렬은 역시 여기가 맛집이긴 합니다.

ps. BLAST에서도 query파일을 fastq 형식을 받을 수 있다고 합니다. 이름하여 Magic-BLAST를 사용하시면되겠습니다. 속도는 잘 모르겠습니다. ㅎㅎ 

@ye._.vely618

 

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

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

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

Local에서 Blast 작업 돌리기

BLAST는 최근 생명공학을 하는데 기본 도구중에
하나가 되었지만 그럼에도 많은 연구자들이 BLAST를
제대로 사용하지 못하고 있는 것이 현실이다.

BLAST를 수행하는데 Database가 NCBI나 다른 여타의 사이트에서
제공안되는 Database를 이용해야만 되는 경우가 발생하면 어떻게 할 것인가?

혹은 BLAST를 해야할 Query가 수십개가 아닌 수백, 수천개라면
어떻게 할 것인가?

이런경우 자신의 컴퓨터에서 BLAST를 수행하게 된다면 원하는 작업을
손쉽고 빠르게 할 수 있다.

BLAST 프로그램은 대부분 연구자들이 알고 있듯이
NCBI에서 다운로드 받을 수 있다.
ftp://ftp.ncbi.nih.gov/blast/executables/release
위의 ftp 주소에 들어가 자신의 플랫폼에 맞는 파일을 다운로드 받으면 일단 BLAST를
수행할 수 있는 준비가 된다.

다운로드 받은 압축 파일(실행파일(*.exe)로 압축되어있는)을 풀면
bin, data, doc 세개의 폴더가 나타난다.
BLAST를 직접 수행하는 실행파일은 bin 폴더안에 있다.

기본적으로 BLAST를 사용하기 위해서 두개의 파일이 필요하다.
blastall과 formatdb이다.
blastall은 일반적인 blast, 즉 blastn, blastp, blastx, tblastx,tblastn를 수행할 때 사용된다.
formatdb는 blast를 할 수 있는 database를 만들어 주는 파일이다.
blast에 사용되는 database는 항상 ncbi나 다른 웹사이트에서 제공해주는 것이 아니기
때문에 자신만의 database를 만들 수 있어야 한다.

formatdb -i INPUT_FILE -p T|F -o T|F

- i INPUT_FILE은 fasta form을 따르는 서열들이 모인 파일이면 문제없다.
- p INPUT_FILE이 DNA서열인지  Protein 서열인지 확인하는 옵션값 protein의 경우 T
- o INPUT_FILE를 paser하는 옵션 NCBI에서 수집한 서열의 경우 -o T를 하여도 문제 없지만, NCBI의 form을 완벽하게 따르지않았다면 F 값을 사용.

blastall -p SELET_PROGRAM -i INPUT_FILE -d DATABASE_FILE -o OUTPUT_FILE  -m  OUTPUT_FORMAT

- p 어떤 blast 프로그램을 사용할지 선택
- i Query가 될 서열, 하나의 서열 혹은 다수의 서열이 하나의 파일에 존재 할 수 있다.
- d blast할 database
-o blast를 수행한 후 결과를 저장할 파일 이름을 지정한다.
-m -m 옵션을 지정하지 않으면 NCBI에서 blast를 수행한 화면을 볼 수 있다. 만약 다른 정보들은 필요 없고, 어느 서열이 어떤 서열과 유사성이 있는지 환인 할 수 있는 정보만 필요하다면 -m의 옵션값을 조절하여 원하는 정보만 저장 할 수 있다. 본인의 경우 -m 8을 많이 애용한다.