bcgsc / mavis

Merging, Annotation, Validation, and Illustration of Structural variants

About

MAVIS is python command-line tool for the post-processing of structural variant calls. The general MAVIS pipeline consists of six main stages

• convert
• cluster
• validate
• annotate
• pairing
• summary

Getting Help

All steps in the MAVIS pipeline are called following the main mavis entry point. The usage menu can be viewed by running without any arguments, or by giving the -h/--help option

mavis -h

Help sub-menus can be found by giving the pipeline step followed by no arguments or the -h options

mavis cluster -h

Common problems and questions are addressed on the wiki. If you have a question or issue that is not answered there (or already a github issue) please submit a github issue to our github page or contact us by email at [email protected]

Getting Started

The simplest way to use MAVIS is via Singularity. The MAVIS docker container used by singularity will take care of installing the aligner as well.

pip install -U setuptools pip wheel
pip install mavis_config  # also installs snakemake

Now you will run mavis via Snakemake as follows

snakemake \
    -j <MAX JOBS> \
    --configfile <YOUR CONFIG> \
    --use-singularity \
    -s Snakefile

For other installation options which do not use docker/singularity see the comprehensive install instructions in the user manual

Citation

If you use MAVIS as a part of your project please cite

Reisle,C. et al. (2018) MAVIS: Merging, Annotation, Validation, and Illustration of Structural variants. Bioinformatics.

mavis's Issues

Split Flanking read support to multiple calls

Sometimes the flanking evidence collected may support multiple events. Currently the evidence is not used if the coverage region (region that the aligned portion flanking reads span at a single breakpoint) is larger than expected.

Instead it should be split into multiple groups of evidence

Drop Breakpoint pairs (with warning) outside chromosomes

For breakpoints which are called passed the end or overlapping past the end of a chromosome we expect errors. Would like to handle these and allow mavis to run past the bad input. Should output a warning documenting what was dropped. Currently the error looks like

[2018-01-08 11:28:30.307058] loaded 439 breakpoint pairs
Traceback (most recent call last):
  File "/gsc/pipelines/mavis/v1.5.1/venv/bin/mavis", line 11, in <module>
    load_entry_point('mavis==1.5.1', 'console_scripts', 'mavis')()
  File "/gsc/pipelines/mavis/v1.5.1/venv/lib/python3.6/site-packages/mavis-1.5.1-py3.6.egg/mavis/main.py", line 622, in main
  File "/gsc/pipelines/mavis/v1.5.1/venv/lib/python3.6/site-packages/mavis-1.5.1-py3.6.egg/mavis/validate/main.py", line 101, in main
  File "/gsc/pipelines/mavis/v1.5.1/venv/lib/python3.6/site-packages/mavis-1.5.1-py3.6.egg/mavis/validate/evidence.py", line 15, in __init__
  File "/gsc/pipelines/mavis/v1.5.1/venv/lib/python3.6/site-packages/mavis-1.5.1-py3.6.egg/mavis/validate/base.py", line 64, in __init__
ValueError: Breakpoint 81194267-81195629 is outside the range of the reference sequence 17 (1-81195210)

Output masked pairs to filtered file

To make the tracking easier of which breakpoint pairs were run through clustering and subsequent stages, all filtered pairs should be output to a file which lists the filtering reason. This will include:

• pairs filtered on non-standard chromosomes
• pairs in other libraries
• pairs overlapping masking regions
• pairs considered uinformative

Suppress output of intermediary files

To reduce the disk space usage we will add options to not output (or clean up) intermediary files that are not required for subsequent pipeline steps. This will include the evidence bams, contig alignment bams, and aligner output files and logs.

braces in config file for v1.6.4 is not compatible with 1.6.6

The brace expansion in 1.6.4 was moved in 1.6.6 and no longer expands braces for arguments in the config file as if the file were generated with 1.6.6 they would already be expanded

Increase caching strict filtering

• cache less reads. Currently any read that indicates an alteration will be cached
• can focus this more on the current event as well as ignoring intron spans for transcriptomes

Add ymax and strand support for overlay diagrams

• currently the overlay function is just using the builtin samtools pileup for the read depth plots.
• want to create own pielup to support filtering reads based on mapping quality, strand, etc.

