일단 Zakour의 다음 자료를 정독하라. detectection_of_recombinational_events_in_bacterial_genomes.pdf
Xabier Didelot는 microevolution 과정에서 recombination이 중요함을 인식하고 이를 검출하는 ClonalFrame이라는 프로그램을 개발하였다(논문 웹사이트 ). 발표 당시에는 multilocus sequence data에 국한되었지만 현재는 박테리아 유전체에도 적용 가능한 ClonalFrameML(GitHub)도 개발되었다. 고전적인 논문 Rapid pneumococcal evolution in response to clinical interventions (Science 2011, PubMed)를 읽어보면 박테리아의 유전체에서 벌어지는 homologous recombination이 얼마나 중요한지를 알 수 있다. 그러나 모든 미생물 종에서 같은 정도로 recombination이 일어나는 것은 아님에 유의하자.
이 페이지에서는 Gubbins(Genealogies Unbiased By recomBinations in Nucleotide Sequences, 논문 웹사이트)을 이용하여 미생물 유전체로부터 recombination이 일어난 곳을 찾아내고 이를 시각화하는 방법을 알아본다.
Snippy(GitHub)는 raw sequencing read와 contig를 전부 사용할 수 있는 유연한 도구이다. 그러나 설치 과정에서 에러가 발견되어 이를 극복하는 과정을 상세히 기술하고자 한다.
Snippy는 BioConda에 포함되어 있어서(최신 버전 4.0_dev) 설치는 매우 쉬웠다. 별도의 환경을 만들지 않고 Base에서 설치하면 된다. 설치를 마친 뒤 모든 프로그램이 다 깔린 상태인지를 snippy –check 명령으로 확인을 하니 vcfirstheader가 없다는 메시지가 나왔다. 이 유틸리티는 vcflib(GitHub)에 포함되어 있으므로 git clone 명령을 사용하여 수작업으로 설치하였다. 그리고 설명 파일에서 소개한 대로 <PathWhereIWantToCloneVcflib>/vcflib/bin/을 PATH 환경변수에 추가하였다.
그런데 snippy의 첫 실행에서 다음과 같은 에러가 발생하였다.
java.lang.RuntimeException: Error parsing property 'ref.AE015924.codonTable'. No such codon table 'Bacterial_and_Plant_Plastid' ...
이에 대해서는 이미 snippy issue 페이지에 154번째 문서로 보고된 상태이다(링크). 2018년 4월부터 시작된 쓰레드이므로 매우 최신의 에러라고 할 수 있다. 이에 대하여 3.2 버전을 쓰라는 조언이 있었다. 더욱 바람직하게는 4.0 안정화 버전이 빨리 나와야 한다.
conda remove snippy conda install -c bioconda -c conda-forge snippy=3.2
이렇게 한 뒤 다시 테스트 러닝을 하였으나 제대로 결과가 나오지 않았다. snps.log 파일을 열어보니 이번에는 다음과 같은 내용이었다.
/data/anaconda2/bin/freebayes-parallel: line 40: ../vcflib/bin/vcfstreamsort: 그런 파일이나 디렉터리가 없습니다 /data/anaconda2/bin/freebayes-parallel: line 40: ../vcflib/scripts/vcffirstheader: 그런 파일이나 디렉터리가 없습니다
/data/anaconda2/bin/freebayes-parallel의 40번째 라인을 찾아가서 위 두 가지 유틸리티 앞에 붙은 상대경로를 지워버렸다. 이렇게 하니 비로소 정상적인 실행이 되었다. 일반적인 실행 방법은 다음과 같다.
$ snippy --cpus 16 --outdir mysnps --ref E_coli.gbk --R1 read_1.fastq.gz --R2 read_2.fastq.gz $ snippy-core --prefix core mysnps1 mysnps2 mysnps2
snippy-core 결과물인 core.full.aln(whole genome SNP alignment including invariant sites)을 gubbins에 이용하면 된다.
Harvest suite는 동일 종에 속하는 미생물 유전체를 신속하게 비교하여 core-genome alignment를 만들고 시각화하는 도구이다.
Parsnp가 기본적으로 생성하는 출력파일은 다음의 세 가지이다. Gingr에서 .ggr 파일을 열면 MFA 파일로 export할 수 있다. XMFA(eXtended Multi FastA) 파일의 포맷은 여기에 설명되어 있다.
parsnp.ggr parsnp.tree parsnp.xmfa
HarvestTools 유틸리티는 좀 더 다양한 포맷간의 전환을 도와준다. Gingr 혹은 harvesttools에서 전환 가능한 파일 포맷에 대한 상세한 내용은 여기를 참고하라.
$ harvesttools -i parsnp.gingr -M multiFASTA.fna # (concatenated LCBs)
이렇게 만들어진 multiFASTA 파일을 gubbins에 입력물로 제공해도 실행에는 문제가 없다. 그러나 LCB(locally collinear block)을 갭 없이 이어붙인 것이라서 원래의 좌표 정보가 더 이상 유효하지 않다. parsnp.xmfa를 정밀하게 parsing하면 reference coordinate에 맞추어서 각 LCB 사이에 알맞은 길이의 gap을 넣어서 연결한 서열을 만드는 것이 가능할지도 모른다. 하지만 그렇게 하느니 처음부터 whole genome alignment를 만드는 snippy를 쓰는 것이 바람직할 것이다.
BioConda 패키지를 이용하여 설치하였다(py35 environment).
$ run_gubbins.py dataset.aln $ gubbins_drawer.py -o result.pdf dataset.final_tree.tre dataset.recombination_predictions.embl $ gubbins_drawer.py -o results_SNPS.pdf dataset.final_tree.tre dataset.branch_base_reconstruction.embl