Source code for cnvlib.autobin

"""Estimate reasonable bin sizes from BAM read counts or depths."""
import logging
import os
import tempfile

import numpy as np
import pandas as pd
from skgenome import tabio, GenomicArray as GA

from . import coverage, samutil
from .antitarget import compare_chrom_names
from .descriptives import weighted_median

[docs]def midsize_file(fnames): """Select the median-size file from several given filenames. If an even number of files is given, selects the file just below the median. """ return sorted(fnames, key=lambda f: os.stat(f).st_size )[len(fnames) // 2 - 1]
[docs]def do_autobin(bam_fname, method, targets=None, access=None, bp_per_bin=100000., target_min_size=20, target_max_size=50000, antitarget_min_size=500, antitarget_max_size=1000000, fasta=None): """Quickly calculate reasonable bin sizes from BAM read counts. Parameters ---------- bam_fname : string BAM filename. method : string One of: 'wgs' (whole-genome sequencing), 'amplicon' (targeted amplicon capture), 'hybrid' (hybridization capture). targets : GenomicArray Targeted genomic regions (for 'hybrid' and 'amplicon'). access : GenomicArray Sequencing-accessible regions of the reference genome (for 'hybrid' and 'wgs'). bp_per_bin : int Desired number of sequencing read nucleotide bases mapped to each bin. Returns ------- 2-tuple of 2-tuples: ((target depth, target avg. bin size), (antitarget depth, antitarget avg. bin size)) """ if method in ('amplicon', 'hybrid'): if targets is None: raise ValueError("Target regions are required for method %r " "but were not provided." % method) if not len(targets): raise ValueError("Target regions are required for method %r " "but were not provided." % method) # Closes over bp_per_bin def depth2binsize(depth, min_size, max_size): if depth: bin_size = int(round(bp_per_bin / depth)) if bin_size < min_size:"Limiting est. bin size %d to given min. %d", bin_size, min_size) bin_size = min_size elif bin_size > max_size:"Limiting est. bin size %d to given max. %d", bin_size, max_size) bin_size = max_size return bin_size samutil.ensure_bam_index(bam_fname) rc_table = samutil.idxstats(bam_fname, drop_unmapped=True, fasta=fasta) read_len = samutil.get_read_length(bam_fname, fasta=fasta)"Estimated read length %s", read_len) # Dispatch if method == 'amplicon': # From BAM index # rc_table = update_chrom_length(rc_table, targets) # tgt_depth = average_depth(rc_table, read_len) # By sampling tgt_depth = sample_region_cov(bam_fname, targets, fasta=fasta) anti_depth = None elif method == 'hybrid': tgt_depth, anti_depth = hybrid(rc_table, read_len, bam_fname, targets, access, fasta) elif method == 'wgs': if access is not None and len(access): rc_table = update_chrom_length(rc_table, access) tgt_depth = average_depth(rc_table, read_len) anti_depth = None # Clip bin sizes to specified ranges tgt_bin_size = depth2binsize(tgt_depth, target_min_size, target_max_size) anti_bin_size = depth2binsize(anti_depth, antitarget_min_size, antitarget_max_size) return ((tgt_depth, tgt_bin_size), (anti_depth, anti_bin_size))
[docs]def hybrid(rc_table, read_len, bam_fname, targets, access=None, fasta=None): """Hybrid capture sequencing.""" # Identify off-target regions if access is None: access = idxstats2ga(rc_table, bam_fname) # Verify BAM chromosome names match those in target BED compare_chrom_names(access, targets) antitargets = access.subtract(targets) # Only examine chromosomes present in all 2-3 input datasets rc_table, targets, antitargets = shared_chroms(rc_table, targets, antitargets) # Deal with targets target_depth = sample_region_cov(bam_fname, targets, fasta=fasta) # Antitargets: subtract captured reads from total target_length = region_size_by_chrom(targets)['length'] target_reads = (target_length * target_depth / read_len).values anti_table = update_chrom_length(rc_table, antitargets) anti_table = anti_table.assign(mapped=anti_table.mapped - target_reads) anti_depth = average_depth(anti_table, read_len) return target_depth, anti_depth
# ---
[docs]def average_depth(rc_table, read_length): """Estimate the average read depth across the genome. Returns ------- float Median of the per-chromosome mean read depths, weighted by chromosome size. """ mean_depths = read_length * rc_table.mapped / rc_table.length return weighted_median(mean_depths, rc_table.length)
[docs]def idxstats2ga(table, bam_fname): return GA(table.assign(start=0, end=table.length) .loc[:, ('chromosome', 'start', 'end')], meta_dict={'filename': bam_fname})
[docs]def sample_region_cov(bam_fname, regions, max_num=100, fasta=None): """Calculate read depth in a randomly sampled subset of regions.""" midsize_regions = sample_midsize_regions(regions, max_num) with tempfile.NamedTemporaryFile(suffix='.bed', mode='w+t') as f: tabio.write(regions.as_dataframe(midsize_regions), f, 'bed4') f.flush() table = coverage.bedcov(, bam_fname, 0, fasta) # Mean read depth across all sampled regions return table.basecount.sum() / (table.end - table.start).sum()
[docs]def sample_midsize_regions(regions, max_num): """Randomly select a subset of up to `max_num` regions.""" sizes = regions.end - regions.start lo_size, hi_size = np.percentile(sizes[sizes > 0], [25, 75]) midsize_regions =[(sizes >= lo_size) & (sizes <= hi_size)] if len(midsize_regions) > max_num: midsize_regions = midsize_regions.sample(max_num, random_state=0xA5EED) return midsize_regions
[docs]def shared_chroms(*tables): """Intersection of DataFrame .chromosome values.""" chroms = tables[0].chromosome.drop_duplicates() for tab in tables[1:]: if tab is not None: new_chroms = tab.chromosome.drop_duplicates() chroms = chroms[chroms.isin(new_chroms)] return [None if tab is None else tab[tab.chromosome.isin(chroms)] for tab in tables]
[docs]def update_chrom_length(rc_table, regions): if regions is not None and len(regions): chrom_sizes = region_size_by_chrom(regions) rc_table = rc_table.merge(chrom_sizes, on='chromosome', how='inner') rc_table['length'] = rc_table['length_y'] # ? rc_table = rc_table.drop(['length_x', 'length_y'], axis=1) return rc_table
[docs]def region_size_by_chrom(regions): chromgroups ='chromosome', sort=False) # sizes = chromgroups.apply(total_region_size) # XXX sizes = [total_region_size(g) for _key, g in chromgroups] return pd.DataFrame({'chromosome': regions.chromosome.drop_duplicates(), 'length': sizes})
[docs]def total_region_size(regions): """Aggregate area of all genomic ranges in `regions`.""" return (regions.end - regions.start).sum()