Add hostname logging

Add logging of the hostname to the mavis output for each stage of the pipeline. This will help in debugging on the cluster should a particular node have an issue

translocations filtered out with homopolymer RNA filter

MAVIS version: v1.6.0

the filter for trans homopolymers does not check event type before filtering, since all translocations have an event size of 0, some are mistakenly filtered out

Drawing Error for single bp insertions

mavis version: 1.6.8

error being reported when mavis attempts to draw single bp insertions

Traceback (most recent call last):
  File "/home/creisle/git/mavis/venv/bin/mavis", line 11, in <module>
    load_entry_point('mavis', 'console_scripts', 'mavis')()
  File "/home/creisle/git/mavis/mavis/main.py", line 634, in main
    annotate_main(**args, start_time=start_time)
  File "/home/creisle/git/mavis/mavis/annotate/main.py", line 245, in main
    drawing, legend = draw(drawing_config, ann, reference_genome, template_metadata, drawings_directory)
  File "/home/creisle/git/mavis/mavis/annotate/main.py", line 51, in draw
    **other_settings
  File "/home/creisle/git/mavis/mavis/illustrate/diagram.py", line 292, in draw_sv_summary_diagram
    reference_genome=reference_genome
  File "/home/creisle/git/mavis/mavis/illustrate/elements.py", line 385, in draw_ustranscript
    config, canvas, tl, labels, colors, mapping, genomic_min=genomic_min, genomic_max=genomic_max
  File "/home/creisle/git/mavis/mavis/illustrate/elements.py", line 136, in draw_transcript_with_translation
    genomic_max=genomic_max
  File "/home/creisle/git/mavis/mavis/illustrate/elements.py", line 95, in draw_exon_track
    pxi = Interval(start, end)
  File "/home/creisle/git/mavis/mavis/interval.py", line 28, in __init__
    raise AttributeError('interval start > end is not allowed', self.start, self.end)
AttributeError: ('interval start > end is not allowed', 1375.4481108961684, 1351.4400384766157)

This is because of using the 'skinny' exons where the inside edge from each pixel range is taken. when the start/end are equal this causes the start to be > end

max_proximity not propogated through annotation stage

The max_proximity setting is not being passed down to calls that the main function of the annotation stage makes and therefore the default is being used.

Gene2 is colored wrong in gene level

When a gene is both overlapping the opposite breakpoint and the gene selected at the current breakpoint, if the current breakpoint is the second breakpoint the colors are incorrect.

Reduce density in overlay read depth plots

There is a lot of data in the read depth plots, and due to this the svgs can be very large. However b/c the diagram is limited to a certain size, a lot of this data overlaps and could be removed without losing anything. This is more effective than the previous binning argument to reduce the file size.

Change blat rank to aligner rank

Calculate the 'aligner' rank for bwa and blat. For blat this will be based on scoring and for bwa this will be based on primary/supplementary as well as mapping qualities

Then this can be used to filter out alignments of the same contig

Call duplications from insertions

Logic mostly in place just not being used. Should compare the untemplated sequence to the reference genome during the contig alignment selection and create a duplication split alignment where possible

Checker reporting parallel run time as sum of libraries vs max

MAVIS version: 1.5.2

Python version: 3.6.1

OS: centos06

Expected Behaviour

The parallel run time summary should be the estimated max when all things are run parallel

Actual Behaviour

The max of each library is summed. Basically assumes that the libraries are not parallel

Steps to Reproduce the Behaviour

run

mavis checker -o /path/to/any/multi/library/run

using inner boundaries for exon pixel mapping

When the exons are drawn, currently the largest possible exon is drawn and sometimes they overlap. It would be better to use the smallest possible values. This variability happens because sometime 1 genomic position will map to a range of pixels.

Improve error message on breakpoint outside chr length

The error message when the input breakpoint is outside the length of the chromosome is not very helpful. It throws an error message when validating split reads. It would be more useful if the error were thrown when the Evidence object is created

Current error message

