Protein sequence clustering과 dN/dS 분석

새로운 과제를 시작하면서 단백질 서열 클러스터링이라는 매우 고전적인 일을 하게 되었다. 단백질 패밀리의 분석, pan-genome 계산, ortholog 추출, dN/dS 분석, 데이터셋의 크기 축소 등 단백질 클러스터링은 매우 여러 분야에 쓰인다. 자료를 조사하면서 지난 2월에 작성했던 글인 Usearch 맛보기도 업데이트하였다. 오늘 작성하는 글은 "How to cluster protein sequences: tools, tips and commands"을 많이 참고하였다.

단백질 서열 클러스터링 기법은 크게 다음의 두 가지로 나눌 수 있다.

• Sequence-based clustering: CD-HIT, UCLUST, USEARCH...
• Graph-based clustering: LAST, legacy BLAST, BLAST+...
• Network clustering: MCL

Graph-based clustering에서는 모든 단백질 서열 사이의 1:1 similarity 정보를 얻어내는 것이 필요하다. BLAST 계열의 프로그램을 사용하는 경우 결과물을 파싱하여 Reciprocal Best BLAST hit(RBB, RBH, or biodirectional best hit BBH)을 추출하여 활용하면 되는 것이다. RBB는 두 genome 사이에서 ortholog candidate를 산출하는 데에도 널리 이용된다.

RBH으로부터 ortholog를 정확하기 추출하기 위한 BLAST option을 제안한 논문이 2008년에 발표된바 있다(PMID: 18042555). 이를 기반으로 하여 다음의 BLASTP 옵션이 제안되었다. 논문에서는 첫번째 조건만을 제안하였다. 이렇게 함으로써 ortholog의 수는 가장 많게, 그러나 error는 최소화할 수 있었다고 한다.

• -F "m S" -s T (combination of soft filtering and Smith-Waterman final alignment)
• E value: 10^-5 또는 10^-6 이하

BLAST 프로그램은 database search와 alignment의 두 단계로 이루어지는데, soft filtering이란 첫번째 단계에서만 low-complexity sequence를 마스킹하는 것을 의미한다(-F "m S"). 기본 조건은 -F T -s F이다. -F "m S"라는 옵션이 있다는 것은 blastall의 모든 옵션이 수록된 페이지(링크; 혹은 명령행에서 blastall이라고만 치면 나오는 텍스트)에서 찾아보기 어렵다.

이상의 옵션은 legacy BLAST에 대한 것이다. 그러면 요즘 권장되는 BLAST+에서는 무엇이라고 옵션을 제공해야 하는가? NCBI의 매뉴얼 페이지(링크)를 참조해 보면 -soft_masking T -use_sw_tback (-evalue 10^-5)일 것으로 추청된다. 

다음의 북 챕터에서는 BLAST+를 이용하여 상동유전자 정보를 찾고, 이로부터 dN/dS 분석일 실시하는 방법을 친절하게 설명하였다. 사용된 Perl script(BLAST+ 결과물에서 RBB를 추출하는 스크립트 포함)는 제1저자인 Daniel Jeffares의 웹사이트에서 입수 가능하다고 해 놓았다.

A Beginners Guide to Estimating the Non-synonymous to Synonymous Rate Ratio of all Protein-Coding Genes in a Genome. Methods Mol Biol. 2015;1201:65-90. doi: 10.1007/978-1-4939-1438-8_4. PMID.

그러나 지금 이 사이트를 방문해보면 figshare에 Perl 스크립트 일부만이 남아있다. 그러나 Daniel은 이를 기반으로 하여 VESPA라는 파이프라인을 개발하여 2017년 PeerJ Comput. Sci.에 발표하였다(VESPA: Very large-scale Evolutionary and Selective Pressure Analyses. 저널 GitHub 매뉴얼). VESPA를 공부하려면 단순히 RBB를 찾으려는 목적이 아니라 더 큰 안목을 가지고 착수해야 할 것이다.

InParanoid를 사용해서 두 단백질 세트 사이에서 ortholog를 찾을 수도 있다. 그런데 Why InParanoid sucks라는 글(2014년도)을 보니 이제는 이 프로그램을 떠나야 하는 것이 아닐까하는 생각도 들기 시작하였다. 알고리즘적으로 특별히 문제가 되는 것은 아니지만 말이다. 어쩌면 blast_rbh.py 스크립트가 바로 내가 찾던 것이 아닌지 모르겠다. 이를 다운로드하여 실행해 보았다. duplicate가 있다는 경고는 나오지만 결과 파일에서 이들이 어느 것인지를 보여주지는 않는다. 상동유전자를 찾는 프로그램으로서 GET_HOMOLOGUES도 유용한 프로그램으로 여겨진다.

$ ls
A.faa  B.faa
$ python ~/script/blast_rbh.py -o example.tab -t blastp -a prot A.faa B.faa 

Done, 10 RBH found
Warning: Sequences with tied best hits found, you may have duplicates/clusters
$ ls
A.faa  B.faa  example.tab
$ head -n 1 example.tab 
#A_id B_id A_length B_length A_qcovhsp B_qcovhsp length pident bitscore

BLAST+ 분석 결과물에서 best hit(RBB는 아님)을 추출한 뒤 mcl로 클러스터링하는 명령행을 정리하여 보았다. 이는 서두에서 언급한 참고문헌에서 발췌한 것이다.

$ makeblastdb -in DB.fasta -dbtype prot -out DB

$ blastp -query example.fasta -db DB -outfmt '6 std qlen slen' -out output_tmp.list -evalue 1.03e-05

$ awk '{if($1!~/^#/ && $3>=50) print $1"\t"$2"\t"($3/100)}' output_tmp.list > hits.list

$ sort -nrk 3 hits.list | sort -k 1,2 -u > best_hit.list

$ mcl hits.list -I 1.8 --abc -o clusters.txt

$ awk '{if(NF>=2) print}' clusters.txt > clusters_after_threshold.txt

맨 아래의 awk 명령어는 singleton을 제거하기 위함이다.

마지막으로 대용량의 pan-genome 분석을 위한 새로운 프로그램을 하나 소개해 본다. 늘 Roary만 편식을 심하게 했었는데 이제는 다른 것에도 눈을 돌려볼 때이다.

BPGA - an ultra-fasta pan-genome analysis pipeline. Sci. Rep. (2016).
Piggy: a rapid, large-scale pan-genome analysis tool for intergenic regions in bacteria. GigaScience (2018)