Running The Code:

Put the "Aquila/bin" in the ".bashrc" file, and source the ".bashrc" file
Or use the fullpath of "Aquila_assembly_based_variants_call.py" and "Aquila_phasing_all_variants.py"

Aquila/bin/Aquila_assembly_based_variants_call.py --assembly_dir Assembly_results_S12878 --out_dir Variants_results_S12878 --ref_file Aquila/source/ref.fa 

*Required parameters

--assembly_dir: "Assembly_results_S12878" is the output folder from Aquila assembly steps. (Here, Aquila will use contig files from "Assembly_results_S12878/Assembly_Contigs_files" to perform variants calling.)

--ref_file: "Aquila/source/ref.fa" is the reference fasta file from GRCh38, you can download it by running "./install.sh"

*Optional parameters

--var_size: cut off size for indel and SV, default = 1.

--num_of_threads: default = 1. To increase "--num_of_threads", large memory is needed (see the below Memory/Time Usage table). Maximum is 23 which corresponds to the number of chromosomes.

--out_dir: default = ./Aquila_Variant_Results

--chr_start --chr_end: if you only want to call variants from some chromosomes or only one chromosome. For example: use "--chr_start 1 --chr_end 5" will call variants from chromsomes 1,2,3,4,5. Use "--chr_start 2 --chr_end 2" will only call variants from chromosome 2.

--clean_flag: default = 0 (keep all intermediate files). Use "--clean_flag 1" to clean all intermediate files.

--all_regions_flag: default = 0 (variants calling only for diploid assemblies). For compound SVs, we randomly select one alternate allele(s) for the "ALT" field of the VCF file. Use "--all_regions_flag 1" for variants calling in all regions. For compound SVs, we split one compound SV into two heterozygous SVs. You can process compound SVs by your own definition based on the coordinates from the final VCF file.

*Final Output

Aquila_final_sorted.vcf: it includes SNPs, indels, and SVs. (This coordinate of this VCF file starts from 0)

:octocat: To evaluate snp/indel/SV from "Aquila_final_sorted.vcf" with GIAB Gold standards, check here in details.

Memory/Time Usage for "--all_regions_flag 1" (This mode is used for reporting variants calling in Aquila paper)

Memory	--num_of_threads	Time for WGS
350GB	4	01:08:14
200GB	1	02:12:44

Memory/Time Usage for "--all_regions_flag 0"

Memory	--num_of_threads	Time for WGS
100GB	1	1-10:24:32
150GB	2	17:35:45
250GB	3	13:25:06

2. Phasing all types of variants :eagle:

Aquila/bin/Aquila_phasing_all_variants.py --assembly_vcf Variants_results_S12878/Aquila_final_sorted.vcf --vcf_file S12878_freebayes.vcf --assembly_dir Assembly_results_S12878

*Required parameters

--assembly_vcf: "Variants_results_S12878/Aquila_final_sorted.vcf", it is generated by the above step.

--vcf_file: "S12878_freebayes.vcf", it is the origin vcf file generated by FreeBayes.

--assembly_dir: "Assembly_results_S12878" is the output folder from Aquila assembly steps. (Here, Aquila will use phasing information from "Assembly_results_S12878/results_phased_probmodel", and high-confidence partitioning points profile from "Assembly_results_S12878/phase_blocks_cut_highconf", to perform variant calling.)

*Optional parameters

--out_dir: default = ./Aquila_Phasing_Results

--chr_start --chr_end: if you only want to call variants from some chromosomes or only one chromosome. For example: use "--chr_start 1 --chr_end 5" will call variants from chromsomes 1,2,3,4,5. Use "--chr_start 2 --chr_end 2" will only call variants from chromosome 2.

--block_len_use": default = 100000, it is the same "--block_len_use" used in assembly steps.

*Final Output

Aquila_all_phased_sorted.vcf: It includes phased SNPs, indels, SVs, all homozygous variants, and a small set of unphased variants. (This coordinate of this VCF file starts from 1)

Memory/Time Usage

Memory	Time for WGS
100GB	05:22:10