[2017-12-14 06:48:24.463580] (1 of 31) gathered evidence for: HgrjuHPaiN5LGqt9txW3Ya (tracking_id: hZWpx6FrgAGzpAdRFqwzWZ)
                             BPP(Breakpoint(1:7660671-7666671L), Breakpoint(22:56571650-56577650R), opposing=False)
                             possible event type(s): {'translocation'}
                             outer window regions:  1:7659938-7666830  22:56571491-56578383
                             inner window regions:  1:7660512-7666830  22:56571491-56577809
Traceback (most recent call last):
  File "/projects/trans_scratch/software/mavis/v1.4.1/venv/bin/mavis", line 11, in <module>
    load_entry_point('mavis==1.4.1', 'console_scripts', 'mavis')()
  File "/projects/trans_scratch/software/mavis/v1.4.1/venv/lib/python3.6/site-packages/mavis-1.4.1-py3.6.egg/mavis/main.py", line 622, in main
  File "/projects/trans_scratch/software/mavis/v1.4.1/venv/lib/python3.6/site-packages/mavis-1.4.1-py3.6.egg/mavis/validate/main.py", line 149, in main
  File "/projects/trans_scratch/software/mavis/v1.4.1/venv/lib/python3.6/site-packages/mavis-1.4.1-py3.6.egg/mavis/validate/base.py", line 783, in load_evidence
  File "/projects/trans_scratch/software/mavis/v1.4.1/venv/lib/python3.6/site-packages/mavis-1.4.1-py3.6.egg/mavis/validate/base.py", line 506, in collect_split_read
  File "/projects/trans_scratch/software/mavis/v1.4.1/venv/lib/python3.6/site-packages/mavis-1.4.1-py3.6.egg/mavis/bam/read.py", line 249, in nsb_align
AttributeError: cannot overlap on an empty sequence

add general vcf convert option

Most of the SV tools which produce vcfs use similar notation, should allow users to upload as vcf type tool vs having to specify the exact tool.

Retain all possible data on grouping calls

For a group of calls merged during summary columns with differing values may not be retain all values after the merge. Currently only retaining specific columns that become a list.

For all other columns we should retain the value when all columns agree

Limit Contig sequences to investigate

• increase unique filtering on contig sequences. Reduce similar contigs being both passed
• Do not send contig sequences for alignment that do not have remapped reads represented from both side/breakpoints of the event

Mark supplementary events

Events calls as part of another event should be marked supplementary. This is b/c a germline event may be found as part of looking for a somatic event. Since the original event is not in the germline we would not expect to find evidence in the germline. Therefore if the supplementary event is present in the germline it may be incorrectly labelled as somatic during the summary stage

Add job count to checker summary

Add summary of the number of jobs run to the checker to make summarizing the run time trends easier

Add frameshift notation to protein hgvs notation

Currently calling frameshifting indels as indels deleting the whole tail of the protein and inserting the new sequence, while not technically incorrect this would be better represented as a frameshift

Resolve untemplated sequence for calls by split reads

Currently no attempt to resolve the untemplated sequence is made when calls are called by split read evidence. While this won't be possible in all cases it may be resolve in cases where 1 or more read has aligned to both breakpoints.

Calculate the breakpoints using all 'double'/target aligned split reads. If the untemplated sequence and any minor breakpoint shifting agrees then adjust the call accordingly

touble finding samtools version

MAVIS version: <1.6.10>

Python version: <3.6.1>

OS: <CentOS 6>

Expected Behaviour

SAMTools check success / no error

Actual Behaviour

I'm finding reports like the following in the output logs in the validate folders

Traceback (most recent call last):
File "/gsc/pipelines/mavis/v1.6.10/venv/bin/mavis", line 11, in
load_entry_point('mavis==1.6.10', 'console_scripts', 'mavis')()
File "/gsc/pipelines/mavis/v1.6.10/venv/lib/python3.6/site-packages/mavis-1.6.10-py3.6.egg/mavis/main.py", line 522, in main
File "/gsc/pipelines/mavis/v1.6.10/venv/lib/python3.6/site-packages/mavis-1.6.10-py3.6.egg/mavis/bam/read.py", line 169, in get_samtools_version
ValueError: unable to parse samtools version number

Steps to Reproduce the Behaviour

