Merge pull request #128 from jbelyeu/master

max coverage limit

Author: ryanlayer

README.md

@@ -16,27 +16,8 @@ substantiate the SV.
   <summary>samplot plot</summary>
 
   ```
-usage: samplot plot [-h] [--marker_size MARKER_SIZE] [-n TITLES [TITLES ...]]
-                    [-r REFERENCE] [-z Z] -b BAMS [BAMS ...] -o OUTPUT_FILE -s
-                    START -e END -c CHROM [-w WINDOW] [-d MAX_DEPTH]
-                    [--minq MINQ] [-t SV_TYPE] [-T TRANSCRIPT_FILE]
-                    [-A ANNOTATION_FILES [ANNOTATION_FILES ...]]
-                    [--coverage_tracktype {stack,superimpose}] [-a]
-                    [-H PLOT_HEIGHT] [-W PLOT_WIDTH] [-q MIN_MQUAL] [-j]
-                    [--start_ci START_CI] [--end_ci END_CI]
-                    [--long_read LONG_READ] [--min_event_size MIN_EVENT_SIZE]
-                    [--xaxis_label_fontsize XAXIS_LABEL_FONTSIZE]
-                    [--yaxis_label_fontsize YAXIS_LABEL_FONTSIZE]
-                    [--legend_fontsize LEGEND_FONTSIZE]
-                    [--annotation_fontsize ANNOTATION_FONTSIZE]
-                    [--common_insert_size] [--hide_annotation_labels]
-                    [--coverage_only] [--same_yaxis_scales] [--zoom ZOOM]
-                    [--debug DEBUG]
-
 optional arguments:
   -h, --help            show this help message and exit
-  --marker_size MARKER_SIZE
-                        Size of marks on pairs and splits (default 3)
   -n TITLES [TITLES ...], --titles TITLES [TITLES ...]
                         Space-delimited list of plot titles. Use quote marks
                         to include spaces (i.e. "plot 1" "plot 2")
@@ -46,36 +27,53 @@ optional arguments:
   -b BAMS [BAMS ...], --bams BAMS [BAMS ...]
                         Space-delimited list of BAM/CRAM file names
   -o OUTPUT_FILE, --output_file OUTPUT_FILE
-                        Output file name
+                        Output file name/type. Defaults to
+                        {type}_{chrom}_{start}_{end}.png
+  --output_dir OUTPUT_DIR
+                        Output directory name. Defaults to working dir.
+                        Ignored if --output_file is set
   -s START, --start START
-                        Start position of region/variant
-  -e END, --end END     End position of region/variant
+                        Start position of region/variant (add multiple for
+                        translocation/BND events)
+  -e END, --end END     End position of region/variant (add multiple for
+                        translocation/BND events)
   -c CHROM, --chrom CHROM
-                        Chromosome
+                        Chromosome (add multiple for translocation/BND events)
   -w WINDOW, --window WINDOW
                         Window size (count of bases to include in view),
                         default(0.5 * len)
   -d MAX_DEPTH, --max_depth MAX_DEPTH
                         Max number of normal pairs to plot
-  --minq MINQ           coverage from reads with MAPQ <= minq plotted in
-                        lighter grey. To disable, pass in negative value
   -t SV_TYPE, --sv_type SV_TYPE
                         SV type. If omitted, plot is created without variant
                         bar
   -T TRANSCRIPT_FILE, --transcript_file TRANSCRIPT_FILE
-                        GFF of transcripts
+                        GFF3 of transcripts
+  --transcript_filename TRANSCRIPT_FILENAME
+                        Name for transcript track
+  --max_coverage_points MAX_COVERAGE_POINTS
+                        number of points to plot in coverage axis (downsampled
+                        from region size for speed)
   -A ANNOTATION_FILES [ANNOTATION_FILES ...], --annotation_files ANNOTATION_FILES [ANNOTATION_FILES ...]
                         Space-delimited list of bed.gz tabixed files of
                         annotations (such as repeats, mappability, etc.)
-  --coverage_tracktype {stack,superimpose}
+  --annotation_filenames ANNOTATION_FILENAMES [ANNOTATION_FILENAMES ...]
+                        Space-delimited list of names for the tracks in
+                        --annotation_files
+  --coverage_tracktype {stack,superimpose,none}
                         type of track to use for low MAPQ coverage plot.
   -a, --print_args      Print commandline arguments
   -H PLOT_HEIGHT, --plot_height PLOT_HEIGHT
                         Plot height
   -W PLOT_WIDTH, --plot_width PLOT_WIDTH
                         Plot width
-  -q MIN_MQUAL, --min_mqual MIN_MQUAL
+  -q INCLUDE_MQUAL, --include_mqual INCLUDE_MQUAL
                         Min mapping quality of reads to be included in plot
+                        (default 1)
+  --separate_mqual SEPARATE_MQUAL
+                        coverage from reads with MAPQ <= separate_mqual
+                        plotted in lighter grey. To disable, pass in negative
+                        value
   -j, --json_only       Create only the json file, not the image plot
   --start_ci START_CI   confidence intervals of SV first breakpoint (distance
                         from the breakpoint). Must be a comma-separated pair
@@ -86,9 +84,10 @@ optional arguments:
   --long_read LONG_READ
                         Min length of a read to be treated as a long-read
                         (default 1000)
+  --ignore_hp           Choose to ignore HP tag in alignment files
   --min_event_size MIN_EVENT_SIZE
                         Min size of an event in long-read CIGAR to include
-                        (default 100)
+                        (default 20)
   --xaxis_label_fontsize XAXIS_LABEL_FONTSIZE
                         Font size for X-axis labels (default 6)
   --yaxis_label_fontsize YAXIS_LABEL_FONTSIZE
@@ -100,10 +99,17 @@ optional arguments:
   --common_insert_size  Set common insert size for all plots
   --hide_annotation_labels
                         Hide the label (fourth column text) from annotation
-                        files, useful for region with many annotations
+                        files, useful for regions with many annotations
   --coverage_only       Hide all reads and show only coverage
+  --max_coverage MAX_COVERAGE
+                        apply a maximum coverage cutoff. Unlimited by default
   --same_yaxis_scales   Set the scales of the Y axes to the max of all
-  --zoom ZOOM           Only show +- zoom amount around breakpoints
+  --marker_size MARKER_SIZE
+                        Size of marks on pairs and splits (default 3)
+  --dpi DPI             Dots per inches (pixel count, default 300)
+  --zoom ZOOM           Only show +- zoom amount around breakpoints, much
+                        faster for large regions. Ignored if region smaller
+                        than --zoom (default 500000)
   --debug DEBUG         Print debug statements
 ```
 </details>

