1. determine bin coverage
    2. determine rRNA gene coverage
    3. compare
4955
def copies(mapping, s2bins, rna, min_rna = 800, mismatches = 0):
    
    cov = {} # cov[scaffold] = [bases, length]
    s2bins, bins2s = parse_s2bins(s2bins)
    rna_cov = parse_rna(rna, s2bins, min_rna)
    s2bins, bins2s = filter_missing_rna(s2bins, bins2s, rna_cov)
    # count bases mapped to scaffolds and rRNA gene regions
    for line in mapping:
        line = line.strip().split()
        # get scaffold lengths
        if line[0].startswith('@'):
            if line[0].startswith('@SQ') is False:
                continue
            s = line[1].split(':')[1]
            l = int(line[2].split(':')[1])
            # check if scaffold is binned
            if s not in s2bins:
                continue
            if s not in cov:
                cov[s] = [0, l]
        # check mismatch threshold
        mm = count_mismatches(line)
        if mm is False or mm > mismatches:
            continue
        # check that scaffold is in bin
        s, bases = line[2], len(line[9])
        if s not in cov:
            continue
        cov[s][0] += bases 
        rna_cov = rna_bases(rna_cov, s, bases, line) 
    print('# mismatches threshold: %s' % (mismatches))
    header = ['#rRNA scaffold', 'rRNA genes >=%sbp on scaffold' % (min_rna), \
            'rRNA coverage', \
            'bin', 'bin info', 'bin coverage', \
            'rRNAs >=%sbp in bin' % (min_rna), \
            'rRNA coverage/bin coverage', \
            'estimated number of copies']
    print('\t'.join(header))
    for bin, scaffolds in list(bins2s.items()):
        rna_count = sum([len(rna_cov[s][2]) for s in scaffolds if s in rna_cov])
        for s in scaffolds:
            if s not in rna_cov:
                continue
            out = []
            counts = rna_cov[s]
            bin_cov = calc_bin_cov(bins2s[bin], cov)
            num_genes = len(counts[2])
            rna_coverage = float(float(counts[0])/float(counts[1]))
            if bin_cov == 0:
                rna_div_bin = 0
            else:
                rna_div_bin = float(rna_coverage/bin_cov)
            est = int(max([rna_count, counts, rna_div_bin]))
            out = [s, num_genes, rna_coverage, bin, bin_cov, rna_count, rna_div_bin, est]
            print('\t'.join([str(i) for i in out]))
876133
def _average_genome_coverage(data, bam_file):
    
    total = sum([c.size for c in ref.file_contigs(dd.get_ref_file(data), data["config"])])
    read_counts = sum(x.aligned for x in bam.idxstats(bam_file, data))
    with pysam.Samfile(bam_file, "rb") as pysam_bam:
        read_size = np.median(list(itertools.islice((a.query_length for a in pysam_bam.fetch()), int(1e7))))
    avg_cov = float(read_counts * read_size) / total
    return avg_cov
254543
def get_coords(data):
    
    for category, vtypes in [("LOH", {"LOSS", "HETEROZYGOSITY"}),
                             ("amplification", {"AMPLIFICATION"})]:
        out = tz.get_in([category, dd.get_genome_build(data)], _COORDS, {})
        priority_file = dd.get_svprioritize(data)
        if priority_file:
            if os.path.basename(priority_file).find("civic") >= 0:
                for chrom, start, end, gene in _civic_regions(priority_file, vtypes, dd.get_disease(data)):
                    out[gene] = (chrom, start, end)
            elif os.path.basename(priority_file).find(".bed") >= 0:
                for line in utils.open_gzipsafe(priority_file):
                    parts = line.strip().split("\t")
                    if len(parts) >= 4:
                        chrom, start, end, gene = parts[:4]
                        out[gene] = (chrom, int(start), int(end))
        yield category, out