ssh numbers
cd /projects/rcorbettprj2/POGHOMD
source /gsc/pipelines/mavis/production/venv/bin/activate
source /gsc/pipelines/mavis/environments/hg19.sh
mavis config -w config.cfg --library A36971 genome diseased False /projects/POG/POG_data/TESTCOLO829STAR/wgs/biop1_t_A36971/merge/hg19a_bwa-mem-0.7.6a/A36971_2_lanes_dupsFlagged.bam --library A36973 genome normal False /projects/POG/POG_data/TESTCOLO829STAR/wgs/blood1_n_A36973/merge/hg19a_bwa-mem-0.7.6a/A36973_1_lane_dupsFlagged.bam --assign A36971 testHOMD.txt --assign A36973 testHOMD.txt 

before and after I submit the jobs to the cluster I get:

which samtools
/projects/trans_scratch/transabyss/trans-ABySS/v1.4.10/bin/samtools

and running samtools without parameters seems to generate the require output for get_samtools_version()

samtools

Program: samtools (Tools for alignments in the SAM format)
Version: 1.2 (using htslib 1.2.1)

Usage:   samtools <command> [options]

Commands:
  -- indexing
         faidx       index/extract FASTA
         index       index alignment
  -- editing

I can't help but wonder if a different samtools is getting run underneath. What do you recommend?

thanks
Richard

error on single bp deletion call by flanking pairs

MAVIS version: 1.6.3

Python version: 3.6

OS: centos06

Error not reported by checker and value error raised on calling single bp deletion by flanking pairs.
Should be ignoring the error but should also be reporting the error if raised in the checker output

calls resolve by split reads should remove strand if the evidence type is unstranded

Was not resetting the strand to NS for unstranded calls. This was resulting in multiple calls being retained and the resolution not being used.

Will also add support for multiple split read resolutions

Use subparsers for pipeline commands

Currently the pipeline commands are manually separated rather than using subparsers, would be better practice to use the built-in subparsers

Flag events in homopolymer regions

RNA reads often produce false positive calls in homopolymer regions due to slippage. To facilitate filtering these later on they should be flagged.

Example deletion

AATAGATCGGATATTTTTGGATCAGGATAGC
AATAGATCGGATATTTT-GGATCAGGATAGC

Move tab (private package) dependency into mavis repository

Currently mavis is dependent on a package that has not been made publicly available. As a stop gap solution the package will be added to the mavis repository

Drop unmapped Alignments

When using bwa as the aligner, sometimes unmapped contigs are returned. These should be dropped from the results as otherwise they will throw key/value errors when retrieving the reference sequence and they do not provide any useful information.

Retain shorted contigs

Filtering on contig length is too strict. Will lower requirements to allow shorted contigs through. Will limit to read length instead of 1.25 * read length

Error loading the annotations for generating a config for transcriptome libraries

MAVIS version: 1.6.7

Python version: 3.6

OS: centos06

Expected Behaviour

Run without errors

Actual Behaviour

                             annotations = None
                             annotations_filename = '/gsc/pipelines/mavis/reference_files/ensembl69_hg19_annotations.json'
...
usage: mavis [-h] [-v]
             {annotate,validate,pipeline,cluster,pairing,summary,checker,config,convert,overlay}
             ...
mavis: error: argument --annotations: is required to gather bam stats for transcriptome libr

Steps to Reproduce the Behaviour

Run the config command with any transcriptome library

output additional summary files

The default summary file is quite large, it would be nice to have other subsets of it output for easier viewing.
To start, one for each library containing only non-synonymous non-supplementary variants which produce a protein

Build fusion products for unresolved untemplated sequence

Currently only building fusion products when the breakpoints are exact and the untemplated sequence has been specified.
Instead, will assume the untemplated sequence is an empty string and proceed adding a column to note the assumption (still only for exact breakpoint calls)

summary dropping single split read calls

problem

Split reads calls are being dropped during the 'group_by_distance' stage when they are single nodes and not paired with any other similar or exact split read calls because the algorithm is creating a sparse adjacency matrix instead of a dense one.

solution

initialize the adjacency matrix

Throw error when config does not exist

The config is generated with no libraries when it does not exist. It would be more helpful if it threw an error instead

Add 'command not found' to checker errors to detect

Add "command not found" to the error parser to ensure that these errors are reported by the checker