samplot/__init__.py

@@ -1,2 +1,2 @@
 #!/usr/bin/env python
-__version__ = "1.1.1"
+__version__ = "1.1.4"

samplot/samplot.py

@@ -177,7 +177,7 @@ def points_in_window(points):
 
 # {{{ def get_tabix_iter(chrm, start, end, datafile):
 def get_tabix_iter(chrm, start, end, datafile):
-    """Gets an iterator from a tabix BED/GFF/GFF3 file
+    """Gets an iterator from a tabix BED/GFF3 file
 
     Used to avoid chrX vs. X notation issues when extracting data from
     annotation files
@@ -277,7 +277,6 @@ def plot_coverage(
     max_coverage,
     tracktype,
     yaxis_label_fontsize,
-    same_yaxis_labels,
     max_coverage_points,
 ):
     """Plots high and low quality coverage for the region
@@ -314,8 +313,8 @@ def plot_coverage(
     cover_y_highqual = np.array(cover_y_highqual)
     cover_y_all = np.array(cover_y_all)
 
-    if max_coverage > 0 and same_yaxis_labels:
-        max_plot_depth = max_coverage+5
+    if max_coverage > 0:
+        max_plot_depth = max_coverage
     elif cover_y_all.max() > 3 * cover_y_all.mean():
         max_plot_depth = max(
             np.percentile(cover_y_all, 99.5), np.percentile(cover_y_all, 99.5)
@@ -2209,7 +2208,7 @@ def gff_file(transcript_file):
             ext = os.path.splitext(fields[0])[1][1:]
         ext = ext.lower()
         if ext not in options:
-            parser.error("transcript file {} is not in GFF/GFF3 format".format(transcript_file))
+            parser.error("transcript file {} is not in GFF3 format".format(transcript_file))
 
         idx_file = transcript_file + ".tbi"
         if not os.path.isfile(idx_file):
@@ -2427,6 +2426,13 @@ def bed_file(annotation_file):
         required=False,
     )
 
+    parser.add_argument(
+        "--max_coverage",
+        default=0,
+        type=int,
+        help="apply a maximum coverage cutoff. Unlimited by default",
+    )
+
     parser.add_argument(
         "--same_yaxis_scales",
         action="store_true",
@@ -2751,7 +2757,8 @@ def get_read_data(
         read_data["all_pairs"] = downsample_pairs(
             max_depth, z_score, read_data["all_pairs"]
         )
-
+    if not same_yaxis_scales:
+        max_coverage = 0
     return read_data, max_coverage
 
 
@@ -2864,7 +2871,6 @@ def plot_samples(
                 max_coverage,
                 coverage_tracktype,
                 yaxis_label_fontsize,
-                same_yaxis_scales,
                 max_coverage_points,
             )
 
@@ -3292,6 +3298,8 @@ def get_transcript_plan(ranges, transcript_file):
 
     for r in ranges:
         itr = get_tabix_iter(r.chrm, r.start, r.end, transcript_file)
+        if not itr:
+            continue
         for row in itr:
             gene_annotation = row.rstrip().split()
 
@@ -3584,6 +3592,8 @@ def plot(parser):
             options.start_ci,
             options.end_ci,
         )
+    if options.max_coverage:
+        max_coverage = options.max_coverage
 
     # Plot each sample
     current_axis_idx = plot_samples(

samplot/samplot_vcf.py

@@ -104,12 +104,17 @@ def get_format_fields(ids, variant):
     returns:
         list
     """
