New: Support added to SALSA for DNAse HiC Data. To run SALSA in this mode please use -e DNASE option.
New: Support added to SALSA to convert scaffolds to .hic format for visualization. You can load this file to Juicebox to visualize HiC contact map.
This code is used to scaffold your assemblies using Hi-C data. This version implements some improvements in the original SALSA algorithm. If you want to use the old version, it can be found in the old_salsa branch.
To use the latest version, first run the following commands:
cd SALSA
make
To run the code, you will need Python 2.7, BOOST libraries and Networkx(version lower than 1.2).
If you consider using this tool, please cite our publication which describes the methods used for scaffolding.
Ghurye, J., Pop, M., Koren, S., Bickhart, D., & Chin, C. S. (2017). Scaffolding of long read assemblies using long range contact information. BMC genomics, 18(1), 527. Link
Ghurye, J., Rhie, A., Walenz, B.P., Schmitt, A., Selvaraj, S., Pop, M., Phillippy, A.M. and Koren, S., 2018. Integrating Hi-C links with assembly graphs for chromosome-scale assembly. bioRxiv, p.261149 Link
For any queries, please either ask on github issue page or send an email to Jay Ghurye ([email protected]).
The new version of SALSA has been designed to consider several use cases depending on the input. Some assemblers output assembly graph as well along with the contig sequences. We provide options to use different information provided by the assembly to use for the scaffolding. Here is the what input options look like
python run_pipeline.py -h
usage: run_pipeline.py [-h] -a ASSEMBLY -l LENGTH -b BED [-o OUTPUT] [-O]
[-c CUTOFF] [-g GFA] [-e ENZYME] [-i ITER] [-x DUP]
[-s EXP] [-m CLEAN] [-f FILTER] [-p PRNT]
SALSA Iterative Pipeline
optional arguments:
-h, --help show this help message and exit
-a ASSEMBLY, --assembly ASSEMBLY
Path to initial assembly, headers must not contain ':'
-l LENGTH, --length LENGTH
Length of contigs at start
-b BED, --bed BED Bed file of alignments sorted by read names
-o OUTPUT, --output OUTPUT
Output directory to put results
-O, --output-original-coords
Run python run_pipeline.py -h to see the help
message for this option.
-c CUTOFF, --cutoff CUTOFF
Minimum contig length to scaffold, default=1000
-g GFA, --gfa GFA GFA file for assembly
-e ENZYME, --enzyme ENZYME
Restriction Enzyme used for experiment
-i ITER, --iter ITER Number of iterations to run, default = 3
-x DUP, --dup DUP File containing duplicated contig information
-s EXP, --exp EXP Expected Genome size of the assembled genome
-m CLEAN, --clean CLEAN
Set this option to "yes" if you want to find
misassemblies in input assembly
-f FILTER, --filter FILTER
Filter bed file for contigs present in the assembly
-p PRNT, --prnt PRNT Set this option to "yes" if you want to output the
scaffolds sequence and agp file for each iteration
To start the scaffolding, first step is to map reads to the assembly. We recommend using BWA or BOWTIE2 aligner to map reads. The read mapping generates a bam file. SALSA requires bed file as the input. This can be done using the bamToBed command from the Bedtools package. Also, SALSA requires bed file to be sorted by the read name, rather than the alignment coordinates. Once you have bam file, you can run following commands to get the bam
file needed as an input to SALSA.
Since Hi-C reads and alignments contain experimental artifacts, the alignments needs some postprocessing. To align and postprocess the alignments, you can use the pipeline released by Arima Genomics which can be found here https://github.com/ArimaGenomics.
bamToBed -i alignment.bam > alignment.bed
sort -k 4 alignment.bed > tmp && mv tmp alignment.bed
SALSA requires contig lengths as an input. You can generate this file using this command on your contig sequence file.
samtools faidx contigs.fasta
This will generate contigs.fasta.fai as an input for -l option. You can use this file for contig lengths while running SALSA.
Hi-C experiments can use different restriction enzymes. We use the enzyme frequency in contigs to normalize the Hi-C interaction frequency. You will need to specify the restriction site for the enzyme which was used for Hi-C experiment while running SALSA in -e option. If multiple enzymes were used, they can specified by separating with comma without space, like -e GATC,AAGCTT. Note that you need to specify the actual sequence of the cutting site for a restriction enzyme and not the enzyme name. For example, if you use MboI in the Hi-C protocol ,then you would specify it as -e GATC. You can also specify DNASE as an enzyme if you use an enzyme-free prep, e.g. Omin-C, as -e DNASE.
This is the minimum input you will require Suppose you only have contig sequences generated. Once you prepare the bed file as described above, the code can be run as follows:
python run_pipeline.py -a contigs.fasta -l contigs.fasta.fai -b alignment.bed -e {Your Enzyme} -o scaffolds
2) I have contig sequences and the alignment bam file but also want to use Hi-C data to correct input assembly errors
We also implemented a method in SALSA that can correct some of the errors in the assembly with Hi-C data. To use this method, you need to run following
python run_pipeline.py -a contigs.fasta -l contigs.fasta.fai -b alignment.bed -e {Your Enzyme} -o scaffolds -m yes
If you want to know what were the locations in the contigs where SALSA found errors, you can look at the input_breaks file in the output directory.
Some assembles output gfa file for the assembly graph. You can use that as an input for SALSA as follows
python run_pipeline.py -a contigs.fasta -l contigs.fasta.fai -b alignment.bed -e {Your Enzyme} -o scaffolds -m yes -g contigs_graph.gfa
We utilize graph to guide the scaffolding, which in turn reduces the errors.
Assemblers sometimes outputs unitigs along with contigs. Usually unitigs are shorter in size compared to contigs and hence contain much fewer errors. Because of this, the error correction with Hi-C data is not usually needed if unitigs are used as an input. Also, assembly graph contributes more to the scaffolding if used with unitigs instead of contigs. If you have all this data, you can run SALSA as follows.
python run_pipeline.py -a unitigs.fasta -l unitigs.fasta.fai -b alignment.bed -e {Your Enzyme} -o scaffolds -m yes -g unitigs_graph.gfa
Note here that using unitigs_tiling.bed is not mandatory to run SALSA in this mode. You can still run SALSA with -u option not set.
SALSA generates a bunch of files in the output folder. SALSA is an iterative algorithm, so it generates files for each iteration. The files you will be interested in are scaffolds_FINAL.fasta, which contains the sequences for the scaffolds generated by the algorithm. Another file which is of interest is scaffolds_FINAL.agp, which is the agp style output for the scaffolds describing the assignment, orientation and ordering of contigs along the scaffolds.
These are the postprocessing options which SALSA provides.
This can happen as unitigs are shorter chunks of sequence. If you have a bed file describing the tiling of unitigs along the contigs, you can use it to close the gaps inserted during the scaffolding. This script can be run as follows:
python stitch.py -c contigs.fasta -u unitigs.fasta -p scaffolds_iteration_x.p -o output_scaffolds.fasta
Here, scaffolds_iteration_x.p is the pickle file for scaffolds. If your code ran for 3 iterations, this file will be present in the output folder as scaffolds_iteration_3.p.
To dig more into the output, we have added an option that can output scaffolds along with the agp file for all intermediate iterations. This option is usually helpful in debugging and exploring the errors. Here is how you can run it:
python run_pipeline.py -a unitigs.fasta -l unitigs.fasta.fai -b alignment.bed -e {Your Enzyme} -o scaffolds -m yes -g unitigs_graph.gfa -u unitigs_tiling.bed -p yes
You need to have a jar file for Juicertools in order to generate .hic file. It can be found here. One you download this file on your machine, set the path to the jar file (including the jar file itself) in the convert.sh script to JUICER_JAR variable. Once that's done, simply run convert.sh script as convert.sh SALSA_OUT_DIR, where SALSA_OUT_DIR is the directory where the SALSA scaffolds and all the intermediate output resides. Please make sure you have latest version of GNU core utils in order to make use of --parallel option in GNU sort. After this script finishes execution, the .hic file will be generated in the SALSA_OUT_DIR as salsa_scaffolds.hic. This file can be loaded to Juicebox and visualized.