annotation error in 1.6.9

MAVIS version: 1.6.9

Python version: 3.6

OS: centos06

A bunch of annotation jobs have died due to this index error:

[2018-02-27 17:46:11.421346] (22 of 38) current annotation ZV7do3tXi5AADpQqgdGwWW-v1-a1 PreTranscript(ENSG00000144642:29322803-30051886, name=ENST00000396583) PreTranscript(ENSG00000144642:29322803-30051886, name=ENST00000396583) duplication
                             drawing attempt: 1 1000px 
[2018-02-27 17:46:11.948482] closing: /projects/trans_scratch/mavis_validation/6_pairing_analysis/POG217/v1.6.9/P00637_normal_genome/annotate/batch-GpJraczAd2efbUpAvToauT-1/annotations.tab
[2018-02-27 17:46:11.951993] closing: /projects/trans_scratch/mavis_validation/6_pairing_analysis/POG217/v1.6.9/P00637_normal_genome/annotate/batch-GpJraczAd2efbUpAvToauT-1/annotations.fusion-cdna.fa
Traceback (most recent call last):
  File "/gsc/pipelines/mavis/v1.6.9/venv/bin/mavis", line 11, in <module>
    load_entry_point('mavis==1.6.9', 'console_scripts', 'mavis')()
  File "/gsc/pipelines/mavis/v1.6.9/venv/lib/python3.6/site-packages/mavis-1.6.9-py3.6.egg/mavis/main.py", line 634, in main
  File "/gsc/pipelines/mavis/v1.6.9/venv/lib/python3.6/site-packages/mavis-1.6.9-py3.6.egg/mavis/annotate/main.py", line 245, in main
  File "/gsc/pipelines/mavis/v1.6.9/venv/lib/python3.6/site-packages/mavis-1.6.9-py3.6.egg/mavis/annotate/main.py", line 51, in draw
  File "/gsc/pipelines/mavis/v1.6.9/venv/lib/python3.6/site-packages/mavis-1.6.9-py3.6.egg/mavis/illustrate/diagram.py", line 292, in draw_sv_summary_diagram
  File "/gsc/pipelines/mavis/v1.6.9/venv/lib/python3.6/site-packages/mavis-1.6.9-py3.6.egg/mavis/illustrate/elements.py", line 395, in draw_ustranscript
  File "/gsc/pipelines/mavis/v1.6.9/venv/lib/python3.6/site-packages/mavis-1.6.9-py3.6.egg/mavis/illustrate/elements.py", line 146, in draw_transcript_with_translation
  File "/gsc/pipelines/mavis/v1.6.9/venv/lib/python3.6/site-packages/mavis-1.6.9-py3.6.egg/mavis/illustrate/elements.py", line 78, in draw_exon_track
  File "/gsc/pipelines/mavis/v1.6.9/venv/lib/python3.6/site-packages/mavis-1.6.9-py3.6.egg/mavis/interval.py", line 353, in convert_ratioed_pos
