"""
calculate coverage across a list of regions
"""
import tempfile
import subprocess
import sys
import pandas as pd
if sys.version_info[0] < 3:
    from StringIO import StringIO
else:
    from io import StringIO
import matplotlib.pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
import seaborn as sns
from bcbio.utils import rbind, file_exists
from bcbio.distributed.transaction import file_transaction

def _calc_regional_coverage(in_bam, chrom, start, end, samplename):
    """
    given a BAM and a region, calculate the coverage for each base in that
    region. returns a pandas dataframe of the format:

    chrom position coverage name

    where the samplename column is the coverage at chrom:position
    """
    region_file = tempfile.NamedTemporaryFile(delete=False).name
    with open(region_file, "w") as out_handle:
        out_handle.write("%s\t%s\t%s\n" % (chrom, start, end))

    cmd = ("bedtools coverage -abam {in_bam} -b {region_file} -d")
    out = subprocess.check_output(cmd.format(**locals()), shell=True)
    names = ["chom", "start", "end", "offset", "coverage"]
    df = pd.io.parsers.read_table(StringIO(out), sep="\t", header=None,
                                  names=names)
    df["sample"] = samplename
    df["chrom"] = chrom
    df["position"] = df["start"] + df["offset"] - 1
    return df[["chrom", "position", "coverage", "sample"]]

def _combine_regional_coverage(in_bams, samplenames, chrom, start, end):
    """
    given a list of bam files, sample names and a region, calculate the
    coverage in the region for each of the samples and return a tidy pandas
    dataframe of the format:

    chrom position coverage name
    """
    in_bams = [dd.get_work_bam(x) for x in samples]
    samplenames = [dd.get_sample_name(x) for x in samples]
    dfs = [_calc_regional_coverage(bam, chrom, start, end, sample) for bam, sample
           in zip(in_bams, samplenames)]
    return rbind(dfs)

def plot_multiple_regions_coverage(samples, regions, out_file):
    """
    given a list of bcbio samples and a list of tuples of regions of the format
    (chrom, start, end), make a plot of the coverage in the regions for
    the set of samples
    """
    if file_exists(out_file):
        return out_file
    in_bams = [dd.get_work_bam(x) for x in samples]
    samplenames = [dd.get_sample_name(x) for x in samples]
    with file_transaction(out_file) as tx_out_file:
        with PdfPages(tx_out_file) as pdf_out:
            sns.despine()
            for chrom, start, end in regions:
                df = _combine_regional_coverage(in_bams, samplenames, chrom,
                                                start, end)
                plot = sns.tsplot(df, time="position", unit="chrom",
                                  value="coverage", condition="sample")
                plt.title("{chrom}:{start}-{end}".format(**locals()))
                pdf_out.savefig(plot.get_figure())
    return out_file

@@ -1691,119 +1691,8 @@
 %22%22%22%0A
-    in_bams = %5Bdd.get_work_bam(x) for x in samples%5D%0A    samplenames = %5Bdd.get_sample_name(x) for x in samples%5D%0A
     
@@ -1834,16 +1834,249 @@
 d(dfs)%0A%0A
+def _bed_to_regions(bed_file):%0A    regions = %5B%5D%0A    with open(bed_file) as in_handle:%0A        for line in in_handle:%0A            tokens = line.split()%0A            regions.append((tokens%5B0%5D, tokens%5B1%5D, tokens%5B2%5D))%0A    return regions%0A%0A
 def plot
@@ -2115,24 +2115,46 @@
 es, 
-regions, out_fil
+out_file, region_bed=None, regions=Non
 e):%0A
@@ -2201,16 +2201,41 @@
 s and a 
+bed file of regions or a 
 list of 
@@ -2240,16 +2240,20 @@
 f tuples
+%0A   
  of regi
@@ -2256,27 +2256,24 @@
 regions 
-of 
 the form
 at%0A    (
@@ -2264,22 +2264,16 @@
 the form
-at%0A   
  (chrom,
@@ -2284,17 +2284,16 @@
 rt, end)
-,
  make a 
@@ -2315,24 +2315,28 @@
 erage in the
+%0A   
  regions for
@@ -2335,20 +2335,16 @@
 ions for
-%0A   
  the set
@@ -2359,24 +2359,170 @@
 les%0A    %22%22%22%0A
+    assert (region_bed or regions), (%22A region bed file or a list of tuples of %22%0A                                     %22regions must be provided%22)%0A
     if file_
@@ -2670,24 +2670,90 @@
 in samples%5D%0A
+    if not regions:%0A        regions = _bed_to_regions(region_bed)%0A
     with fil
@@ -3310,16 +3310,44 @@
 gure())%0A
+                plt.close()%0A
     retu

Add option to provide a BED file of regions to plot.