544628
function (collector, options) {
	var summaries = [];
	var thresholds = options.thresholds;

	var files = collector.files();
	files.forEach(function(file) {
		summaries.push(istanbul.utils.summarizeFileCoverage(
				collector.fileCoverageFor(file)));
	});

	var finalSummary = istanbul.utils.mergeSummaryObjects.apply(null,
			summaries);

	grunt.util._.each(thresholds, function (threshold, metric) {
		var actual = finalSummary[metric];
		if(!actual) {
			grunt.warn('unrecognized metric: ' + metric);
		}
		if(actual.pct < threshold) {
			grunt.warn('expected ' + metric + ' coverage to be at least '
					+ threshold + '% but was ' + actual.pct + '%');
		}
	});
}
1032671
def _determine_stream_spreadLB(self,simple=_USESIMPLE,
                                   ro=None,vo=None,
                                   R0=None,Zsun=None,vsun=None):
        """"""
        if not hasattr(self,'_allErrCovs'):
            self._determine_stream_spread(simple=simple)
        if ro is None:
            ro= self._ro
        if vo is None:
            vo= self._vo
        if R0 is None:
            R0= self._R0
        if Zsun is None:
            Zsun= self._Zsun
        if vsun is None:
            vsun= self._vsun
        allErrCovsLB= numpy.empty_like(self._allErrCovs)
        obs= [R0,0.,Zsun]
        obs.extend(vsun)
        obskwargs= {}
        obskwargs['ro']= ro
        obskwargs['vo']= vo
        obskwargs['obs']= obs
        self._ErrCovsLBScale= [180.,90.,
                               self._progenitor.dist(**obskwargs),
                               numpy.fabs(self._progenitor.vlos(**obskwargs)),
                               numpy.sqrt(self._progenitor.pmll(**obskwargs)**2.
                                          +self._progenitor.pmbb(**obskwargs)**2.),
                               numpy.sqrt(self._progenitor.pmll(**obskwargs)**2.
                                          +self._progenitor.pmbb(**obskwargs)**2.)]
        allErrCovsEigvalLB= numpy.empty((len(self._thetasTrack),6))
        allErrCovsEigvecLB= numpy.empty_like(self._allErrCovs)
        eigDir= numpy.array([numpy.array([1.,0.,0.,0.,0.,0.]) for ii in range(6)])
        for ii in range(self._nTrackChunks):
            tjacXY= bovy_coords.galcenrect_to_XYZ_jac(*self._ObsTrackXY[ii])
            tjacLB= bovy_coords.lbd_to_XYZ_jac(*self._ObsTrackLB[ii],
                                               degree=True)
            tjacLB[:3,:]/= ro
            tjacLB[3:,:]/= vo
            for jj in range(6):
                tjacLB[:,jj]*= self._ErrCovsLBScale[jj]
            tjac= numpy.dot(numpy.linalg.inv(tjacLB),tjacXY)
            allErrCovsLB[ii]=\
                numpy.dot(tjac,numpy.dot(self._allErrCovsXY[ii],tjac.T))
            #Eigen decomposition for interpolation
            teig= numpy.linalg.eig(allErrCovsLB[ii])
            #Sort them to match them up later
            sortIndx= numpy.argsort(teig[0])
            allErrCovsEigvalLB[ii]= teig[0][sortIndx]
            #Make sure the eigenvectors point in the same direction
            for jj in range(6):
                if numpy.sum(eigDir[jj]*teig[1][:,sortIndx[jj]]) < 0.:
                    teig[1][:,sortIndx[jj]]*= -1.
                eigDir[jj]= teig[1][:,sortIndx[jj]]
            allErrCovsEigvecLB[ii]= teig[1][:,sortIndx]
        self._allErrCovsLBUnscaled= allErrCovsLB
        #Interpolate the allErrCovsLB covariance matrices along the interpolated track
        #Interpolate the eigenvalues
        interpAllErrCovsEigvalLB=\
            [interpolate.InterpolatedUnivariateSpline(self._thetasTrack,
                                                      allErrCovsEigvalLB[:,ii],
                                                      k=3) for ii in range(6)]
        #Now build the interpolated allErrCovsXY using slerp
        interpolatedAllErrCovsLB= numpy.empty((len(self._interpolatedThetasTrack),
                                               6,6))
        interpolatedEigval=\
            numpy.array([interpAllErrCovsEigvalLB[ii](self._interpolatedThetasTrack) for ii in range(6)]) #6,ninterp
        #Interpolate in chunks
        interpolatedEigvec= numpy.empty((len(self._interpolatedThetasTrack),
                                         6,6))
        for ii in range(self._nTrackChunks-1):
            slerpOmegas=\
                [numpy.arccos(numpy.sum(allErrCovsEigvecLB[ii,:,jj]*allErrCovsEigvecLB[ii+1,:,jj])) for jj in range(6)]
            slerpts= (self._interpolatedThetasTrack-self._thetasTrack[ii])/\
                (self._thetasTrack[ii+1]-self._thetasTrack[ii])
            slerpIndx= (slerpts >= 0.)*(slerpts <= 1.)
            for jj in range(6):
                for kk in range(6):
                    interpolatedEigvec[slerpIndx,kk,jj]=\
                    (numpy.sin((1-slerpts[slerpIndx])*slerpOmegas[jj])*allErrCovsEigvecLB[ii,kk,jj]
                     +numpy.sin(slerpts[slerpIndx]*slerpOmegas[jj])*allErrCovsEigvecLB[ii+1,kk,jj])/numpy.sin(slerpOmegas[jj])
        for ii in range(len(self._interpolatedThetasTrack)):
            interpolatedAllErrCovsLB[ii]=\
                numpy.dot(interpolatedEigvec[ii],
                          numpy.dot(numpy.diag(interpolatedEigval[:,ii]),
                                    interpolatedEigvec[ii].T))
        self._interpolatedAllErrCovsLBUnscaled= interpolatedAllErrCovsLB
        #Also calculate the (l,b,..) -> (X,Y,..) Jacobian at all of the interpolated and not interpolated points
        trackLogDetJacLB= numpy.empty_like(self._thetasTrack)
        interpolatedTrackLogDetJacLB=\
            numpy.empty_like(self._interpolatedThetasTrack)
        for ii in range(self._nTrackChunks):
            tjacLB= bovy_coords.lbd_to_XYZ_jac(*self._ObsTrackLB[ii],
                                                degree=True)
            trackLogDetJacLB[ii]= numpy.log(numpy.linalg.det(tjacLB))
        self._trackLogDetJacLB= trackLogDetJacLB
        for ii in range(len(self._interpolatedThetasTrack)):
            tjacLB=\
                bovy_coords.lbd_to_XYZ_jac(*self._interpolatedObsTrackLB[ii],
                                            degree=True)
            interpolatedTrackLogDetJacLB[ii]=\
                numpy.log(numpy.linalg.det(tjacLB))
        self._interpolatedTrackLogDetJacLB= interpolatedTrackLogDetJacLB
        return None
