[TaxonKit] BLAST결과에서 전체 계통 분류 정보 가져오기

Bioinformatics/기타

[TaxonKit] BLAST결과에서 전체 계통 분류 정보 가져오기

김해김씨99대손 2024. 8. 14. 15:55

🟦 BLAST 옵션의 단점

Local BLAST의 결과에서 우리가 확인할 수 있는 정보는 Tax id와 Phylum, Species이름 정보 등이다. 그러나 전체 계통 정보를 얻는 옵션은 기본적으로 추가되어있지 않다.

나는 기본적으로 blast 수행 시 아래와 같은 output format옵션을 선택한다.

-outfmt "7 delim=, qacc sacc evalue bitscore qcovus pident sscinames"

결과는 아래와 같이

"query/ 접근번호/E-value/bitscore/query coverage/ identity/ 종속명"을 나타냅니다.

>> outfmt의 옵션 전체 보기

약어	내용
qseqid	Query Seq-id (쿼리 시퀀스 ID)
qgi	Query GI (쿼리 GI)
qacc	Query accession (쿼리 액세션)
qaccver	Query accession.version (쿼리 액세션.버전)
qlen	Query sequence length (쿼리 시퀀스 길이)
sseqid	Subject Seq-id (대상 시퀀스 ID)
sallseqid	All subject Seq-id(s), separated by a ';' (모든 대상 시퀀스 ID, ';'로 구분)
sgi	Subject GI (대상 GI)
sallgi	All subject GIs (모든 대상 GI)
sacc	Subject accession (대상 액세션)
saccver	Subject accession.version (대상 액세션.버전)
sallacc	All subject accessions (모든 대상 액세션)
slen	Subject sequence length (대상 시퀀스 길이)
qstart	Start of alignment in query (쿼리에서 정렬 시작 위치)
qend	End of alignment in query (쿼리에서 정렬 끝 위치)
sstart	Start of alignment in subject (대상에서 정렬 시작 위치)
send	End of alignment in subject (대상에서 정렬 끝 위치)
qseq	Aligned part of query sequence (정렬된 쿼리 시퀀스 부분)
sseq	Aligned part of subject sequence (정렬된 대상 시퀀스 부분)
evalue	Expect value (기대값)
bitscore	Bit score (비트 점수)
score	Raw score (원시 점수)
length	Alignment length (정렬 길이)
pident	Percentage of identical matches (동일 일치 비율)
nident	Number of identical matches (동일 일치 수)
mismatch	Number of mismatches (불일치 수)
positive	Number of positive-scoring matches (긍정 점수 일치 수)
gapopen	Number of gap openings (갭 열기 수)
gaps	Total number of gaps (총 갭 수)
ppos	Percentage of positive-scoring matches (긍정 점수 일치 비율)
frames	Query and subject frames separated by a '/' (쿼리 및 대상 프레임, '/'로 구분)
qframe	Query frame (쿼리 프레임)
sframe	Subject frame (대상 프레임)
btop	Blast traceback operations (BTOP) (BLAST 추적 작업)
staxid	Subject Taxonomy ID (대상 계통학적 ID)
ssciname	Subject Scientific Name (대상 학명)
scomname	Subject Common Name (대상 일반명)
sblastname	Subject Blast Name (대상 BLAST 이름)
sskingdom	Subject Super Kingdom (대상 초왕국)
staxids	Unique Subject Taxonomy ID(s), separated by a ';' (고유 대상 계통학적 ID, ';'로 구분)
sscinames	Unique Subject Scientific Name(s), separated by a ';' (고유 대상 학명, ';'로 구분)
scomnames	Unique Subject Common Name(s), separated by a ';' (고유 대상 일반명, ';'로 구분)
sblastnames	Unique Subject Blast Name(s), separated by a ';' (고유 대상 BLAST 이름, ';'로 구분)
sskingdoms	Unique Subject Super Kingdom(s), separated by a ';' (고유 대상 초왕국, ';'로 구분)
stitle	Subject Title (대상 제목)
salltitles	All Subject Title(s), separated by a '<>' (모든 대상 제목, '<>'로 구분)
sstrand	Subject Strand (대상 스트랜드)
qcovs	Query Coverage Per Subject (대상별 쿼리 커버리지)
qcovhsp	Query Coverage Per HSP (HSP별 쿼리 커버리지)
qcovus	Query Coverage Per Unique Subject (고유 대상별 쿼리 커버리지, blastn 전용)

이를 해결하고자 검색해 본 결과 아래와 같은 문의글을 찾을 수 있었다.

