Haplomerger #

Haplomerger는 polymorphic 한 diploid 어셈블리 유전체의 allele relationship을 reconstructing 하는 프로그램이다. diploid 에 해당하는 어떤 종에 대한 assembly를 수행했을 때, 예상했던 contig의 수보다 많게 assembly되어 contig가 실제보다 더 잘게 분할되었을 경우를 대비하여 diploid 어셈블리 유전체를 self-alignment를 통해 overlap 되는 contig의 수를 줄일 수 있다.

논문 #

Genome Res. 2012 Aug;22(8):1581-8. doi: 10.1101/gr.133652.111. Epub 2012 May 3.

Haplomerger pipeline #

Haplomerger는 다음과 같은 파이프라인으로 구성되어 있다.

[출처 : Genome Res. 2012 Aug;22(8):1581-8. doi: 10.1101/gr.133652.111. Epub 2012 May 3.]

파이프라인에 대한 내용을 간단히 살펴보면, Haplomerger의 첫 번째 단계는 diploid assembly 유전체의 각 contig들을 pairwise하게 self-alignment를 수행하게 된다. 먼저 assembly유전체 파일을 두 개의 copies로 duplication 시킨 후에 하나는 query로, 다른 하나는 target으로 이름을 정한 후에 local alignment프로그램인 LASTZ를 통해 all-against-all alignment를 수행하게 된다.

[출처 : Genome Res. 2012 Aug;22(8):1581-8. doi: 10.1101/gr.133652.111. Epub 2012 May 3.]

query와 target의 alignment에 대한 두 가지 시나리오가 존재한다. 먼저 query A가 target B에는 alignment가 되지만, query B가 target A에 alignment되지 않거나, query B가 target A에는 alignment가 되지만 query A에는 Hit되지 않는 경우 등을 통해 각 contigs끼리 잇는 모든 alignment에 scoring이 매겨지게 되고 각 contig에 Highest scored alignment들을 잇는 DGA(Diploid Genome Assembly) 그래프가 형성되게 된다.
다음으로 DGA 그래프에서 bubble을 형성하거나, 중복되는 edges 들을 제거하여 각 simplified한 graph가 형성되게 하고, 그 그래프를 연결하는 scaffold connection 과정을 거치게 됩니다. 마지막으로 연결된 scaffold들로 만들어진 diploid 유전체의 각 allele pair에 representative한 allele sequence를 선택하여 reference를 출력하게 된다. 또한 Haplomerger는 representative allele 내에 있는 각 N gap에 대해서 alternative allele sequence로 대체하는 과정을 거쳐서 alternative를 출력하게 된다.

명령어 #

Haplomerger는 다음과 같은 7개의 shell script로 구성되어 있다.(hm.batchA-G).

hm.batchA.initiation_and_all_lastz (create environment and run Lastz)
hm.batchB.chainNet_and_netToMaf (run chainNet on lastz alignments)
hm.batchC.haplomerger (run haplomerger)
hm.batchD.refine_removing_tandem_assemblies
hm.batchE.refine_haplomerger_updating_nodes_and_portions
hm.batchF.refine_haplomerger_connections_and_Ngap_fillings
hm.batchG.refine_unincorporated_sequences

이 스크립트들은 batch 알파벳 순서대로 차례차례 수행되어야 하므로 InstallDir폴더로 이동하여 project_template이라는 폴더를 Haplomerger를 수행하고자 하는 위치로 복사하여 사용하는 것이 가장 쉬운 방법이다. 또한 Merging하기를 원하는 diploid file의 경우 또한 name을 genome.fa.gz 형태로 바꾸어주게 되면 아래의 내용과 같은 run_all.batch 프로그램을 이용하여 한번에 처리가 가능하다.

#!/bin/bash

#### =============================================================
#### set the path to executables  // 환경설정
#### =============================================================
#PATH=/path/to/perl/bin/:/path/to/chainNet/:/path/to/lastz/:/path/to/gapCloser/:/path/to/sspace/:$PATH
#echo $PATH

#### =============================================================
#### lift the file handle limit (require root privilege) 
#### =============================================================
ulimit -n 655350
반드시 이 설정을 해주어야 한다.
보통 리눅스의 경우, 한 프로세스의 파일 열기를 1,024개로 제한해 놓았는데, 이 설정을 해주지 않으면 Haplomerger 수행에 에러를 발생할 수 있다.
(can't open genome.seq/~~.nib: 등과 같은 에러가 발생한다)

#### =============================================================
#### step 1: break the misjoins and output the new diploid assembly
#### =============================================================
#### default input:        ./{assembly_name}.fa.gz
#### default final output: ./{assembly_name}_mp.fa.gz

./hm.batchA1.initiation_and_all_lastz  genome
./hm.batchA2.chainNet_and_netToMaf     genome
./hm.batchA3.misjoin_processing        genome

#### =============================================================
#### step 2: separate/merge two alleles and output haploid assemblies
#### =============================================================
#### default input:        ./{assembly_name}_mp.fa.gz
#### default final output: ./{assembly_name}_mp_ref.fa.gz
#### default final output: ./{assembly_name}_mp_alt.fa.gz

./hm.batchB1.initiation_and_all_lastz             genome_mp
./hm.batchB2.chainNet_and_netToMaf                genome_mp
./hm.batchB3.haplomerger                          genome_mp
./hm.batchB4.refine_unpaired_sequences            genome_mp
./hm.batchB5.merge_paired_and_unpaired_sequences  genome_mp

#### =============================================================
#### =============================================================
#### step 3: further scaffold the obtained haploid assemblies
#### =============================================================
#### default input:        ./{assembly_name}_mp_ref.fa.gz
#### default final output: ./{assembly_name}_mp_ref_ss.fa.gz
#### default final output: ./{assembly_name}_mp_alt_ss.fa.gz

#./hm.batchC1.hierarchical_scaffolding                      genome_mp_ref
#./hm.batchC2.update_reference_and_alternative_assemblies   genome_mp_ref

#### =============================================================
#### step 4: remove tandem errors from haploid assemblies
#### =============================================================
#### default input:        ./{assembly_name}_mp_ref_ss.fa.gz
#### default final output: ./{assembly_name}_mp_ref/alt_ss_rt.fa.gz

#./hm.batchD1.initiation_and_all_lastz     genome_mp_ref/alt_ss
#./hm.batchD2.chainNet_and_netToMaf        genome_mp_ref/alt_ss
#./hm.batchD3.remove_tandem_assemblies     genome_mp_ref/alt_ss

#### =============================================================
#### step 5: try to close Ngaps in haploid assemblies
#### =============================================================
#### default input:        ./{assembly_name}_mp_ref_ss_rt.fa.gz
#### default final output: ./{assembly_name}_mp_ref/alt_ss_rt_gf.fa.gz

#./hm.batchE1.wrapper_for_gapCloser_v1.12  genome_mp_ref/alt_ss_rt