SMRT analysis의 HGAP(RS_HGAP_assembly.2 or 3)에서는 1번의 resequencing을 실행하지만 error가 남아있는 경우가 많다. 이것을 그대로 GenBank에 제출하면 pseudogene이 많다는 이유로 반송되기도 한다. 보통 전체 유전자의 10% 이상인 경우 문제가 있다고 제출자에게 되돌려 보내는 것으로 생각된다. 이런 일이 항상 벌어지는 것은 아니지만, 경험상 HGAP으로 얻은 assembly는 추가적인 교정을 하는 것이 바람직하다고 생각한다. 물론 assembly error까지를 전부 점검할 수는 없는 노릇이다. 여기서 말하는 것은 염기 수준의 에러 교정을 의미한다.

일루미나 read가 있다면 적당한 도구를 이용하여 매핑한 뒤 consensus를 추출하거나, 아니면 pilon(논문 GitHub)으로 후처리를 하는 것을 권장한다. 이것은 PacBio long read의 교정이 아니라 조립을 마친 결과물의 교정임을 잊지 말자. 외부 데이터가 없더라고 PacBio long read를 이용하여 반복적으로 교정을 할 수 있다.

HGAP이 끝난 뒤 DATA→Assembly→Polished Assembly(FASTA)를 다운로드한다. 파일 이름은 다음과 같이 고친다. 꼭 이러한 방식을 따라야 한다는 것은 아니다. 다만 일관성있게 하는 것이 중요하다.

(수정 전)_das_smrtanalysis_userdata_jobs_016_016491_data_polished_assembly.fasta.gz
(압축 해제 후 수정) jobs_016_016491_data_polished_assembly.fasta

이를 Reference sequence로 등록한다(“Manage Reference Sequences”). Name은 “균주명 polished_01” 정도로 하여 업로드한 뒤 이를 template로 하여 RS_Resequencing.1을 실행한다. Reference sequence ID는 균주명_polished_id와 같이 공백과 dash는 전부 밑줄로 바뀌어 있을 것이다. 작업이 끝나면 결과창에서 Analysis→Variants를 확인해 보자. 만약 contig 각각에 대해서 1개 정도의 variant가 발견되어 수정되었다면 이것으로 작업을 종결할 수도 있으나 안전을 위해서 두번 resequencing을 하는 것을 표준으로 삼는다. 이번에는 DATA→ANALYSIS→Consensus(FASTA)를 다운로드하여 압축을 해제한 뒤 파일명을 수정한다.

(수정 전) _das_smrtanalysis_userdata_jobs_016_016493_data_consensus.fasta.gz
(압축 해제 후 수정) jobs_016_016493_data_consensus.fasta

다시 이를 Reference sequence로 등록한다. Name은 “균주명 consensus_01”로 한다. Job ID의 맨 끝자리 숫자가 1→3→5로 증가하고 있다.

각 서열의 ID는 다음과 같이 변한다. (3)을 취하여 최종 서열로 사용하면 된다.

(1) scf7180000000030|quiver
(2) scf7180000000030|quiver|quiver
(3) scf7180000000030|quiver|quiver|quiver

$ midir bowtie2
$ bowtie2-build REFERENCE.fasta bowtie2/refernece
$ bowtie2 -p 24 -x bowtie2/reference -1 ../SRR5032362_1.fastq -2 ../SRR5032362_2.fastq -S mapping.sam
$ samtools view -b -S -o mapping.bam mapping.sam
$ samtools sort [ -T /tmp/aln.sorted ] -o mapping.sorted.bam mapping.bam
$ samtools index mapping.sorted.bam
$ pilon --genome 06.fixstart.fasta --frags mapping.sorted.bam

PacBio long read를 이용하여 염색체 수준의 조립물을 얻는 것은 이제 아주 일상적인 일이 되었다. 세균의 염색체는 원형이므로, 가장 이상적인 경우에는 양 말단에 서로 겹치는 염기서열이 출현하는 linear contig를 얻게 된다. 이를 조정하여 중복을 제거하고, 복제 원점(보통 dnaA를 첫번째 유전자로 삼음)을 기준으로 서열을 조절하는 일이 과제로 남는다(참조: 단일 contig 서열의 후처리)

매우 다행스럽게도 이러한 후처리 작업을 자동으로 실시하는 프로그램인 circlator(Genome Biology 2015 16:294)가 개발되었다. 입력물은 error-corrected long read와 assembly(FASTA)만 주어지면 된다. 프로그램에 대한 자세한 설명은 GitHub site를 참조하라.

http://sanger-pathogens.github.io/circlator/

pip3를 이용하거나, 혹은 최신 소스를 받아서 설치하면 된다. de novo assembler는 canu 혹은 SPAdes(3.7.1 권장)가 필요하다. 사전에 필요한 프로그램을 구성하는 것이 귀찮으면 소스 패키지 안에 들어있는 install_dependencies.sh 스크립트를 실행하는 것으로 충분하다. 그러면 build directory 안에 다음의 서브디렉토리가 생기고 그 하위에 필요한 프로그램이 전부 빌드된다.

MUMmer3.23     SPAdes-3.7.1-Linux     bwa-0.7.12          
canu-1.4       prodigal-2.6.2         samtools-1.3

• https://github.com/sanger-pathogens/circlator/wiki
• (인내심이 부족한 사람들을 위해) https://github.com/sanger-pathogens/circlator/wiki/Brief-instructions

입력물은 contig sequence file과 corrected reads(fasta or fastq)이며, 기본 동작은 all task이다. 이를 실행하면 progcheck, mapreads, bam2reads, assemble, mergem clean 및 fixstart가 단계적으로 이루어진다. 기본 사용법은 다음과 같다. 첫줄의 $PATH 환경변수 설정은 각 상황에 맞게 알아서 실행하라.

(필요한 경우) $ pyenv global 3.5.1
$ PATH=/usr/local/apps/canu/Linux-amd64/bin:/usr/local/apps/SPAdes-3.7.1-Linux/bin:$PATH
$ circlator progcheck (optional)
$ circlator all [options] <assembly.fasta> <corrected_reads.fasta> <output directory>