1318795
def parse_reports(self):
    

    # General stats - genome_results.txt
    self.qualimap_bamqc_genome_results = dict()
    for f in self.find_log_files('qualimap/bamqc/genome_results'):
        parse_genome_results(self, f)
    self.qualimap_bamqc_genome_results = self.ignore_samples(self.qualimap_bamqc_genome_results)

    # Coverage - coverage_histogram.txt
    self.qualimap_bamqc_coverage_hist = dict()
    for f in self.find_log_files('qualimap/bamqc/coverage', filehandles=True):
        parse_coverage(self, f)
    self.qualimap_bamqc_coverage_hist = self.ignore_samples(self.qualimap_bamqc_coverage_hist)

    # Insert size - insert_size_histogram.txt
    self.qualimap_bamqc_insert_size_hist = dict()
    for f in self.find_log_files('qualimap/bamqc/insert_size', filehandles=True):
        parse_insert_size(self, f)
    self.qualimap_bamqc_insert_size_hist = self.ignore_samples(self.qualimap_bamqc_insert_size_hist)

    # GC distribution - mapped_reads_gc-content_distribution.txt
    self.qualimap_bamqc_gc_content_dist = dict()
    self.qualimap_bamqc_gc_by_species = dict()  # {'HUMAN': data_dict, 'MOUSE': data_dict}
    for f in self.find_log_files('qualimap/bamqc/gc_dist', filehandles=True):
        parse_gc_dist(self, f)
    self.qualimap_bamqc_gc_by_species = self.ignore_samples(self.qualimap_bamqc_gc_by_species)

    num_parsed = max(
        len(self.qualimap_bamqc_genome_results),
        len(self.qualimap_bamqc_coverage_hist),
        len(self.qualimap_bamqc_insert_size_hist),
        len(self.qualimap_bamqc_gc_content_dist)
    )
    # Go no further if nothing found
    if num_parsed == 0:
        return 0

    try:
        covs = config.qualimap_config['general_stats_coverage']
        assert type(covs) == list
        assert len(covs) > 0
        covs = [str(i) for i in covs]
        log.debug("Custom Qualimap thresholds: {}".format(", ".join([i for i in covs])))
    except (AttributeError, TypeError, AssertionError):
        covs = [1, 5, 10, 30, 50]
        covs = [str(i) for i in covs]
        log.debug("Using default Qualimap thresholds: {}".format(", ".join([i for i in covs])))
    self.covs = covs

    # Make the plots for the report
    report_sections(self)

    # Set up the general stats table
    general_stats_headers(self)

    # Return the number of reports we found
    return num_parsed
