Merge pull request #111 from jbelyeu/master

Samplot minor fixes, mostly in VCF

Author: ryanlayer

.circleci/setup.sh

@@ -38,7 +38,13 @@ if [[ ! -d $WORKSPACE/anaconda ]]; then
 
     mkdir -p $WORKSPACE/anaconda/conda-bld/$tag-64
 
-    # step 2: setup channels
+    #step 2: downgrade python3.8 to 3.7
+    if (( $pythonversion == 3 ))
+    then
+        conda install -y python=3.7
+    fi
+
+    # step 3: setup channels
     conda config --system --add channels defaults
     conda config --system --add channels r
     conda config --system --add channels bioconda

samplot/__init__.py

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

samplot/samplot.py

@@ -15,6 +15,8 @@
 import matplotlib.ticker as ticker
 import numpy as np
 import pysam
+import warnings
+warnings.filterwarnings('ignore', 'FixedFormatter should only be used together with FixedLocator')
 from matplotlib.offsetbox import AnchoredText
 
 
@@ -275,7 +277,8 @@ def plot_coverage(
     max_coverage,
     tracktype,
     yaxis_label_fontsize,
-    same_yaxis_labels=False,
+    same_yaxis_labels,
+    max_coverage_points,
 ):
     """Plots high and low quality coverage for the region
 
@@ -290,7 +293,14 @@ def plot_coverage(
 
     for i in range(len(ranges)):
         r = ranges[i]
-        for pos in range(r.start, r.end + 1):
+        region_len = r.end-r.start
+        downsample = 1
+        if region_len > max_coverage_points:
+            downsample = int(region_len / max_coverage_points)
+
+        for i,pos in enumerate(range(r.start, r.end + 1)):
+            if i%downsample !=  0: 
+                continue
             cover_x.append(map_genome_point_to_range_points(ranges, r.chrm, pos))
             if r.chrm in coverage and pos in coverage[r.chrm]:
                 cover_y_all.append(coverage[r.chrm][pos][0] + coverage[r.chrm][pos][1])
@@ -351,12 +361,14 @@ def plot_coverage(
         ax2.fill_between(
             cover_x, cover_y_lowqual, bottom_fill, color="grey", step="pre", alpha=0.15
         )
+    ## tracktype==None also allowed
 
     # number of ticks should be 6 if there's one hp, 3 otherwise
     tick_count = 5 if hp_count == 1 else 2
     tick_count = max(int(max_plot_depth / tick_count), 1)
 
     # set axis parameters
+    #ax2.yaxis.set_major_locator(ticker.FixedLocator(tick_count))
     ax2.yaxis.set_major_locator(ticker.MultipleLocator(tick_count))
     ax2.tick_params(axis="y", colors="grey", labelsize=yaxis_label_fontsize)
     ax2.spines["top"].set_visible(False)
@@ -2133,7 +2145,7 @@ def bam_file(bam):
         "-s",
         "--start",
         type=int,
-        help="Start position of region/variant",
+        help="Start position of region/variant (add multiple for translocation/BND events)",
         action="append",
         required=True,
     )
@@ -2142,15 +2154,15 @@ def bam_file(bam):
         "-e",
         "--end",
         type=int,
-        help="End position of region/variant",
+        help="End position of region/variant (add multiple for translocation/BND events)",
         action="append",
         required=True,
     )
 
     parser.add_argument(
         "-c",
         "--chrom", type=str,
-        help="Chromosome",
+        help="Chromosome (add multiple for translocation/BND events)",
         action="append",
         required=True
     )
@@ -2210,6 +2222,14 @@ def gff_file(transcript_file):
         required=False,
         type=str,
     )
+    
+    parser.add_argument(
+        "--max_coverage_points",
+        help="number of points to plot in coverage axis (downsampled from region size for speed)",
+        required=False,
+        type=int,
+        default=1000,
+    )
 
     def bed_file(annotation_file):
         if not os.path.isfile(annotation_file):
@@ -2250,7 +2270,7 @@ def bed_file(annotation_file):
         "--coverage_tracktype",
         type=str,
         help="type of track to use for low MAPQ " + "coverage plot.",
-        choices=["stack", "superimpose"],
+        choices=["stack", "superimpose", "none"],
         default="stack",
         required=False,
     )
@@ -2779,6 +2799,7 @@ def plot_samples(
     yaxis_label_fontsize,
     annotation_files,
     transcript_file,
+    max_coverage_points,
     max_coverage,
     marker_size,
 ):
@@ -2833,6 +2854,7 @@ def plot_samples(
                 coverage_tracktype,
                 yaxis_label_fontsize,
                 same_yaxis_scales,
+                max_coverage_points,
             )
 
             if len(curr_linked_reads) > 0:
@@ -2913,6 +2935,7 @@ def plot_samples(
             curr_ax.tick_params(axis="y", labelsize=yaxis_label_fontsize)
             # if there's one hp, 6 ticks fit. Otherwise, do 3
             tick_count = 6 if len(hps) == 1 else 3
+            #curr_ax.yaxis.set_major_locator(ticker.FixedLocator(tick_count))
             curr_ax.yaxis.set_major_locator(ticker.LinearLocator(tick_count))
             curr_ax.tick_params(axis="both", length=0)
             curr_ax.set_xticklabels([])
@@ -2950,7 +2973,7 @@ def plot_samples(
             else:
                 sys.stderr.write("Ranges greater than 2 are not supported\n")
                 sys.exit(1)
-
+            
             curr_ax.set_xticklabels(labels, fontsize=xaxis_label_fontsize)
             chrms = [x.chrm for x in ranges]
             curr_ax.set_xlabel("Chromosomal position on " + "/".join(chrms), fontsize=8)
@@ -3565,6 +3588,7 @@ def plot(parser):
         options.yaxis_label_fontsize,
         options.annotation_files,
         options.transcript_file,
+        options.max_coverage_points,
         max_coverage,
         marker_size,
     )
@@ -3598,7 +3622,7 @@ def plot(parser):
     matplotlib.rcParams["agg.path.chunksize"] = 100000
     plt.tight_layout(pad=0.8, h_pad=0.1, w_pad=0.1)
     try:
-        plt.savefig(output_file)
+        plt.savefig(output_file, dpi=100)
     except Exception as e:
         print(
             "Failed to save figure " + output_file

samplot/samplot_vcf.py

@@ -440,30 +440,56 @@ def vcf(parser):
 
     out_file = open(args.command_file, "w")
 
-    for variant in vcf:
+    for var_count,variant in enumerate(vcf):
+        try:
+            translocation_chrom = variant.info.get("CHR2")
+        except:
+            translocation_chrom = None
         svtype = variant.info.get("SVTYPE", "SV")
         if args.important_regions:
             if not var_in_important_regions(
                 important_regions, variant.chrom, variant.start, variant.stop
             ):
+                if args.debug:
+                    print("Skipping {} at {}:{}-{}, outside important_regions coordinates".format(
+                        svtype, variant.chrom, variant.start, variant.stop),file=sys.stderr)
                 continue
 
         if svtype in ("INS"):
+            if args.debug:
+                print("Skipping {} at {}:{}-{}, INS type not supported".format(
+                    svtype, variant.chrom, variant.start, variant.stop),file=sys.stderr)
             continue
-        if variant.stop - variant.start > args.max_mb * 1000000:
+        if args.max_mb and (variant.stop - variant.start > args.max_mb * 1000000):
+            if args.debug:
+                print("Skipping {} at {}:{}-{}, variant length greater than max_mb".format(
+                    svtype, variant.chrom, variant.start, variant.stop),file=sys.stderr)
             continue
-        if variant.stop - variant.start < args.min_bp:
+        
+        if (variant.stop - variant.start < args.min_bp) and translocation_chrom is None:
+            if args.debug:
+                print("Skipping {} at {}:{}-{}, variant length shorter than min_bp".format(
+                    svtype, variant.chrom, variant.start, variant.stop),file=sys.stderr)
             continue
 
         gts = [s.get("GT", (None, None)) for s in variant.samples.values()]
 
         if sum(None in g for g in gts) >= args.min_call_rate * len(vcf_samples):
+            if args.debug:
+                print("Skipping {} at {}:{}-{}, call rate of variant below min_call_rate".format(
+                    svtype, variant.chrom, variant.start, variant.stop),file=sys.stderr)
             continue
         if args.max_hets:
             # requisite hets/hom-alts
             if sum(sum(x) >= 1 for x in gts if not None in x) > args.max_hets:
+                if args.debug:
+                    print("Skipping {} at {}:{}-{}, more than max_hets heterozygotes".format(
+                        svtype, variant.chrom, variant.start, variant.stop),file=sys.stderr)
                 continue
         if not any(sum(x) > 0 for x in gts if not None in x):
+            if args.debug:
+                print("Skipping {} at {}:{}-{}, no samples have non-ref genotypes".format(
+                    svtype, variant.chrom, variant.start, variant.stop),file=sys.stderr)
             continue
 
         test_idxs = [i for i, gt in enumerate(gts) if not None in gt and sum(gt) > 0]
@@ -493,7 +519,10 @@ def vcf(parser):
         for i in idxs:
             if vcf_samples[i] in names_to_bams:
                 variant_samples.append(vcf_samples[i])
-        if len(variant_samples) <= 0:
+        if len(variant_samples) == 0:
+            if args.debug:
+                print("Skipping {} at {}:{}-{}, no samples with matched alignment files have variant".format(
+                    svtype, variant.chrom, variant.start, variant.stop),file=sys.stderr)
             continue
 
         bams = [names_to_bams[s] for s in variant_samples]
@@ -510,6 +539,9 @@ def vcf(parser):
                     is_dn.append(sample.id)
 
         if len(is_dn) <= 0 and args.dn_only:
+            if args.debug:
+                print("Skipping {} at {}:{}-{}, dn_only selected and no de novos found".format(
+                    svtype, variant.chrom, variant.start, variant.stop),file=sys.stderr)
             continue
 
         # save these for the html.
@@ -581,6 +613,7 @@ def vcf(parser):
 
         data_dict = {
             "chrom": variant.chrom,
+            "chrom2": translocation_chrom,
             "start": variant.start,
             "end": variant.stop,
             "svtype": svtype,
@@ -590,15 +623,17 @@ def vcf(parser):
         }
         if annotations:
             data_dict["overlaps"] = get_overlap(
-                annotations, variant.chrom, variant.start, variant.stop
+                annotations, variant.chrom, variant.start, variant.stop, translocation_chrom
             )
         if dn_row != "":
             data_dict["dn"] = ",".join(is_dn)
+        if translocation_chrom is None:
+            template = "{svtype}_{chrom}_{start}_{end}.{itype}"
+        else:
+            template = "{svtype}_{chrom}_{start}_{chrom2}_{end}.{itype}"
         fig_path = os.path.join(
             args.out_dir,
-            "{svtype}_{chrom}_{start}_{end}.{itype}".format(
-                itype=args.output_type, **data_dict
-            ),
+            template.format(itype=args.output_type, **data_dict),
         )
         tabledata.append(data_dict)
 
@@ -633,8 +668,22 @@ def vcf(parser):
             else:
                 title_list.append(variant_sample)
 
+        template = ("samplot plot {extra_args} -z {z} -n {titles}"
+                + "{cipos} {ciend} {svtype} -c {chrom} -s {start} "
+                + "-e {end} -o {fig_path} -d {downsample} -b {bams}")
+        start = variant.start
+        stop = variant.stop
+        start2 = None
+        stop2 = None
+        if translocation_chrom is not None:
+            template += " -c {chrom2} -s {start2} -e {end2}"
+            start2 = stop
+            stop2 = start2 + 1
+            stop = start + 1
+        template += "\n"
+
         out_file.write(
-            "samplot plot {extra_args} -z {z} -n {titles} {cipos} {ciend} {svtype} -c {chrom} -s {start} -e {end} -o {fig_path} -d {downsample} -b {bams}\n".format(
+            template.format(
                 extra_args=" ".join(pass_through_args),
                 bams=" ".join(bams),
                 titles=" ".join(title_list),
@@ -644,11 +693,16 @@ def vcf(parser):
                 svtype="-t " + svtype if svtype != "SV" else "",
                 fig_path=fig_path,
                 chrom=variant.chrom,
-                start=variant.start,
-                end=variant.stop,
+                start=start,
+                end=stop,
                 downsample=args.downsample,
+                chrom2=translocation_chrom,
+                start2=start2,
+                end2=stop2,
             )
         )
+    if args.debug:
+        print("VCF entry count:",var_count+1 ,file=sys.stderr)
 
     if args.command_file:
         out_file.close()
@@ -750,7 +804,6 @@ def add_vcf(parent_parser):
         "--max_mb",
         type=int,
         help="skip variants longer than this many megabases",
-        default=1,
     )
     parser.add_argument(
         "--min_bp",
@@ -802,7 +855,12 @@ def add_vcf(parent_parser):
         default=False,
         action="store_true",
     )
-
+    parser.add_argument(
+        "--debug",
+        help="prints out the reason each skipped variant entry is skipped",
+        default=False,
+        action="store_true",
+    )
     parser.set_defaults(func=vcf)