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
For any queries, please either ask on github issue page or send an email to Jay Ghurye (jayg@cs.umd.edu).
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]
[-c CUTOFF] [-g GFA] [-u UNITIGS] [-t TENX] [-e ENZYME]
[-i ITER] [-x DUP] [-s EXP] [-m CLEAN]
SALSA Iterative Pipeline
optional arguments:
-h, --help show this help message and exit
-a ASSEMBLY, --assembly ASSEMBLY
Path to initial assembly
-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
-c CUTOFF, --cutoff CUTOFF
Minimum contig length to scaffold, default=1000
-g GFA, --gfa GFA GFA file for assembly
-u UNITIGS, --unitigs UNITIGS
The tiling of unitigs to contigs in bed format
-t TENX, --tenx TENX 10x links tab separated file, sorted by last columnls
-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
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 sored 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 which enzyme 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
.
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 -u unitigs_tiling.bed
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