https://www.seqanswers.com/forum/bioinformatics/bioinformatics-aa/35570-which-blast-output-format-can-give-full-taxonomony-information

안타깝게도 이에 대한 해결책과 본인이 만든 perl 스크립트를 공유해 준 jameslz 님이 있지만, 스크립트 파일은 현재 다운 불가능이다.

나중에 perl좀 배워봐야겠다.

🟦 TaxonKit 소개

이를 극복하기 위하여 새로운 도구를 사용해 보았다. 바로 TaxonKit - NCBI Taxonomy Toolkit이다.

- 공식 홈페이지: https://bioinf.shenwei.me/taxonkit/

- 논문: Shen, W., & Ren, H. (2021). TaxonKit: A practical and efficient NCBI taxonomy toolkit. Journal of genetics and genomics = Yi chuan xue bao, 48(9), 844–850. https://doi.org/10.1016/j.jgg.2021.03.006 (2024.08.15 기준 인용 232회)

- Github: https://github.com/shenwei356/taxonkit

Shen, W., & Ren, H. (2021). TaxonKit: A practical and efficient NCBI taxonomy toolkit. Journal of genetics and genomics = Yi chuan xue bao , 48 (9), 844–850. https://doi.org/10.1016/j.jgg.2021.03.006

이 도구에서 제공하는 기능을 추리자면 아래와 같다.

1) taxonkit을 사용하여 tax id로 계통 정보를 추출한다

2) reformat으로 추출한 계통 정보를 subfamily, subgenus, stain 등 제거하고 깔끔하게 kingdom, phylum, class, order, family, Genus, Species만 추출 후 정리해 준다.

3) name2taxid를 사용하여 반대로 species이름으로 tax id를 추출할 수 있다.

🟦 설치

- https://bioinf.shenwei.me/taxonkit/download/

- 설치 방법은 1) 직접 다운로드 2) conda사용 3) home brew를 사용할 수 있다.

이 중에서 conda가 가장 최신 버전을 지원해 줌으로 conda로 설치해 보자.

conda install -c bioconda taxonkit
taxonkit -h # 설치 완료

그다음에는 계통 정보를 추출할 NCBI의 taxonomy database를 다운로드하여 보자.

cd ~/.taxonkit # 없다면 mkdir -p $HOME/.taxonkit 로 디렉토리 생성하기
wget ftp://ftp.ncbi.nih.gov/pub/taxonomy/taxdump.tar.gz
tar -xvzf taxdump.tar.gz # 압축풀기

🟦 TaxonKit를 사용한 blast결과 정리

1) blast 돌리기

기존 blast 서열 output옵션에 tax id를 추가로 얻기 위하여 staxids를 추가해 보자,

blastn -db /data/Reference/BLASTDB/ITS_RefSeq_Fungi  \
       -query sequence.fasta \
       -task blastn \
       -dust no \
       -outfmt "7 delim=, qacc staxids sacc evalue bitscore qcovus pident sscinames" \
       -max_target_seqs 1 >   blast-output
grep -v '#'  blast-output > blast-output_summary

그렇다면 결과물의 두 번째 열에서 tax id를 찾을 수 있을 것이다.

2) tax_id 추출하기

두 번째 열에 존재하는 tax id를 추출하고 이를 저장해 보자.

awk -F',' '{print $2}'blast-output_summary > blast-output_taxid  # tax_id만 추출
sed 's/.*;//' blast-output_taxid > blast-output_taxid2           # 간혹 tax_id가 ;로 두 숫자가 이어진 경우가 있음, 이때 두 번째 숫자만 추출

3) TaxonKit를 사용하여 전체 계통 정보 확인하기

taxonkit를 사용하여 lineage 정보를 추출하고, 이를 reformat을 사용하여 정리해 준다.

# tax_id에서 계통 정보 전체 추출하기 
taxonkit lineage blast-output_taxid2 | awk '$2!=""' > blast_lineage.txt

# 추출된 정보 정리하기 
cat blast_lineage.txt   \
   | taxonkit reformat -F   \
   | csvtk -H -t cut -f 1,3     \
   | csvtk -H -t sep -f 2 -s ';' -R     \
   | csvtk add-header -t -n taxid,kindom,phylum,class,order,family,genus,species     \
   | csvtk pretty -t

taxonkit reformat -F에서 -F옵션은 빈칸을 "unclassified + 상위 tax level"으로 바꾸어 주는 기능입니다.

아래와 같은 결과물을 얻을 수 있습니다.