-    fields = list()
-    for i in ids:
-        fields.append(
-            ["%s=%s" % (i, flatten(j.get(i, ""))) for j in variant.samples.values()]
-        )
-    return zip_lists(fields)
+    sample_format = []
+    for i,sample_fields in enumerate(variant.samples.values()):
+        for field_id in ids:
+            sample_field_val = flatten(sample_fields.get(field_id, ""))
+            if sample_field_val:
+                if len(sample_format) < i+1:
+                    sample_format.append("")
+                else:
+                    sample_format[i] += " "
+                sample_format[i] += "{}={}".format(field_id,sample_field_val)
+    return sample_format
 
 
 def get_format_title(samples, ids, variant):
@@ -175,8 +180,8 @@ def get_overlap(
             [i.split("\t")[2].lower() for i in tabix.fetch(chrom, start, end)]
         )
     except IndexError:
-        # probably not a gff or gtf
-        print("Invalid annotation file specified for --gff")
+        # probably not a gff3
+        print("Invalid annotation file specified for --gff3")
         overlaps = None
     except ValueError:
         if fix_chr:
@@ -415,8 +420,8 @@ def vcf(parser):
     vcf_samples_list = list(vcf_samples)
 
     annotations = None
-    if args.gff:
-        annotations = pysam.TabixFile(args.gff)
+    if args.gff3:
+        annotations = pysam.TabixFile(args.gff3)
 
     filters = [to_exprs(f) for f in args.filter]
 
@@ -721,7 +726,7 @@ def vcf(parser):
             html_template.render(
                 data=tabledata,
                 plot_type=args.output_type,
-                gff="true" if annotations else "false",
+                gff3="true" if annotations else "false",
                 denovo="true" if dn_row else "false",
             ),
             file=fh,
@@ -845,7 +850,7 @@ def add_vcf(parent_parser):
         help="comma separated list of FORMAT fields to include in sample plot title",
     )
     parser.add_argument(
-        "--gff",
+        "--gff3",
         help="genomic regions (.gff with .tbi in same directory) used when building HTML table and table filters",
     )
     parser.add_argument(

samplot/templates/samplot_vcf.html

@@ -201,7 +201,7 @@ <h5>SV Overlaps</h5>
 <script>
     const data = {{ data|tojson }}
     const plot_type = "{{plot_type}}"
-    const annotation = {{ gff }}
+    const annotation = {{ gff3 }}
     const denovo = {{ denovo }}
 
     dc.config.defaultColors(d3.schemeSet1)