66493
def use_vqsr(algs, call_file=None):
    
    from bcbio.variation import vcfutils
    vqsr_callers = set(["gatk", "gatk-haplotype"])
    vqsr_sample_thresh = 50
    vqsr_supported = collections.defaultdict(int)
    coverage_intervals = set([])
    for alg in algs:
        callers = alg.get("variantcaller")
        if isinstance(callers, six.string_types):
            callers = [callers]
        if not callers:  # no variant calling, no VQSR
            continue
        if "vqsr" in (alg.get("tools_off") or []):  # VQSR turned off
            continue
        for c in callers:
            if c in vqsr_callers:
                if "vqsr" in (alg.get("tools_on") or []):  # VQSR turned on:
                    vqsr_supported[c] += 1
                    coverage_intervals.add("genome")
                # Do not try VQSR for gVCF inputs
                elif call_file and vcfutils.is_gvcf_file(call_file):
                    pass
                else:
                    coverage_intervals.add(alg.get("coverage_interval", "exome").lower())
                    vqsr_supported[c] += 1
    if len(vqsr_supported) > 0:
        num_samples = max(vqsr_supported.values())
        if "genome" in coverage_intervals or num_samples >= vqsr_sample_thresh:
            return True
    return False
1311611
def coverage_lineplot (self):
        
        # Add line graph to section
        data = list()
        data_labels = list()
        if len(self.rna_seqc_norm_high_cov) > 0:
            data.append(self.rna_seqc_norm_high_cov)
            data_labels.append({'name': 'High Expressed'})
        if len(self.rna_seqc_norm_medium_cov) > 0:
            data.append(self.rna_seqc_norm_medium_cov)
            data_labels.append({'name': 'Medium Expressed'})
        if len(self.rna_seqc_norm_low_cov) > 0:
            data.append(self.rna_seqc_norm_low_cov)
            data_labels.append({'name': 'Low Expressed'})
        pconfig = {
            'id': 'rna_seqc_mean_coverage_plot',
            'title': 'RNA-SeQC: Gene Body Coverage',
            'ylab': '% Coverage',
            'xlab': "Gene Body Percentile (5' -> 3')",
            'xmin': 0,
            'xmax': 100,
            'tt_label': "<strong>{point.x}% from 5'</strong>: {point.y:.2f}",
            'data_labels': data_labels
        }
        if len(data) > 0:
            self.add_section (
                name = 'Gene Body Coverage',
                anchor = 'rseqc-rna_seqc_mean_coverage',
                helptext = 'The metrics are calculated across the transcripts with tiered expression levels.',
                plot = linegraph.plot(data, pconfig)
            )
389632
def determine_coverage(cls, coverage_file):
        
        lines = []
        source_file = 'ERROR'
        for line in coverage_file:
            m = title_re.match(line)
            if m:
                if m.group(2) == '100':
                    return ('', [])
                source_file = m.group(1)
                continue
            m = source_re.match(line)
            if m:
                lines.append(int(m.group(1)))
                continue
            if end_re.match(line):
                break
        line_ranges = cls.make_ranges(lines)
        return (source_file, line_ranges)
1140752
def get_orfs(_orf_length = 100, prediction = @prediction)
      '-' if prediction.type != 'nucleotide'

      seq = Bio::Sequence::NA.new(prediction.raw_sequence)

      result = {}
      key = 0
      (1..6).each do |f|
        s = seq.translate(f)
        f = -1 if f == 4
        f = -2 if f == 5
        f = -3 if f == 6
        s.scan(/(\w{30,})/) do |_orf|
          orf_start = $LAST_MATCH_INFO.offset(0)[0] + 1
          orf_end   = $LAST_MATCH_INFO.offset(0)[1] + 1
          coverage = (((orf_end - orf_start) / s.length.to_f) * 100).ceil
          # reduce the orf_end and the translated length by 2% to increase the
          #   width between ORFs on the plot
          chopping = s.length * 0.02
          orf_end = (orf_end.to_f - chopping).floor
          translated_length = (s.length - chopping).ceil
          key += 1
          result[key] = { frame: f, orf_start: orf_start, orf_end: orf_end,
                          coverage: coverage,
                          translated_length: translated_length }
        end
      end
      result
    end
