Replace mapping with de novo assembly

[donkirkby]

We currently use bowtie2 to map reads to a large set of reference sequences. For most samples, it works well. However, we have had some problems with reference drift (#290), calling HCV subtypes (#436), insertion and deletion positions (#398), and samples that produce different results when you rerun the mapping (#405).

We'd like to experiment with using de novo assembly instead of mapping.

• [x] use de novo assembly in micall pipeline
• [x] soft clipping is lost, possibly because of contig names
• [x] run one iteration of the remap step to map all the reads onto all the contigs
• [x] use Smith-Waterman Gotoh to align all the contigs onto the references
• [x] display contigs aligned with references
• [x] display contigs that extend beyond the references, with 1693P1Y04608-1E10-HIV_S22 as an example
• [x] display labels on short contigs by using an uncoloured track
• [x] display contigs grouped by reference, with 2569P1Y02945-1F6-HIV_S28 as an example
• [x] sort contigs by reference, then by descending coverage
• [x] add a dashed banner to the "Partial Blast Results" section, with dashes 100 bases long
• [x] smooth out coverage plots to reduce file sizes
• [x] show deletions as gaps in tracks
• [x] show insertions as green blocks
• [x] show partial segments in the banner above partial blast results, with 90314A-HCV_S49 as an example
• [x] use Smith-Waterman to align all the primers onto all the contigs moved to #478.
• [x] combine counts from all the contigs onto the coordinate references
• [x] V3LOOP doesn't build a contig, because the amplicon is shorter than 540 bases (HLA will probably have the same problem.)
• [x] dashes for partial blast results should begin at zero, not left border, with 61832A-HCV-43562A-GRrerun-HLA-B-POSA1211-5-2-V3LOOP_S2 as an example
• [x] some insertions show up to the left of the contig (not considering the offset for negative positions?), with 61832A-HCV-43562A-GRrerun-HLA-B-POSA1211-5-2-V3LOOP_S2 as an example
• [x] add HLA gene map
• [x] support mapping and assembly versions in a single Singularity image
• [x] make MiCall watcher call resistance after denovo steps
• [x] let MiCall watcher run denovo without standard pipeline
• [x] update columns in cascade.csv
• [x] use denovo pipeline as a backup for denovo combined pipeline in Kive watcher
• [x] download contig coverage files
• [x] add nuc_detail.csv
• [x] combine amino_detail.csv and nuc_detail.csv by seed groups, not by seeds. For example, sample 1693-1IN2C2-HIV_S16 from the 09-Aug-2019.M01841 run. Not needed after embedding contigs in ref.
• [x] check if assembly is stable - sample 1693-1IN2C2-HIV_S16 from the 09-Aug-2019.M01841 run had two very different sets of contigs. Difference between workstation and cluster? Basespace and Kive?
• [x] rename contig_coverage files to genome_coverage, and produce them from both the denovo version and the mapped version
• [x] collect genome_coverage.svg files in a separate folder from the other coverage maps
• [x] combine deletions, partial deletions, insertions, and low quality counts from the detail files
• [x] try aligning contigs in the reverse direction to see which aligns better, or look at blast result. HIV0836-P2K10-HIV_S26 in the 30-Aug-2019.M04401 run is an example.
• [x] embed contigs into a complete reference, to try and map regions that didn't assemble
• [x] lots of samples have reversed sections near the ends - report reversed contigs with partials
• [x] de novo assembly is very slow for some samples IVA seems better than savage.
• [x] bring back G2P analysis
• [x] should G2P continue to use merged reads, or should it switch to aligned reads? moved to #481
• [x] should we try to report V3LOOP overlap again? moved to #481
• [x] remove the built index of genotypes from source control
• [x] add other viruses to index
• [x] don't show separate contigs for -partial or -reversed suffixes.
• [x] show deletions in contigs - 1497DEL-HIV_S119 from 7 Dec 2018 is an example
• [x] make microtest samples work with de novo assembly
• [x] add tests for microtest samples, or an option to verify their results in micall_basespace.py
• [x] should we make the contig coverage diagram match what we used to cut up the gene regions? 73051ANS5A1-HCV-NS5a_S89 from 15-Jul-2016.M01841 is an example where they don't match. moved to #479.
• [x] denovo assembly sometimes uses more than one thread, specifically smalt uses up to 8. Either reduce CPU usage, or increase number of threads requested from Slurm.
• [x] use BLAST results to assemble contigs into a full reference? Haven't found any clear cases where they should be combined. HIV3428P100IN200-C19-HIV-S51 from 20 Sep 2019 run is the closest, but it looks like one contig has primer at the end. Samples HIV0887-P2D21-HIV_S3 and HIV0887-P2C12-HIV_S32 from 30 Aug 2019 looks even better, but have very little overlap. Some of the HCV samples look more promising: 73060A-HCV_S46 from 15 Jul 2016, for example. Moved to issue #484.
• [x] display BLAST results on genome coverage diagram: find the top scoring result, then display all others on the same reference, numbering each section
• [x] put arrow on top of BLAST results in reference landmarks
• [x] adjust BLAST arrows when diagram has a horizontal offset, as in JRCCC3-GP160NEF-HIV_S83 from 27 Sep 2019
• [x] bring back amino_details.csv and combination into amino.csv
• [x] look at BLAST results for recombinant samples from 27 Sep 2019 - are there any where the best blast hit for one region matches a different reference from the best last hit for another region? HCV example: 73087A-HCV_S2 from 15 Jul 2016 has a gap in BLAST matches from C to NS2 in HCV-1b, but has a big match to HCV-2k from 5' to NS2, and 73088A-HCV_S15 is similar.
• [x] Reference arrows aren't aligned right vertically (B5-VIR5POL-HIV_S12) or horizontally (B4-VIR4POL-HIV_S8), examples are from 27 Sep 2019
• [x] compare assembled and remapped versions with previous release
• [x] deal with HIV references that don't reach 5' and 3' ends or bring back the refs we removed Moved to issue #484.
• [x] check for similar problems with other seed groups Moved to issue #484.
• [x] remove Savage?
• [x] remove second version of consensus, with lower minimum coverage?

[donkirkby]

@jeff-k has been using de novo assembly to look at several samples that got strange results with the current MiCall pipeline. It looks like one advantage of the technique will be that we can distinguish between these two scenarios:

• high error rates (generates many, unrelated contigs)
• unknown sequence like a chimera or new strain (generates a small number of related contigs, but the contigs don't map well to known references)

With the current MiCall pipeline, both of those scenarios just look like lousy mapping with gaps in coverage.

We propose this plan for a full MiCall pipeline that includes de novo assembly:

• [x] trim adapter sequences
• [x] use de novo assembly to generate a set of contigs
• [x] combine the contigs into one or more consensus sequences using prelim_map and remap as we combine the reads today
• [x] do a final remap with all of the reads against the consensus sequences
• [x] report the reads aligned against the consensus sequences
• [x] count the coverage and report the mixtures at each position
• [x] generate coverage maps
• [x] generate resistance interpretations

As you can see, this just affects the prelim_map and remap steps. Because they're only running on a small number of contigs, they should be much faster.

One risk is that the de novo assembly step might be much slower in some cases than the remap step is currently.