git clone https://github.com/maiziex/Aquila.git
Aquila utilizes Python3 (+ numpy, pysam, and scipy), SAMtools, and minimap2. To be able to execute the above programs by typing their name on the command line, the program executables must be in one of the directories listed in the PATH environment variable (".bashrc").
Or you could just run "./install.sh" to check their availability and install them if not, but make sure you have installed "python3", "conda" and "wget" first.
cd Aquila
chmod +x install.sh
./install.sh
After running "./install.sh", a folder "source" would be download, it includes human GRCh38 reference fasta file, or you could also just download it by yourself from the corresponding official websites.
Put the "Aquila/bin" in the ".bashrc" file, and source the ".bashrc" file
Or just use the fullpath of "Aquila_step1.py" and "Aquila_step2.py"
Aquila/bin/Aquila_step1.py --bam_file possorted_bam.bam --vcf_file S12878_freebayes.vcf --sample_name S12878 --out_dir Assembly_results_S12878 --uniq_map_dir Aquila/Uniqness_map
--bam_file: "possorted_bam.bam" is bam file generated from barcode-awere aligner like "Lonranger align". How to get bam file, you can also check here.
--vcf_file: "S12878_freebayes.vcf" is VCF file generated from variant caller like "FreeBayes". How to get vcf file, you can also check here.
--uniq_map_dir: "Aquila/Uniqness_map" is the uniqness file for GRCh38 you can download by "./install.sh".
--out_dir, default = ./Asssembly_results. You can define your own folder, for example "Assembly_results_S12878".
--num_threads, default = 8. It's recommended not to change this setting unless large memory node could be used (2*memory capacity(it suggests for assembly below)), then try to use "--num_threads 12".
--chr_start, --chr_end: if you only want to assembly some chromosomes or only one chromosome. For example: use "--chr_start 1 --chr_end 5" will assemble chromsomes 1,2,3,4,5. Use "--chr_start 2 --chr_end 2" will only assemlby chromosome 2.
To use the above option "--chr_start, --chr_end", it is recommended to run the below command first to save more time later.
python Aquila/bin/Aquila_step0_sortbam.py --bam_file possorted_bam.bam --out_dir Assembly_results_S12878 --num_threads_for_bwa_mem 20
| Coverage | Memory | Time for chr1 on a single node |
|---|---|---|
| 60X | 100GB | 13:12:35 |
| 90X | 150GB | 1-01:08:38 |
| Coverage | Memory | Time for WGS on a single node |
|---|---|---|
| 60X | 450GB | 2-14:27:37 |
| 90X | 500GB | 3-12:34:12 |
Aquila/bin/Aquila_step2.py --out_dir Assembly_results_S12878 --num_threads 30 --reference Aquila/source/ref.fa
--out_dir, default = ./Asssembly_results, make sure it's the same as "--out_dir" from Step1 if you want to define your own output directory name.
--num_threads, default = 30, this determines the number of files assembled simultaneously by SPAdes.
--chr_start, --chr_end: if you only want to assembly some chromosomes or only one chromosome. For example: use "--chr_start 1 --chr_end 2"
| Coverage | Memory | Time for chr1 on a single node | --num_threads | --num_threads_spades |
|---|---|---|---|---|
| 60X | 100GB | 09:50:43 | 30 | 20 |
| 90X | 100GB | 13:25:08 | 40 | 20 |
| Coverage | Memory | Time for WGS on a single node | --num_threads | --num_threads_spades |
|---|---|---|---|---|
| 60X | 100GB | 3-12:16:27 | 40 | 20 |
| 90X | 100GB | 4-15:00:00 | 40 | 20 |
If your hard drive storage is limited, it is suggested to quily clean some data by running "Aquila_clean.py". Or you can keep them for some analysis.
Aquila/bin/Aquila_clean.py --out_dir Assembly_results_S12878
Aquila outputs an overall contig file “Aquila_Contig_chr*.fasta” for each chromosome, and one contig file for each haplotype: “Aquila_Contig_chr*_hp1.fasta” and “Aquila_Contig_chr*hp2.fasta”. For each contig, the header, for an instance, “>36_PS39049620:39149620_hp1” includes contig number “36”, phase block start coordinate “39049620”, phase block end coordinate “39149620”, and haplotype number “1”. Within the same phase block, the haplotype number “hp1” and “hp2” are arbitrary for maternal and paternal haplotypes. For some contigs from large phase blocks, the headers are much longer and complex, for an instance, “>56432_PS176969599:181582362_hp1 merge177969599:178064599_hp1-177869599:177969599_hp1”. “56” denotes contig number, “176969599” denotes the start coordinate of the final big phase block, “181582362” denotes the end coordinate of the final big phase block, and “hp1” denotes the haplotype “1”. “177969599:178064599_hp1” and “177869599:177969599_hp1” mean that this contig is concatenated from minicontigs in small chunk (start coordinate: 177969599, end coordinate: 178064599, and haplotype: 1) and small chunk (start coordinate: 177869599, end coordinate: 177969599, and haplotype: 1).
For example, you can use "Assemlby_results_S12878" as input directory to generate a VCF file which includes SNPs, small Indels and SVs, and the phased profile of all of them.
Please check check Assembly_based_variants_call_and_phasing for details.
1. Download hg19 reference from 10x Genomics website
wget http://xinzhouneuroscience.org/wp-content/uploads/2019/06/Uniqness_map_hg19.tar.gz
If you want to run Aquila for other diploid species with high quality reference genomes, to generate "Uniqness_map" for Aquila, check the details of hoffmanMappability to get the corresponding "k100.umap.bed.gz", then run "Aquila/bin/Get_uniqnessmap_for_Aquila.py" to get the final "Uniqness_map" folder to run Aquila.
Or you can use our "Aquila_uniqmap" to generate the "Uniqness_map" folder to run Aquila, check How_to_get_Umap for details.
Aquila/bin/Aquila_step1_multilibs.py --bam_file_list ./S24385_Lysis_2/Longranger_align_bam/S24385_lysis_2/outs/possorted_bam.bam,./S24385_Lysis_2H/Longranger_align_bam/S24385_lysis_2H/outs/possorted_bam.bam --vcf_file_list ./S24385_lysis_2/Freebayes_results/S24385_lysis_2_grch38_ref_freebayes.vcf,./S24385_lysis_2H/Freebayes_results/S24385_lysis_2H_grch38_ref_freebayes.vcf --sample_name_list S24385_lysis_2,S24385_lysis_2H --out_dir Assembly_results_merged --uniq_map_dir Aquila/Uniqness_map
--bam_file: "possorted_bam.bam" is bam file generated from barcode-awere aligner like "Lonranger align". Each bam file is seperately by comma (",").
--vcf_file: "S12878_freebayes.vcf" is VCF file generated from variant caller like "FreeBayes". Each VCF file is seperately by comma (",").
--sample_name: S24385_lysis_2,S24385_lysis_2H are the sample names you can define. Each sample name is seperately by comma (",").
--num_threads, default = 8. It's recommended not to change this setting unless large memory node could be used (2*memory capacity(it suggests for assembly below)), then try to use "--num_threads 12".
--chr_start, --chr_end: if you only want to assembly some chromosomes or only one chromosome. For example: use "--chr_start 1 --chr_end 5" will assemble chromsomes 1,2,3,4,5. Use "--chr_start 2 --chr_end 2" will only assemlby chromosome 2.
To use the above option "--chr_start, --chr_end", it is recommended to run the below command first to save more time later.
python Aquila/bin/Aquila_step0_sortbam_hybrid.py --bam_file_list ./S24385_Lysis_2/Longranger_align_bam/S24385_lysis_2/outs/possorted_bam.bam,./S24385_Lysis_2H/Longranger_align_bam/S24385_lysis_2H/outs/possorted_bam.bam --out_dir Assembly_results_merged --num_threads_for_bwa_mem 10 --sample_name_list S24385_lysis_2,S24385_lysis_2H
Aquila/bin/Aquila_step2.py --out_dir Assembly_results_merged --num_threads 30 --reference Aquila/source/ref.fa