IndexError: (1, 'is outside mapped range', {Interval(8, 446): Interval(0.0, 45.035099575560835, type=<class 'float'>), Interval(447, 153432): Interval(45.035099575560835, 59.06732425775306, type=<class 'float'>), Interval(153433, 153605): Interval(59.06732425775306, 91.35914951759821, type=<class 'float'>), Interval(153606, 207130): Interval(91.35914951759821, 102.76989838425489, type=<class 'float'>), Interval(207131, 207189): Interval(102.76989838425489, 129.60032036593617, type=<class 'float'>), Interval(207190, 305803): Interval(129.60032036593617, 142.19947195372657, type=<class 'float'>), Interval(305804, 305895): Interval(142.19947195372657, 170.6108264632974, type=<class 'float'>), Interval(305896, 458409): Interval(170.6108264632974, 184.63061072577804, type=<class 'float'>), Interval(458410, 458567): Interval(184.63061072577804, 216.20383029112796, type=<class 'float'>), Interval(458568, 481599): Interval(216.20383029112796, 226.81088060213844, type=<class 'float'>), Interval(481600, 481679): Interval(226.81088060213844, 254.64735055611308, type=<class 'float'>), Interval(481680, 587547): Interval(254.64735055611308, 267.4376945265047, type=<class 'float'>), Interval(587548, 587654): Interval(267.4376945265047, 296.5676547305708, type=<class 'float'>), Interval(587655, 599468): Interval(296.5676547305708, 306.87903421886006, type=<class 'float'>), Interval(599469, 599507): Interval(306.87903421886006, 332.7513152745477, type=<class 'float'>), Interval(599508, 602851): Interval(332.7513152745477, 342.83945248538566, type=<class 'float'>), Interval(602852, 602898): Interval(342.83945248538566, 369.09498991147075, type=<class 'float'>), Interval(602899, 616068): Interval(369.09498991147075, 379.44210924960936, type=<class 'float'>), Interval(616069, 616165): Interval(379.44210924960936, 408.0929989906786, type=<class 'float'>), Interval(616166, 618394): Interval(408.0929989906786, 418.1517483456301, type=<class 'float'>), Interval(618395, 618445): Interval(418.1517483456301, 444.5989139569139, type=<class 'float'>), Interval(618446, 654775): Interval(444.5989139569139, 455.5564572792051, type=<class 'float'>), Interval(654776, 654886): Interval(455.5564572792051, 484.87804566846995, type=<class 'float'>), Interval(654887, 706832): Interval(484.87804566846995, 496.2471771141, type=<class 'float'>), Interval(706833, 706913): Interval(496.2471771141, 524.1315541143742, type=<class 'float'>), Interval(706914, 709771): Interval(524.1315541143742, 534.2068819100007, type=<class 'float'>), Interval(709772, 709899): Interval(534.2068819100007, 564.3428900863601, type=<class 'float'>), Interval(709900, 722545): Interval(564.3428900863601, 574.6761984500731, type=<class 'float'>), Interval(722546, 729085): Interval(574.6761984500731, 911.9921875, type=<class 'float'>)})

An example log is here: /projects/trans_scratch/mavis_validation/6_pairing_analysis/POG217/v1.6.9/P00637_normal_genome/annotate/batch-GpJraczAd2efbUpAvToauT-1/slurm-MA_P00637_batch-GpJraczAd2efbUpAvToauT-1-823770.log

speed up pairing algorithm

Improve pairing algorithm speed.

Current naive algorithm can be improved to not perform unnecessary comparisons

Will use the same method used for speeding up clustering by sorting the list twice by break1 and break2 and then looking locally only (neighbors in the list) until distance is exceeded.

Add BreakDancer support

Add conversion method for breakdancer (https://github.com/genome/breakdancer) to the supported tools

Simplify validate options

• remove unnecessary options
• group options that are always the same

typo in user manual

Lists the environment variable for the --dgv_annotation option as MAVIS_DGV_ANNOTATIONS when it should not have the s and should be MAVIS_DGV_ANNOTATION

Improve checker python error detection

checker is ignoring when the last line of the log file contains the following error type

/.../mavis: line 51: 128612 Aborted                 (core dumped) mavis "$@"

Will need to add 'aborted' to the list of words to look for in order to capture this

Add transcriptome-specific overrides for mapping quality and read limit settings

Currently there is no good way to set different default values for transcriptomes and genomes for validation settings. Will add transcriptome-specific override settings for min_mapping_quality and fetch_reads_limit validation settings since these are of particular importance to transcriptomes.

Account for untemplated sequence in pairing distance

Increase the maximum allowed distance for pairing by breakpoint when there is untemplated sequence.

Stack reference transcripts in fusion diagrams

Currently when there are too many exons to fit within the specified width, mavis defaults to a gene-level only diagram. This is considerably less useful than the full diagram, but if the width is increased to much the diagram becomes very wide and hard to read

Instead, an option should be added to stack the reference transcripts so that you don't need to fit two across

Clustering bug

Grouping events which we would not expect to be grouped.

{('4', '4', 157002046, 157002358, 157002046, 157002358, 'L', 'R', '?', '?', False, False), ('4', '4', 156935061, 156935245, 156935061, 156935245, 'L', 'R', '?', '?', False, False)}

breakpoint pair keys were bring grouped

Add mail option for slurm/sge job

Add mail options for slurm and sge job submission scripts