| jupytext | kernelspec | ||||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
|
|
This guide provides an introdution into how to use bioframe to perform genomic interval operations. For the full list of genomic interval operations, see the :ref:`API reference<API_ops>`.
The following modules are used in this guide:
import itertools
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
import bioframe as bf
import bioframe.vis
The core objects in bioframe are pandas DatFrames of genomic intervals, or BedFrames. These can either be defined directly with :py:class:`pandas.DataFrame`:
df1 = pd.DataFrame([
['chr1', 1, 5],
['chr1', 3, 8],
['chr1', 8, 10],
['chr1', 12, 14]],
columns=['chrom', 'start', 'end']
)
Or via functions in :mod:`bioframe.core.construction`, e.g.:
df2 = bioframe.from_any(
[['chr1', 4, 8],
['chr1', 10, 11]],
name_col='chrom')
Or ingested from datasets and databases with functions in :mod:`bioframe.io.fileops` and :mod:`bioframe.io.resources`.
BedFrames satisfy the following properties:
- chrom, start, end columns
- columns have valid dtypes (object/string/categorical, int/pd.Int64Dtype(), int/pd.Int64Dtype())
- for each interval, if any of chrom, start, end are null, then all are null
- all starts <= ends.
Whether a dataframe satisfies these properties can be checked with :func:`bioframe.core.checks.is_bedframe`:
bioframe.is_bedframe(df2)
See :mod:`bioframe.core.checks` for other functions that test properties of BedFrames and :ref:`Technical Notes<Technical_Notes>` for detailed definitions.
:func:`bioframe.core.construction.sanitize_bedframe` attempts to modfiy a DataFrame such that it satisfies bedFrame requirements.
:py:mod:`bioframe.vis` provides plotting utilities for intervals:
bf.vis.plot_intervals(df1, show_coords=True, xlim=(0,16))
plt.title('bedFrame1 intervals');
bf.vis.plot_intervals(df2, show_coords=True, xlim=(0,16), colors='lightpink')
plt.title('bedFrame2 intervals');
Calculating the overlap between two sets of genomic intervals is a crucial genomic interval operation.
Using :func:`bioframe.overlap`, we can see the two dataframes defined above, ``df1`` and ``df2``, contain two pairs of overlapping intervals:
overlapping_intervals = bf.overlap(df1, df2, how='inner', suffixes=('_1','_2'))
display(overlapping_intervals)
for i, reg_pair in overlapping_intervals.iterrows():
bf.vis.plot_intervals_arr(
starts = [reg_pair.start_1,reg_pair.start_2],
ends = [reg_pair.end_1,reg_pair.end_2],
colors = ['skyblue', 'lightpink'],
levels = [2,1],
xlim = (0,16),
show_coords = True)
plt.title(f'overlapping pair #{i}')
Note that we passed custom suffixes for the outputs (defaults are suffixes=("","_")),
as well as a custom overlap mode (how='inner'). The default overlap mode, how='left' returns each interval in df1 whether or not it overlaps an interval in df2.
overlapping_intervals = bf.overlap(df1, df2)
display(overlapping_intervals)
It is often useful to find overlapping intervals within a single set of genomic intervals. In `bioframe`, this is achieved with :func:`bioframe.cluster`. This function returns a DataFrame where subsets of overlapping intervals are assigned to the same group, reported in a new column.
To demonstrate the usage of :func:`bioframe.cluster`, we use the same ``df1`` as above:
df1 = pd.DataFrame([
['chr1', 1, 5],
['chr1', 3, 8],
['chr1', 8, 10],
['chr1', 12, 14],
],
columns=['chrom', 'start', 'end']
)
bf.vis.plot_intervals(df1, show_coords=True, xlim=(0,16))
Cluster returns a DataFrame where each interval is assigned to a group:
df_annotated = bf.cluster(df1, min_dist=0)
display(df_annotated)
bf.vis.plot_intervals(df_annotated, labels=df_annotated['cluster'], show_coords=True, xlim=(0,16))
Note that using min_dist=0 and min_dist=None give different results, as the latter only clusters overlapping intervals and not adjacent intervals:
df_annotated = bf.cluster(df1, min_dist=None)
display(df_annotated)
bf.vis.plot_intervals(df_annotated, labels=df_annotated['cluster'], show_coords=True, xlim=(0,16))
Extending the minimum distance to two (min_dist=2) makes all intervals part of the same cluster "0":
df_annotated = bf.cluster(df1, min_dist=2)
display(df_annotated)
bf.vis.plot_intervals(df_annotated, labels=df_annotated['cluster'], show_coords=True, xlim=(0,16))
Instead of returning cluster assignments, :func:`bioframe.merge` returns a new dataframe of merged genomic intervals. As with :func:`bioframe.cluster`, using ``min_dist=0`` and ``min_dist=None`` gives different results.
If ``min_dist=0``, this returns a dataframe of two intervals:
df_merged = bf.merge(df1, min_dist=0)
display(df_merged)
bf.vis.plot_intervals(df_merged, show_coords=True, xlim=(0,16))
If min_dist=None, this returns a dataframe of three intervals:
df_merged = bf.merge(df1, min_dist=None)
display(df_merged)
bf.vis.plot_intervals(df_merged, show_coords=True, xlim=(0,16))
In genomics, it is often useful not only to find features that overlap, but also features that are nearby along the genome. In bioframe, this is achieved using :func:`bioframe.closest`.
closest_intervals = bf.closest(df1, df2, suffixes=('_1','_2'))
display(closest_intervals)
for i, reg_pair in closest_intervals.iterrows():
bf.vis.plot_intervals_arr(
starts = [reg_pair.start_1,reg_pair.start_2],
ends = [reg_pair.end_1,reg_pair.end_2],
colors = ['skyblue', 'lightpink'],
levels = [2,1],
xlim = (0,16),
show_coords = True)
plt.title(f'closest pair #{i}')
By default, :func:`bioframe.closest` reports overlapping intervals. This can be modified by passing ``ignore_overlap=True``. Note the closest pair #2 and #3, which did not overlap, remain the same:
closest_intervals = bf.closest(df1, df2, ignore_overlaps=True, suffixes=('_1','_2'))
for i, reg_pair in closest_intervals.iterrows():
bf.vis.plot_intervals_arr(
starts = [reg_pair.start_1,reg_pair.start_2],
ends = [reg_pair.end_1,reg_pair.end_2],
colors = ['skyblue', 'lightpink'],
levels = [2,1],
xlim = (0,16),
show_coords = True)
plt.title(f'closest pair #{i}')
Closest intervals within a single DataFrame can be found simply by passing a single dataframe to :func:`bioframe.closest`. The number of closest intervals to report per query interval can be adjusted with ``k``.
bf.closest(df1, k=2)
Closest intervals upstream of the features in df1 can be found by ignoring downstream and overlaps.
Upstream/downstream direction is defined by genomic coordinates by default (smaller coordinate is upstream).
bf.closest(df1, df2,
ignore_overlaps=True,
ignore_downstream=True)
If the features in df1 have direction (e.g., genes have transcription direction), then the definition of upstream/downstream
direction can be changed to the direction of the features by `direction_col`:
df1["strand"] = np.where(np.random.rand(len(df1)) > 0.5, "+", "-")
bf.closest(df1, df2,
ignore_overlaps=True,
ignore_downstream=True,
direction_col='strand')
For two sets of genomic features, it is often useful to calculate the number of basepairs covered and the number of overlapping intervals. While these are fairly straightforward to compute from the output of :func:`bioframe.overlap` with :func:`pandas.groupby` and column renaming, since these are very frequently used, they are provided as core bioframe functions.
df1_coverage = bf.coverage(df1, df2)
display(df1_coverage)
df1_coverage_and_count = bf.count_overlaps(df1_coverage, df2)
display(df1_coverage_and_count)
This plot shows the coverage and number of overlaps for intervals in df1 by df2:
bf.vis.plot_intervals(
df1_coverage_and_count,
show_coords=True, xlim=(0,16),
labels = [f'{cov} bp, {n} intervals'
for cov, n in zip(df1_coverage_and_count.coverage, df1_coverage_and_count['count'])])
bf.vis.plot_intervals(df2, show_coords=True, xlim=(0,16), colors='lightpink')
Bioframe has two functions for computing differences between sets of intervals: at the level of basepairs and at the level of whole intervals.
Basepair-level subtraction is performed with :func:`bioframe.subtract`:
subtracted_intervals = bf.subtract(df1, df2)
display(subtracted_intervals)
bf.vis.plot_intervals(subtracted_intervals, show_coords=True, xlim=(0,16))
Interval-level differences are calculated with :func:`bioframe.setdiff`:
setdiff_intervals = bf.setdiff(df1, df2)
display(setdiff_intervals)
bf.vis.plot_intervals(setdiff_intervals, show_coords=True, xlim=(0,16))
:func:`bioframe.expand` enables quick resizing of intervals.
Expand supports additive resizing, with ``pad``.
Note that unless subsequently trimmed (with :func:`bioframe.trim`),
expanded intervals can have negative values:
expanded_intervals = bf.expand(df1, pad=2)
display(expanded_intervals)
bf.vis.plot_intervals(expanded_intervals, show_coords=True, xlim=(0,16))
Expand also supports multiplicative resizing, with ``scale``. Note that ``scale=0`` resizes all intervals to their midpoints:
expanded_intervals = bf.expand(df1, scale=0)
display(expanded_intervals)
bf.vis.plot_intervals(expanded_intervals, show_coords=True, xlim=(0,16))
Certain interval operations are often used relative to a set of reference intervals, whether those are chromosomes, scaffolds, or sub-intervals of either. Bioframe formalizes this with the concept of a `genomic view`, implemented as pandas dataframes, termed viewFrames, that satisfy the following:
- all requirements for bedFrames, including columns for ('chrom', 'start', 'end')
- it has an additional column, ``view_name_col``, with default 'name'
- entries in the view_name_col are unique
- intervals are non-overlapping
- it does not contain null values
Importanly a `genomic view` specifies a global coordinate axis, i.e. a conversion from a genomic coordinate system to a single axis. See :ref:`Technical Notes<Technical_Notes>` for more details.
Bioframe provides a function to assign intervals to corresponding intervals in a view, :func:`bioframe.assign_view`, a function to construct views from various input formats, :func:`bioframe.core.construction.make_viewframe`, and a function check that viewframe requirements are met, :func:`bioframe.core.checks.is_viewframe`.
The following genomic interval operations make use of views, though also have useful default behavior if no view is provided: :func:`bioframe.complement`, :func:`bioframe.trim`, :func:`bioframe.sort_bedframe`.
Equally important to finding which genomic features overlap is finding those that do not. :func:`bioframe.complement` returns a BedFrame of intervals not covered by any intervals in an input BedFrame.
Complments are defined relative to a `genomic view`. Here this is provided as a dictionary with a single chromosome of length 15:
df_complemented = bf.complement(df1, view_df={'chr1':15})
display(df_complemented)
bf.vis.plot_intervals(df_complemented, show_coords=True, xlim=(0,16), colors='lightpink')
If no view is provided, complement is calculated per unique chromosome in the input with right limits of np.iinfo(np.int64).max.
df_complemented = bf.complement(df1)
display(df_complemented)
Certain regions are often best avoided for genomic analyses. :func:`bioframe.trim` trims intervals to a specified view. Intervals falling outside of view regions have their filled with null values.
view_df = pd.DataFrame(
[
["chr1", 0, 4, "chr1p"],
["chr1", 5, 9, "chr1q"],
["chrX", 0, 50, "chrX"],
["chrM", 0, 10, "chrM"]],
columns=["chrom", "start", "end", "name"],
)
trimmed_intervals = bf.trim(df1, view_df)
display(trimmed_intervals)
Note that the last interval of df1 fell beyond 'chr1q' and is now null, and the last interval now ends at 9 instead of 10.
bf.vis.plot_intervals(trimmed_intervals, show_coords=True, xlim=(0,16))
If no view is provided, this function trims intervals at zero to avoid negative values.
If no view is provided, :func:`bioframe.sort_bedframe` sorts by ("chrom", "start", "end") columns:
df_unsorted = pd.DataFrame([
['chrM', 3, 8],
['chrM', 1, 5],
['chrX', 12, 14],
['chrX', 8, 10]],
columns=['chrom', 'start', 'end']
)
display( bf.sort_bedframe(df_unsorted) )
Views enable a specifying a sort order on a set of intervals. This flexibility is useful when the desired sorting order is non-lexicographical, e.g. with chrM after autosomes and chrX:
display( bf.sort_bedframe(df_unsorted, view_df) )
Since bioFrame operates directly with pandas DataFrames, all typical selection and slicing operations are directly relevant.
Bioframe also provides a function :func:`bioframe.select` that enables selecting interval subsets using UCSC string format:
display( bioframe.select(df_unsorted,'chrX:8-14') )
+++
Genomic analyses often deal with dataframes with inhomogeneously named columns. Bioframe offers a way to set the default column names that are most convenient for your analyses.
Default bedframe column names are stored in bioframe.core.specs_rc.
bf.core.specs._rc
If the dataframes we wish to work with have ['CHROMOSOME', 'LEFT', 'RIGHT'], we can either pass cols to operations in bioframe.ops:
df1_diff_colnames = pd.DataFrame([
['chr1', 1, 5],
['chr1', 3, 8]],
columns=['CHROMOSOME', 'LEFT', 'RIGHT']
)
df2_diff_colnames = pd.DataFrame([
['chr1', 4, 8],
['chr1', 10, 11]],
columns=['CHROMOSOME', 'LEFT', 'RIGHT']
)
bf.overlap(
df1_diff_colnames, df2_diff_colnames,
cols1=['CHROMOSOME', 'LEFT', 'RIGHT'],
cols2=['CHROMOSOME', 'LEFT', 'RIGHT'],
)
Or, we can update the default column names:
with bf.core.specs.update_default_colnames(['CHROMOSOME', 'LEFT', 'RIGHT']):
display(bf.overlap(
df1_diff_colnames, df2_diff_colnames,
))
# setting colnames back to default.
bf.core.specs.update_default_colnames(['chrom', 'start', 'end'])
bf.core.specs._rc