LS-BSR(large scale BSR)은 대량의 유전체 서열을 대상으로 하여 정해진 마커 유전자(염기 혹은 아미노산 서열)이 존재하는지 여부를 확인하는 프로그램이다. BSR(BLAST Score Ratio)란 (검색 대상 유전체에서 찾아진 hit의 BLAST bit score)/(자기 자신에 대한 BLAST bit score)로 표현되는 비율로서 0에서 1 사이의 값을 가지며, BSL이 1에 가까울수록 유전체 내에 해당 마커 유전자가 존재할 가능성이 높다고 볼 수 있다. 유전체 서열은 annotation이 되지 않은 FASTA 파일로 주어져야 하며, 마커 유전자를 제공하지 않으면 gene prediction을 먼저 수행하여 이를 상호 비교하여 작업을 개시한다. 따라서 pan-genome과 관련한 계산치를 제공한다는 측면에서 Roary와 비교할 수 있으나 성격은 좀 다르다. 두 프로그램의 비교에 대해서는 Some thoughts on comparing Roary to LS BSR을 참조하기 바란다. 결과물로 주어지는 BSR matrix는 취향에 맞는 third party program을 이용하여 시각화할 수 있다.

꽤 많은 prerequisite가 있는 응용프로그램이지만 bioconda 채널에서 지원하지는 않는다. 따라서 약간 번거롭게 설치해야 한다. 다음의 링크를 참조하라.

실습용 서버에 프로그램 설치하기(ls_bsr environment)

유전체 서열 파일은 반드시 fasta 확장자를 가져야 하며, 파일이 위치한 디렉토리를 -d DIRECTORY 옵션으로 지정해야 한다. 마커 서열 파일은 염기(.fasta) 혹은 아미노산(.pep) 수록 정보에 따라 확장자를 정해 놓아야 한다. -x PREFIX는 결과 파일에 부여되는 prefix를 지정할 때 쓴다. 유전체 염기서열은 ani_r_exercise.zip 압축파일에 포함된 accessions.txt의 것을 ncbi-genome-download로 다운로드하여 사용하자. 편의를 위하여 유전체 염기서열 파일을 업로드하였다(51_paenibacillus_genomes.zip).마커 유전자 염기서열은 markers.zip의 압축을 풀어서 생성되는 markers.fa파일을 markers.fasta로 이름을 변경하여 사용한다. 서열 ID의 형식이 'gene|locus_tag|protein_id'라서 매우 긴 편인 데다가 비어있는 필드는 '-'로 채워져 있지만, ls_bsr의 작동에는 이상이 없었다. 실습에 쓰인 유전체 및 마커 유전자 염기서열은 나의 2019년 논문 Chronicle of a soli bacterium: Paenibacillus polymyxa E681 as a tiny guardian of plant and human health(링크)의 그림 2를 작성하기 위해 쓰였던 것이다.

$ python ~/apps/LS-BSR/ls_bsr.py -p 4 -d genomes -g markers.fasta -x RUN
LOG: 2023/08/16 17:33:42 - Testing paths of dependencies
/home/hyjeong/miniconda3/envs/ls_bsr/bin/tblastn
citation: Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, and Lipman DJ. 1997. Gapped BLAST and 
PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res 25:3389-3402
LOG: 2023/08/16 17:33:42 - Using pre-compiled set of predicted genes
LOG: 2023/08/16 17:33:42 - using tblastn
LOG: 2023/08/16 17:33:43 - starting BLAST
LOG: 2023/08/16 17:36:42 - BLAST done
LOG: 2023/08/16 17:36:42 - Duplicate searching turned off
LOG: 2023/08/16 17:36:42 - starting matrix building
LOG: 2023/08/16 17:36:42 - matrix built
LOG: 2023/08/16 17:36:43 - all Done
$ ls -l RUN*
-rw-rw-r-- 1 hyjeong hyjeong 19196  8월 16 17:36 RUN_bsr_matrix.txt
-rw-rw-r-- 1 hyjeong hyjeong   816  8월 16 17:36 RUN_names.txt
-rw-rw-r-- 1 hyjeong hyjeong   212  8월 16 17:36 RUN_run_parameters.txt

3개 결과 파일을 묶어서 업로드해 두었다(ls-bsr_results_20230816.zip).

마커 유전자를 제공하지 않으면 de novo gene prediction을 먼저 실시한 다음 클러스터링을 수행하여 conserved + unique gene의 세트를 만든 뒤 이를 대상으로 BSR matrix를 생성한다. 클러스터링에 사용할 프로그램은 usearch, vsearch 또는 cd-hit 중에서 고를 수 있으며(-c CLUSTER_METHOD) 전부 PATH에 있어야 한다.

path_to_LS-BSR/tools 디렉토리에 있는 20개 부속 스크립트의 사용법과 pan-genome analysis를 위한 응용법, 결과물의 시각화 방법 등에 대해서는 매뉴얼을 참조하라. 부속 스크립트는 전부 BSR matrix를 입력물로 이용한다. 예를 들어 isolate_uniques_BSR.py는 단일 게놈이 유일하게 갖고 있는 CDS를 출력해 준다. 꽤 쓸 만한 스크립트가 많이 있는데 미처 다 테스트를 해 보지는 못하였다.

다음의 R 코드는 gplots의 heatmap.2() 함수를 이용하여 BSR 매트릭스를 시각화하는 사례이다. 보기 좋은 그림을 얻기 위해서는 부수적인 작업이 필요하다.

> library('gplots')
> bsr=read.table("bsr_matrix.txt", sep='\t', header=T, row.names=1)
> bsr.t = t(bsr)
> View(bsr.t) # 생략
> dim(bsr.t)
[1] 51 48
> heatmap.2(as.matrix(bsr.t), col=greenred(10), trace="none", scale="none", dendrogram=c("row"))

gplots 패키지에서 heatmap.2() 함수로 그린 plot의 레이아웃을 조절하는 것은 또 다른 문제이다. 위 R 코드를 실행하면 다음 그림과 같이 column/raw label의 일부가 잘린 그림을 얻게 된다.

gplots보다는 ComplexHeatmap 패키지를 이용하는 것이 더 낫다고 생각한다. 자세한 것은 별도로 작성한 글인 R을 이용하여 ANI matrix 자료 다루기-Heatmap 그리기를 참고하기 바란다.