95639
def effective_genome_size(fasta, read_length, nb_cores, tmpdir="/tmp"):
    # type: (str, int, int, str) -> None
    """"""

    idx = Fasta(fasta)

    genome_length = sum([len(c) for c in idx])

    logging.info("Temporary directory: " + tmpdir)
    logging.info("File analyzed: " + fasta)
    logging.info("Genome length: " + str(genome_length))
    print("File analyzed: ", fasta)
    print("Genome length: ", genome_length)

    chromosomes = ", ".join([c.name for c in idx])

    if "_" in chromosomes:
        print("Warning. The following chromosomes are part of your genome:\n",
              chromosomes.replace(">", "") + "\n",
              file=sys.stderr)
        print(
            "You probably want to remove all chromosomes in your fasta containing '_' for the effective genome size computation to be accurate.",
            file=sys.stderr)

    # if tmpdir is None:
    #     try:
    #         tmpdir = os.environ['TMPDIR']
    #     except KeyError:
    #         tmpdir = '/tmp'
    output_file = os.path.join(tmpdir, '{1}.jf'.format(read_length,
                                                       basename(fasta)))
    atexit.register(
        lambda: call("rm {output_file}".format(output_file=output_file), shell=True))

    call(
        "jellyfish count -t {nb_cores} -m {read_length} -s {genome_length} -L 1 -U 1 --out-counter-len 1 --counter-len 1 {fasta} -o {output_file}".format(
            **vars()),
        shell=True)

    stats = check_output("jellyfish stats {output_file}".format(
        output_file=output_file),
                         shell=True)

    unique_kmers = int(stats.split()[1])

    effective_genome_size = unique_kmers / genome_length

    logging.info("Number unique {read_length}-mers: ".format(
        read_length=read_length) + str(unique_kmers))
    logging.info("Effective genome size: " + str(effective_genome_size))
    print("Number unique {read_length}-mers: ".format(read_length=read_length),
          unique_kmers)
    print("Effective genome size: ", effective_genome_size)
    assert effective_genome_size < 1, "Something wrong happened, effective genome size over 1!"
947460
def estimate_percent_contamination(contamination_report_file):
    
    contam_levels = list()
    with open(contamination_report_file) as csvfile:
        reader = csv.DictReader(csvfile)
        for row in reader:
            lowest_count = 99999
            base_counts = row['Bases'].split(';')
            for count in base_counts:
                num_bases = int(count.split(':')[1])
                if num_bases < lowest_count:
                    lowest_count = num_bases
            total_coverage = int(row['Coverage'])
            contam_levels.append(lowest_count*100/total_coverage)
    return '%.2f' % (np.mean(contam_levels)), '%.2f' % np.std(contam_levels)
1000316
def _calculate_sv_coverage_cnvkit(data, work_dir):
    
    from bcbio.variation import coverage
    from bcbio.structural import annotate
    out_target_file = os.path.join(work_dir, "%s-target-coverage.cnn" % dd.get_sample_name(data))
    out_anti_file = os.path.join(work_dir, "%s-antitarget-coverage.cnn" % dd.get_sample_name(data))
    if ((not utils.file_exists(out_target_file) or not utils.file_exists(out_anti_file)) and
          (dd.get_align_bam(data) or dd.get_work_bam(data))):
        target_cov = coverage.run_mosdepth(data, "target", tz.get_in(["regions", "bins", "target"], data))
        anti_cov = coverage.run_mosdepth(data, "antitarget", tz.get_in(["regions", "bins", "antitarget"], data))
        target_cov_genes = annotate.add_genes(target_cov.regions, data, max_distance=0)
        out_target_file = _add_log2_depth(target_cov_genes, out_target_file, data)
        out_anti_file = _add_log2_depth(anti_cov.regions, out_anti_file, data)
    return out_target_file, out_anti_file
246430
def count(args):
    
    p = OptionParser(count.__doc__)
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    coveragefile, fastafile = args

    countsfile = coveragefile.split(".")[0] + ".bin"
    if op.exists(countsfile):
        logging.error("`{0}` file exists. Remove before proceed."\
                .format(countsfile))
        return

    fastasize, sizes, offsets = get_offsets(fastafile)
    logging.debug("Initialize array of uint8 with size {0}".format(fastasize))
    ar = np.zeros(fastasize, dtype=np.uint8)

    update_array(ar, coveragefile, sizes, offsets)

    ar.tofile(countsfile)
    logging.debug("Array written to `{0}`".format(countsfile))
707808
def post_build_geo_coverage(self):
        """"""

        spaces = set()
        grains = set()

        def resolve(term):
            places = list(self.library.search.search_identifiers(term))

            if not places:
                raise BuildError(
                    "Failed to find space identifier '{}' in full text identifier search".format(term))

            return places[0].vid

        if self.metadata.about.space:  # From the bundle metadata
            spaces.add(resolve(self.metadata.about.space))

        if self.metadata.about.grain:  # From the bundle metadata
            grains.add(self.metadata.about.grain)

        if self.identity.bspace:  # And from the bundle name
            spaces.add(resolve(self.identity.bspace))

        # From all of the partitions
        for p in self.partitions.all:
            if 'geo_coverage' in p.record.data:
                for space in p.record.data['geo_coverage']:
                    spaces.add(space)

            if 'geo_grain' in p.record.data:
                for grain in p.record.data['geo_grain']:
                    grains.add(grain)

        def conv_grain(g):
            """Some grain are expressed as summary level names, not gvids."""
            try:
                c = GVid.get_class(g)
                return b(c().summarize())
            except NotASummaryName:
                return g

        self.metadata.coverage.geo = sorted(spaces)
        self.metadata.coverage.grain = sorted(conv_grain(g) for g in grains)

        self.metadata.write_to_dir()
400019
def _combine_coverages(items, work_dir, input_backs=None):
    
    out_file = os.path.join(work_dir, "sample_coverages.txt")
    if not utils.file_exists(out_file):
        with file_transaction(items[0], out_file) as tx_out_file:
            with open(tx_out_file, 'w') as out_f:
                for data in items:
                    cov_file = tz.get_in(["depth", "bins", "seq2c"], data)
                    with open(cov_file) as cov_f:
                        out_f.write(cov_f.read())
                if input_backs:
                    for input_back in input_backs:
                        with open(input_back) as in_handle:
                            for line in in_handle:
                                if len(line.split()) >= 4:
                                    out_f.write(line)
    return out_file
339880
def calculate(bam_file, data, sv_bed):
    
    params = {"min": dd.get_coverage_depth_min(data)}
    variant_regions = dd.get_variant_regions_merged(data)
    if not variant_regions:
        variant_regions = _create_genome_regions(data)
    # Back compatible with previous pre-mosdepth callable files
    callable_file = os.path.join(utils.safe_makedir(os.path.join(dd.get_work_dir(data), "align",
                                                                 dd.get_sample_name(data))),
                                 "%s-coverage.callable.bed" % (dd.get_sample_name(data)))
    if not utils.file_uptodate(callable_file, bam_file):
        vr_quantize = ("0:1:%s:" % (params["min"]), ["NO_COVERAGE", "LOW_COVERAGE", "CALLABLE"])
        to_calculate = [("variant_regions", variant_regions,
                         vr_quantize, None, "coverage_perbase" in dd.get_tools_on(data)),
                        ("sv_regions", bedutils.clean_file(sv_bed, data, prefix="svregions-"),
                         None, None, False),
                        ("coverage", bedutils.clean_file(dd.get_coverage(data), data, prefix="cov-"),
                         None, DEPTH_THRESHOLDS, False)]
        depth_files = {}
        for target_name, region_bed, quantize, thresholds, per_base in to_calculate:
            if region_bed:
                cur_depth = {}
                depth_info = run_mosdepth(data, target_name, region_bed, quantize=quantize, thresholds=thresholds,
                                          per_base=per_base)
                for attr in ("dist", "regions", "thresholds", "per_base"):
                    val = getattr(depth_info, attr, None)
                    if val:
                        cur_depth[attr] = val
                depth_files[target_name] = cur_depth
                if target_name == "variant_regions":
                    callable_file = depth_info.quantize
    else:
        depth_files = {}
    final_callable = _subset_to_variant_regions(callable_file, variant_regions, data)
    return final_callable, depth_files