import collections
import os
import sqlite3
import shutil
import warnings
from gffutils import bins
from gffutils import helpers
from gffutils import constants
from gffutils import merge_criteria as mc
from gffutils.feature import Feature
from gffutils.exceptions import FeatureNotFoundError
def assign_child(parent, child):
"""
Helper for add_relation()
Sets 'Parent' attribute to parent['ID']
Parameters
----------
parent : Feature
Parent Feature
:param parent: parent Feature
child : Feature
Child Feature
Returns
-------
Child Feature
"""
child.attributes["Parent"] = parent["ID"]
return child
# Reusable constant for FeatureDB.merge()
no_children = tuple()
def _finalize_merge(feature, feature_children):
"""
Helper for FeatureDB.merge() to update source and assign children property
Parameters
----------
feature : Feature
feature to finalise
feature_children
list of children to assign
Returns
-------
feature, modified
"""
if len(feature_children) > 1:
feature.source = ",".join(set(child.source for child in feature_children))
feature.children = feature_children
else:
feature.children = no_children
return feature
[docs]
class FeatureDB(object):
# docstring to be filled in for methods that call out to
# helpers.make_query()
_method_doc = """
limit : string or tuple
Limit the results to a genomic region. If string, then of the form
"seqid:start-end"; if tuple, then (seqid, start, end).
strand : "-" | "+" | "."
Limit the results to one strand
featuretype : string or tuple
Limit the results to one or several featuretypes.
order_by : string or tuple
Order results by one or many fields; the string or tuple items must
be in: 'seqid', 'source', 'featuretype', 'start', 'end', 'score',
'strand', 'frame', 'attributes', 'extra'.
reverse : bool
Change sort order; only relevant if `order_by` is not None. By
default, results will be in ascending order, so use `reverse=True`
for descending.
completely_within : bool
If False (default), a single bp overlap with `limit` is sufficient
to return a feature; if True, then the feature must be completely
within `limit`. Only relevant when `limit` is not None.
"""
[docs]
def __init__(
self,
dbfn,
default_encoding="utf-8",
keep_order=False,
pragmas=constants.default_pragmas,
sort_attribute_values=False,
text_factory=str,
):
"""
Connect to a database created by :func:`gffutils.create_db`.
Parameters
----------
dbfn : str
Path to a database created by :func:`gffutils.create_db`.
text_factory : callable
Optionally set the way sqlite3 handles strings. Default is
str
default_encoding : str
When non-ASCII characters are encountered, assume they are in this
encoding.
keep_order : bool
If True, all features returned from this instance will have the
order of their attributes maintained. This can be turned on or off
database-wide by setting the `keep_order` attribute or with this
kwarg, or on a feature-by-feature basis by setting the `keep_order`
attribute of an individual feature.
Default is False, since this includes a sorting step that can get
time-consuming for many features.
sort_attributes_values : bool
If True, then in cases where there are multiple values for an
attribute then ensure they appear in the same order every time.
This is typically only used for testing, in cases where it is
important to have consistent ordering.
pragmas : dict
Dictionary of pragmas to use when connecting to the database. See
http://www.sqlite.org/pragma.html for the full list of
possibilities, and constants.default_pragmas for the defaults.
These can be changed later using the :meth:`FeatureDB.set_pragmas`
method.
Notes
-----
`dbfn` can also be a subclass of :class:`_DBCreator`, useful for
when :func:`gffutils.create_db` is provided the ``dbfn=":memory:"``
kwarg.
"""
# Since specifying the string ":memory:" will actually try to connect
# to a new, separate (and empty) db in memory, we can alternatively
# pass in a sqlite connection instance to use its existing, in-memory
# db.
from gffutils import create
if isinstance(dbfn, create._DBCreator):
self.conn = dbfn.conn
self.dbfn = dbfn.dbfn
elif isinstance(dbfn, sqlite3.Connection):
self.conn = dbfn
self.dbfn = dbfn
# otherwise assume dbfn is a string.
elif dbfn == ":memory:":
raise ValueError(
"cannot connect to memory db; please provide the connection"
)
else:
if not os.path.exists(dbfn):
raise ValueError("Database file %s does not exist" % dbfn)
self.dbfn = dbfn
self.conn = sqlite3.connect(self.dbfn)
if text_factory is not None:
self.conn.text_factory = text_factory
self.conn.row_factory = sqlite3.Row
self.default_encoding = default_encoding
self.keep_order = keep_order
self.sort_attribute_values = sort_attribute_values
c = self.conn.cursor()
# Load some meta info
# TODO: this is a good place to check for previous versions, and offer
# to upgrade...
c.execute(
"""
SELECT version, dialect FROM meta
"""
)
version, dialect = c.fetchone()
self.version = version
self.dialect = helpers._unjsonify(dialect)
# Load directives from db
c.execute(
"""
SELECT directive FROM directives
"""
)
self.directives = [directive[0] for directive in c if directive]
# Load autoincrements so that when we add new features, we can start
# autoincrementing from where we last left off (instead of from 1,
# which would cause name collisions)
c.execute(
"""
SELECT base, n FROM autoincrements
"""
)
self._autoincrements = collections.defaultdict(int, c)
self.set_pragmas(pragmas)
if not self._analyzed():
warnings.warn(
"It appears that this database has not had the ANALYZE "
"sqlite3 command run on it. Doing so can dramatically "
"speed up queries, and is done by default for databases "
"created with gffutils >0.8.7.1 (this database was "
"created with version %s) Consider calling the analyze() "
"method of this object." % self.version
)
[docs]
def set_pragmas(self, pragmas):
"""
Set pragmas for the current database connection.
Parameters
----------
pragmas : dict
Dictionary of pragmas; see constants.default_pragmas for a template
and http://www.sqlite.org/pragma.html for a full list.
"""
self.pragmas = pragmas
c = self.conn.cursor()
c.executescript(";\n".join(["PRAGMA %s=%s" % i for i in self.pragmas.items()]))
self.conn.commit()
def _feature_returner(self, **kwargs):
"""
Returns a feature, adding additional database-specific defaults
"""
kwargs.setdefault("dialect", self.dialect)
kwargs.setdefault("keep_order", self.keep_order)
kwargs.setdefault("sort_attribute_values", self.sort_attribute_values)
return Feature(**kwargs)
def _analyzed(self):
res = self.execute(
"""
SELECT name FROM sqlite_master WHERE type='table'
AND name='sqlite_stat1';
"""
)
return len(list(res)) == 1
[docs]
def schema(self):
"""
Returns the database schema as a string.
"""
c = self.conn.cursor()
c.execute(
"""
SELECT sql FROM sqlite_master
"""
)
results = []
for (i,) in c:
if i is not None:
results.append(i)
return "\n".join(results)
def __getitem__(self, key):
if isinstance(key, Feature):
key = key.id
c = self.conn.cursor()
try:
c.execute(constants._SELECT + " WHERE id = ?", (key,))
except sqlite3.ProgrammingError:
c.execute(
constants._SELECT + " WHERE id = ?",
(key.decode(self.default_encoding),),
)
results = c.fetchone()
# TODO: raise error if more than one key is found
if results is None:
raise FeatureNotFoundError(key)
return self._feature_returner(**results)
[docs]
def count_features_of_type(self, featuretype=None):
"""
Simple count of features.
Can be faster than "grep", and is faster than checking the length of
results from :meth:`gffutils.FeatureDB.features_of_type`.
Parameters
----------
featuretype : string
Feature type (e.g., "gene") to count. If None, then count *all*
features in the database.
Returns
-------
The number of features of this type, as an integer
"""
c = self.conn.cursor()
if featuretype is not None:
c.execute(
"""
SELECT count() FROM features
WHERE featuretype = ?
""",
(featuretype,),
)
else:
c.execute(
"""
SELECT count() FROM features
"""
)
results = c.fetchone()
if results is not None:
results = results[0]
return results
[docs]
def features_of_type(
self,
featuretype,
limit=None,
strand=None,
order_by=None,
reverse=False,
completely_within=False,
):
"""
Returns an iterator of :class:`gffutils.Feature` objects.
Parameters
----------
{_method_doc}
"""
query, args = helpers.make_query(
args=[],
limit=limit,
featuretype=featuretype,
order_by=order_by,
reverse=reverse,
strand=strand,
completely_within=completely_within,
)
for i in self._execute(query, args):
yield self._feature_returner(**i)
# TODO: convert this to a syntax similar to itertools.groupby
[docs]
def iter_by_parent_childs(
self,
featuretype="gene",
level=None,
order_by=None,
reverse=False,
completely_within=False,
):
"""
For each parent of type `featuretype`, yield a list L of that parent
and all of its children (`[parent] + list(children)`). The parent will
always be L[0].
This is useful for "sanitizing" a GFF file for downstream tools.
Additional kwargs are passed to :meth:`FeatureDB.children`, and will
therefore only affect items L[1:] in each yielded list.
"""
# Get all the parent records of the requested feature type
parent_recs = self.all_features(featuretype=featuretype)
# Return a generator of these parent records and their
# children
for parent_rec in parent_recs:
unit_records = [parent_rec] + list(self.children(parent_rec.id))
yield unit_records
[docs]
def all_features(
self,
limit=None,
strand=None,
featuretype=None,
order_by=None,
reverse=False,
completely_within=False,
):
"""
Iterate through the entire database.
Returns
-------
A generator object that yields :class:`Feature` objects.
Parameters
----------
{_method_doc}
"""
query, args = helpers.make_query(
args=[],
limit=limit,
strand=strand,
featuretype=featuretype,
order_by=order_by,
reverse=reverse,
completely_within=completely_within,
)
for i in self._execute(query, args):
yield self._feature_returner(**i)
[docs]
def featuretypes(self):
"""
Iterate over feature types found in the database.
Returns
-------
A generator object that yields featuretypes (as strings)
"""
c = self.conn.cursor()
c.execute(
"""
SELECT DISTINCT featuretype from features
"""
)
for (i,) in c:
yield i
def _relation(
self,
id,
join_on,
join_to,
level=None,
featuretype=None,
order_by=None,
reverse=False,
completely_within=False,
limit=None,
):
# The following docstring will be included in the parents() and
# children() docstrings to maintain consistency, since they both
# delegate to this method.
"""
Parameters
----------
id : string or a Feature object
level : None or int
If `level=None` (default), then return all children regardless
of level. If `level` is an integer, then constrain to just that
level.
{_method_doc}
Returns
-------
A generator object that yields :class:`Feature` objects.
"""
if isinstance(id, Feature):
id = id.id
other = """
JOIN relations
ON relations.{join_on} = features.id
WHERE relations.{join_to} = ?
""".format(
**locals()
)
args = [id]
level_clause = ""
if level is not None:
level_clause = "relations.level = ?"
args.append(level)
query, args = helpers.make_query(
args=args,
other=other,
extra=level_clause,
featuretype=featuretype,
order_by=order_by,
reverse=reverse,
limit=limit,
completely_within=completely_within,
)
# modify _SELECT so that only unique results are returned
query = query.replace("SELECT", "SELECT DISTINCT")
for i in self._execute(query, args):
yield self._feature_returner(**i)
[docs]
def children(
self,
id,
level=None,
featuretype=None,
order_by=None,
reverse=False,
limit=None,
completely_within=False,
):
"""
Return children of feature `id`.
{_relation_docstring}
"""
return self._relation(
id,
join_on="child",
join_to="parent",
level=level,
featuretype=featuretype,
order_by=order_by,
reverse=reverse,
limit=limit,
completely_within=completely_within,
)
[docs]
def parents(
self,
id,
level=None,
featuretype=None,
order_by=None,
reverse=False,
completely_within=False,
limit=None,
):
"""
Return parents of feature `id`.
{_relation_docstring}
"""
return self._relation(
id,
join_on="parent",
join_to="child",
level=level,
featuretype=featuretype,
order_by=order_by,
reverse=reverse,
limit=limit,
completely_within=completely_within,
)
def _execute(self, query, args):
self._last_query = query
self._last_args = args
c = self.conn.cursor()
c.execute(query, tuple(args))
return c
[docs]
def execute(self, query):
"""
Execute arbitrary queries on the db.
.. seealso::
:class:`FeatureDB.schema` may be helpful when writing your own
queries.
Parameters
----------
query : str
Query to execute -- trailing ";" optional.
Returns
-------
A sqlite3.Cursor object that can be iterated over.
"""
c = self.conn.cursor()
return c.execute(query)
def analyze(self):
"""
Runs the sqlite ANALYZE command to potentially speed up queries
dramatically.
"""
self.execute("ANALYZE features")
self.conn.commit()
[docs]
def region(
self,
region=None,
seqid=None,
start=None,
end=None,
strand=None,
featuretype=None,
completely_within=False,
):
"""
Return features within specified genomic coordinates.
Specifying genomic coordinates can be done in a flexible manner
Parameters
----------
region : string, tuple, or Feature instance
If string, then of the form "seqid:start-end". If tuple, then
(seqid, start, end). If :class:`Feature`, then use the features
seqid, start, and end values.
This argument is mutually exclusive with start/end/seqid.
*Note*: By design, even if a feature is provided, its strand will
be ignored. If you want to restrict the output by strand, use the
separate `strand` kwarg.
strand : + | - | . | None
If `strand` is provided, then only those features exactly matching
`strand` will be returned. So `strand='.'` will only return
unstranded features. Default is `strand=None` which does not
restrict by strand.
seqid, start, end, strand
Mutually exclusive with `region`. These kwargs can be used to
approximate slice notation; see "Details" section below.
featuretype : None, string, or iterable
If not None, then restrict output. If string, then only report
that feature type. If iterable, then report all featuretypes in
the iterable.
completely_within : bool
By default (`completely_within=False`), returns features that
partially or completely overlap `region`. If
`completely_within=True`, features that are completely within
`region` will be returned.
Notes
-----
The meaning of `seqid`, `start`, and `end` is interpreted as follows:
====== ====== ===== ======================================
seqid start end meaning
====== ====== ===== ======================================
str int int equivalent to `region` kwarg
None int int features from all chroms within coords
str None int equivalent to [:end] slice notation
str int None equivalent to [start:] slice notation
None None None equivalent to FeatureDB.all_features()
====== ====== ===== ======================================
If performance is a concern, use `completely_within=True`. This allows
the query to be optimized by only looking for features that fall in the
precise genomic bin (same strategy as UCSC Genome Browser and
BEDTools). Otherwise all features' start/stop coords need to be
searched to see if they partially overlap the region of interest.
Examples
--------
- `region(seqid="chr1", start=1000)` returns all features on chr1 that
start or extend past position 1000
- `region(seqid="chr1", start=1000, completely_within=True)` returns
all features on chr1 that start past position 1000.
- `region("chr1:1-100", strand="+", completely_within=True)` returns
only plus-strand features that completely fall within positions 1 to
100 on chr1.
Returns
-------
A generator object that yields :class:`Feature` objects.
"""
# Argument handling.
if region is not None:
if (seqid is not None) or (start is not None) or (end is not None):
raise ValueError(
"If region is supplied, do not supply seqid, "
"start, or end as separate kwargs"
)
if isinstance(region, str):
toks = region.split(":")
if len(toks) == 1:
seqid = toks[0]
start, end = None, None
else:
seqid, coords = toks[:2]
if len(toks) == 3:
strand = toks[2]
start, end = coords.split("-")
elif isinstance(region, Feature):
seqid = region.seqid
start = region.start
end = region.end
strand = region.strand
# otherwise assume it's a tuple
else:
seqid, start, end = region[:3]
# e.g.,
# completely_within=True..... start >= {start} AND end <= {end}
# completely_within=False.... start < {end} AND end > {start}
if completely_within:
start_op = ">="
end_op = "<="
else:
start_op = "<"
end_op = ">"
end, start = start, end
args = []
position_clause = []
if seqid is not None:
position_clause.append("seqid = ?")
args.append(seqid)
if start is not None:
start = int(start)
if end is not None:
end = int(end)
# See #129
if start and end and not completely_within:
position_clause.append(
"""(
({region_start} <= start AND {region_end} >= start) OR
({region_start} >= start AND {region_end} <= end) OR
({region_start} <= end AND {region_end} >= end)
)""".format(
region_start=start, region_end=end
)
)
else:
if start:
position_clause.append("start %s ?" % start_op)
args.append(start)
if end:
position_clause.append("end %s ?" % end_op)
args.append(end)
position_clause = " AND ".join(position_clause)
# Only use bins if we have defined boundaries and completely_within is
# True. Otherwise you can't know how far away a feature stretches
# (which means bins are not computable ahead of time)
_bin_clause = ""
if (start is not None) and (end is not None) and completely_within:
if start <= bins.MAX_CHROM_SIZE and end <= bins.MAX_CHROM_SIZE:
_bins = list(bins.bins(start, end, one=False))
# See issue #45
if len(_bins) < 900:
_bin_clause = " or ".join(["bin = ?" for _ in _bins])
_bin_clause = "AND ( %s )" % _bin_clause
args += _bins
query = " ".join([constants._SELECT, "WHERE ", position_clause, _bin_clause])
# Add the featuretype clause
if featuretype is not None:
if isinstance(featuretype, str):
featuretype = [featuretype]
feature_clause = " or ".join(["featuretype = ?" for _ in featuretype])
query += " AND (%s) " % feature_clause
args.extend(featuretype)
if strand is not None:
strand_clause = " and strand = ? "
query += strand_clause
args.append(strand)
c = self.conn.cursor()
self._last_query = query
self._last_args = args
self._context = {
"start": start,
"end": end,
"seqid": seqid,
"region": region,
}
c.execute(query, tuple(args))
for i in c:
yield self._feature_returner(**i)
[docs]
def interfeatures(
self,
features,
new_featuretype=None,
merge_attributes=True,
numeric_sort=False,
dialect=None,
attribute_func=None,
update_attributes=None,
):
"""
Construct new features representing the space between features.
For example, if `features` is a list of exons, then this method will
return the introns. If `features` is a list of genes, then this method
will return the intergenic regions.
Providing N features will return N - 1 new features.
This method purposefully does *not* do any merging or sorting of
coordinates. So nested or overlapping features may not behave as you
might expect. You may want to use :meth:`FeatureDB.merge` first, and
when selecting features use the `order_by` kwarg, e.g.,
`db.features_of_type('gene', order_by=('seqid', 'start'))`.
Parameters
----------
features : iterable of :class:`feature.Feature` instances
Sorted, merged iterable
new_featuretype : string or None
The new features will all be of this type, or, if None (default)
then the featuretypes will be constructed from the neighboring
features, e.g., `inter_exon_exon`.
merge_attributes : bool
If True, new features' attributes will be a merge of the neighboring
features' attributes. This is useful if you have provided a list of
exons; the introns will then retain the transcript and/or gene
parents as a single item. Otherwise, if False, the attribute will
be a comma-separated list of values, potentially listing the same
gene ID twice.
numeric_sort : bool
If True, then merged attributes that can be cast to float will be
sorted by their numeric values (but will still be returned as
string). This is useful, for example, when creating introns between
exons and the exons have exon_number attributes as an integer.
Using numeric_sort=True will ensure that the returned exons have
merged exon_number attribute of ['9', '10'] (numerically sorted)
rather than ['10', '9'] (alphabetically sorted).
attribute_func : callable or None
If None, then nothing special is done to the attributes. If
callable, then the callable accepts two attribute dictionaries and
returns a single attribute dictionary. If `merge_attributes` is
True, then `attribute_func` is called before `merge_attributes`.
This could be useful for manually managing IDs for the new
features.
update_attributes : dict
After attributes have been modified and merged, this dictionary can
be used to replace parts of the attributes dictionary.
Returns
-------
A generator that yields :class:`Feature` objects
"""
def _init_interfeature(f):
"""
Used to initialize a new interfeature that is ready to be updated
in-place.
"""
keys = [
"id",
"seqid",
"source",
"featuretype",
"start",
"end",
"score",
"strand",
"frame",
"attributes",
"bin",
]
d = dict(zip(keys, f.astuple()))
d["source"] = "gffutils_derived"
return d
def _prep_for_yield(d):
"""
Finalize the interfeature by adjusting coords, recalculating the
bin, and creating a feature using self._feature_returner.
If start is greater than stop (which happens when trying to get
interfeatures for overlapping features), then return None.
"""
d["start"] += 1
d["end"] -= 1
new_bin = bins.bins(d["start"], d["end"], one=True)
d["bin"] = new_bin
if d["start"] > d["end"]:
return None
new_feature = self._feature_returner(**d)
# concat list of ID to create uniq IDs because feature with
# multiple values for their ID are no longer permitted since v0.11
if "ID" in new_feature.attributes and len(new_feature.attributes["ID"]) > 1:
new_id = "-".join(new_feature.attributes["ID"])
new_feature.attributes["ID"] = [new_id]
return new_feature
# If not provided, use a no-op function instead.
if not attribute_func:
def attribute_func(a):
return a
for i, f in enumerate(features):
# First feature: initialize an interfeature and continue to the next.
if i == 0:
interfeature = _init_interfeature(f)
last_feature = f
nfeatures = 1
continue
# Yield the last interfeature (if we saw at least 2 features) and
# start a new interfeature on this chrom.
if f.chrom != last_feature.chrom:
if nfeatures > 1:
new_feature = _prep_for_yield(interfeature)
if new_feature:
yield new_feature
interfeature = _init_interfeature(f)
last_feature = f
nfeatures = 1
continue
# Otherwise, we've already seen a feature on this chrom so
# this is the second.
nfeatures += 1
# Adjust the interfeature dict in-place with coords...
interfeature["start"] = last_feature.stop
interfeature["end"] = f.start
# ...featuretype
if new_featuretype is None:
interfeature["featuretype"] = "inter_%s_%s" % (
last_feature.featuretype,
f.featuretype,
)
else:
interfeature["featuretype"] = new_featuretype
# ...strand
if last_feature.strand != f.strand:
interfeature["strand"] = "."
else:
interfeature["strand"] = f.strand
# and attributes
if merge_attributes:
new_attributes = helpers.merge_attributes(
attribute_func(last_feature.attributes),
attribute_func(f.attributes),
numeric_sort=numeric_sort,
)
else:
new_attributes = {}
if update_attributes:
new_attributes.update(update_attributes)
interfeature["attributes"] = new_attributes
# Ready to yield
new_feature = _prep_for_yield(interfeature)
if new_feature:
yield new_feature
nfeatures = 1
last_feature = f
[docs]
def delete(self, features, make_backup=True, **kwargs):
"""
Delete features from database.
features : str, iterable, FeatureDB instance
If FeatureDB, all features will be used. If string, assume it's the
ID of the feature to remove. Otherwise, assume it's an iterable of
Feature objects. The classes in gffutils.iterators may be helpful
in this case.
make_backup : bool
If True, and the database you're about to update is a file on disk,
makes a copy of the existing database and saves it with a .bak
extension.
Returns
-------
FeatureDB object, with features deleted.
"""
if make_backup:
if isinstance(self.dbfn, str):
shutil.copy2(self.dbfn, self.dbfn + ".bak")
c = self.conn.cursor()
query1 = """
DELETE FROM features WHERE id = ?
"""
query2 = """
DELETE FROM relations WHERE parent = ? OR child = ?
"""
if isinstance(features, FeatureDB):
features = features.all_features()
if isinstance(features, str):
features = [features]
if isinstance(features, Feature):
features = [features]
for feature in features:
if isinstance(feature, str):
_id = feature
else:
_id = feature.id
c.execute(query1, (_id,))
c.execute(query2, (_id, _id))
self.conn.commit()
return self
[docs]
def update(self, data, make_backup=True, **kwargs):
"""
Update the on-disk database with features in `data`.
WARNING: If you used any non-default kwargs for gffutils.create_db when
creating the database in the first place (especially
`disable_infer_transcripts` or `disable_infer_genes`) then you should
use those same arguments here. The exception is the `force` argument
though -- in some cases including that can truncate the database.
WARNING: If you are creating features from the database and writing
immediately back to the database, you could experience deadlocks. See
the help for `create_introns` for some different options for avoiding
this.
The returned object is the same FeatureDB, but since it is pointing to
the same database and that has been just updated, the new features can
now be accessed.
Parameters
----------
data : str, iterable, FeatureDB instance
If FeatureDB, all data will be used. If string, assume it's
a filename of a GFF or GTF file. Otherwise, assume it's an
iterable of Feature objects. The classes in gffutils.iterators may
be helpful in this case.
make_backup : bool
If True, and the database you're about to update is a file on disk,
makes a copy of the existing database and saves it with a .bak
extension.
kwargs :
Additional kwargs are passed to gffutils.create_db; see the help
for that function for details.
Returns
-------
Returns self but with the underlying database updated.
"""
from gffutils import create
from gffutils import iterators
if make_backup:
if isinstance(self.dbfn, str):
shutil.copy2(self.dbfn, self.dbfn + ".bak")
# get iterator-specific kwargs
_iterator_kwargs = {}
for k, v in kwargs.items():
if k in constants._iterator_kwargs:
_iterator_kwargs[k] = v
# Handle all sorts of input
data = iterators.DataIterator(data, **_iterator_kwargs)
if not data._peek:
return self
kwargs["_autoincrements"] = self._autoincrements
if self.dialect["fmt"] == "gtf":
if "id_spec" not in kwargs:
kwargs["id_spec"] = {"gene": "gene_id", "transcript": "transcript_id"}
db = create._GTFDBCreator(
data=data, dbfn=self.dbfn, dialect=self.dialect, **kwargs
)
elif self.dialect["fmt"] == "gff3":
if "id_spec" not in kwargs:
kwargs["id_spec"] = "ID"
db = create._GFFDBCreator(
data=data, dbfn=self.dbfn, dialect=self.dialect, **kwargs
)
else:
raise ValueError
db._populate_from_lines(data)
db._update_relations()
# Note that the autoincrements gets updated here
db._finalize()
# Read it back in directly from the stored autoincrements table
self._autoincrements.update(db._autoincrements)
return self
[docs]
def add_relation(self, parent, child, level, parent_func=None, child_func=None):
"""
Manually add relations to the database.
Parameters
----------
parent : str or Feature instance
Parent feature to add.
child : str or Feature instance
Child feature to add
level : int
Level of the relation. For example, if parent is a gene and child
is an mRNA, then you might want level to be 1. But if child is an
exon, then level would be 2.
parent_func, child_func : callable
These optional functions control how attributes are updated in the
database. They both have the signature `func(parent, child)` and
must return a [possibly modified] Feature instance. For example,
we could add the child's database id as the "child" attribute in
the parent::
def parent_func(parent, child):
parent.attributes['child'] = child.id
and add the parent's "gene_id" as the child's "Parent" attribute::
def child_func(parent, child):
child.attributes['Parent'] = parent['gene_id']
Returns
-------
FeatureDB object with new relations added.
"""
if isinstance(parent, str):
parent = self[parent]
if isinstance(child, str):
child = self[child]
c = self.conn.cursor()
c.execute(
"""
INSERT INTO relations (parent, child, level)
VALUES (?, ?, ?)""",
(parent.id, child.id, level),
)
if parent_func is not None:
parent = parent_func(parent, child)
self._update(parent, c)
if child_func is not None:
child = child_func(parent, child)
self._update(child, c)
self.conn.commit()
return self
def _update(self, feature, cursor):
values = [list(feature.astuple()) + [feature.id]]
cursor.execute(constants._UPDATE, *tuple(values))
def _insert(self, feature, cursor):
"""
Insert a feature into the database.
"""
try:
cursor.execute(constants._INSERT, feature.astuple())
except sqlite3.ProgrammingError:
cursor.execute(constants._INSERT, feature.astuple(self.default_encoding))
[docs]
def create_introns(
self,
exon_featuretype="exon",
grandparent_featuretype="gene",
parent_featuretype=None,
new_featuretype="intron",
merge_attributes=True,
numeric_sort=False,
):
"""
Create introns from existing annotations.
Parameters
----------
exon_featuretype : string
Feature type to use in order to infer introns. Typically `"exon"`.
grandparent_featuretype : string
If `grandparent_featuretype` is not None, then group exons by
children of this featuretype. If `granparent_featuretype` is
"gene" (default), then introns will be created for all first-level
children of genes. This may include mRNA, rRNA, ncRNA, etc. If
you only want to infer introns from one of these featuretypes
(e.g., mRNA), then use the `parent_featuretype` kwarg which is
mutually exclusive with `grandparent_featuretype`.
parent_featuretype : string
If `parent_featuretype` is not None, then only use this featuretype
to infer introns. Use this if you only want a subset of
featuretypes to have introns (e.g., "mRNA" only, and not ncRNA or
rRNA). Mutually exclusive with `grandparent_featuretype`.
new_featuretype : string
Feature type to use for the inferred introns; default is
`"intron"`.
merge_attributes : bool
Whether or not to merge attributes from all exons. If False then no
attributes will be created for the introns.
numeric_sort : bool
If True, then merged attributes that can be cast to float will be
sorted by their numeric values (but will still be returned as
string). This is useful, for example, when creating introns between
exons and the exons have exon_number attributes as an integer.
Using numeric_sort=True will ensure that the returned exons have
merged exon_number attribute of ['9', '10'] (numerically sorted)
rather than ['10', '9'] (alphabetically sorted).
Returns
-------
A generator object that yields :class:`Feature` objects representing
new introns
Notes
-----
The returned generator can be passed directly to the
:meth:`FeatureDB.update` method to permanently add them to the
database. However, this needs to be done carefully to avoid deadlocks
from simultaneous reading/writing.
When using `update()` you should also use the same keyword arguments
used to create the db in the first place (with the exception of `force`).
Here are three options for getting the introns back into the database,
depending on the circumstances.
**OPTION 1: Create list of introns.**
Consume the `create_introns()` generator completely before writing to
the database. If you have sufficient memory, this is the easiest
option::
db.update(list(db.create_introns(**intron_kwargs)), **create_kwargs)
**OPTION 2: Use `WAL <https://sqlite.org/wal.html>`__**
The WAL pragma enables simultaneous read/write. WARNING: this does not
work if the database is on a networked filesystem, like those used on
many HPC clusters.
::
db.set_pragmas({"journal_mode": "WAL"})
db.update(db.create_introns(**intron_kwargs), **create_kwargs)
**OPTION 3: Write to intermediate file.**
Use this if you are memory limited and using a networked filesystem::
with open('tmp.gtf', 'w') as fout:
for intron in db.create_introns(**intron_kwargs):
fout.write(str(intron) + "\n")
db.update(gffutils.DataIterator('tmp.gtf'), **create_kwargs)
"""
if (grandparent_featuretype and parent_featuretype) or (
grandparent_featuretype is None and parent_featuretype is None
):
raise ValueError(
"exactly one of `grandparent_featuretype` or "
"`parent_featuretype` should be provided"
)
if grandparent_featuretype:
def child_gen():
for gene in self.features_of_type(grandparent_featuretype):
for child in self.children(gene, level=1):
yield child
elif parent_featuretype:
def child_gen():
for child in self.features_of_type(parent_featuretype):
yield child
for child in child_gen():
exons = self.children(
child, level=1, featuretype=exon_featuretype, order_by="start"
)
for intron in self.interfeatures(
exons,
new_featuretype=new_featuretype,
merge_attributes=merge_attributes,
numeric_sort=numeric_sort,
dialect=self.dialect,
):
yield intron
def create_splice_sites(
self,
exon_featuretype="exon",
grandparent_featuretype="gene",
parent_featuretype=None,
merge_attributes=True,
numeric_sort=False,
):
"""
Create splice sites from existing annotations.
Parameters
----------
exon_featuretype : string
Feature type to use in order to infer splice sites. Typically `"exon"`.
grandparent_featuretype : string
If `grandparent_featuretype` is not None, then group exons by
children of this featuretype. If `granparent_featuretype` is
"gene" (default), then splice sites will be created for all first-level
children of genes. This may include mRNA, rRNA, ncRNA, etc. If
you only want to infer splice sites from one of these featuretypes
(e.g., mRNA), then use the `parent_featuretype` kwarg which is
mutually exclusive with `grandparent_featuretype`.
parent_featuretype : string
If `parent_featuretype` is not None, then only use this featuretype
to infer splice sites. Use this if you only want a subset of
featuretypes to have splice sites (e.g., "mRNA" only, and not ncRNA or
rRNA). Mutually exclusive with `grandparent_featuretype`.
merge_attributes : bool
Whether or not to merge attributes from all exons. If False then no
attributes will be created for the splice sites.
numeric_sort : bool
If True, then merged attributes that can be cast to float will be
sorted by their numeric values (but will still be returned as
string). This is useful, for example, when creating splice sites between
exons and the exons have exon_number attributes as an integer.
Using numeric_sort=True will ensure that the returned exons have
merged exon_number attribute of ['9', '10'] (numerically sorted)
rather than ['10', '9'] (alphabetically sorted).
Returns
-------
A generator object that yields :class:`Feature` objects representing
new splice sites
Notes
-----
The returned generator can be passed directly to the
:meth:`FeatureDB.update` method to permanently add them to the
database, e.g., ::
db.update(db.create_splice sites())
"""
if (grandparent_featuretype and parent_featuretype) or (
grandparent_featuretype is None and parent_featuretype is None
):
raise ValueError(
"exactly one of `grandparent_featuretype` or "
"`parent_featuretype` should be provided"
)
if grandparent_featuretype:
def child_gen():
for gene in self.features_of_type(grandparent_featuretype):
for child in self.children(gene, level=1):
yield child
elif parent_featuretype:
def child_gen():
for child in self.features_of_type(parent_featuretype):
yield child
# Two splice features need to be created for each interleave
for side in ["left", "right"]:
for child in child_gen():
exons = self.children(
child, level=1, featuretype=exon_featuretype, order_by="start"
)
# get strand
strand = child.strand
new_featuretype = "splice_site"
if side == "left":
if strand == "+":
new_featuretype = "five_prime_cis_splice_site"
elif strand == "-":
new_featuretype = "three_prime_cis_splice_site"
if side == "right":
if strand == "+":
new_featuretype = "three_prime_cis_splice_site"
elif strand == "-":
new_featuretype = "five_prime_cis_splice_site"
for splice_site in self.interfeatures(
exons,
new_featuretype=new_featuretype,
merge_attributes=merge_attributes,
numeric_sort=numeric_sort,
dialect=self.dialect,
):
if side == "left":
splice_site.end = splice_site.start + 1
if side == "right":
splice_site.start = splice_site.end - 1
# make ID uniq by adding suffix
splice_site.attributes["ID"] = [
new_featuretype + "_" + splice_site.attributes["ID"][0]
]
yield splice_site
def _old_merge(self, features, ignore_strand=False):
"""
DEPRECATED, only retained here for backwards compatibility. Please use
merge().
Merge overlapping features together.
Parameters
----------
features : iterator of Feature instances
ignore_strand : bool
If True, features on multiple strands will be merged, and the final
strand will be set to '.'. Otherwise, ValueError will be raised if
trying to merge features on differnt strands.
Returns
-------
A generator object that yields :class:`Feature` objects representing
the newly merged features.
"""
# Consume iterator up front...
features = list(features)
if len(features) == 0:
raise StopIteration
# Either set all strands to '+' or check for strand-consistency.
if ignore_strand:
strand = "."
else:
strands = [i.strand for i in features]
if len(set(strands)) > 1:
raise ValueError(
"Specify ignore_strand=True to force merging " "of multiple strands"
)
strand = strands[0]
# Sanity check to make sure all features are from the same chromosome.
chroms = [i.chrom for i in features]
if len(set(chroms)) > 1:
raise NotImplementedError("Merging multiple chromosomes not " "implemented")
chrom = chroms[0]
# To start, we create a merged feature of just the first feature.
current_merged_start = features[0].start
current_merged_stop = features[0].stop
# We don't need to check the first one, so start at feature #2.
for feature in features[1:]:
# Does this feature start within the currently merged feature?...
if feature.start <= current_merged_stop + 1:
# ...It starts within, so leave current_merged_start where it
# is. Does it extend any farther?
if feature.stop >= current_merged_stop:
# Extends further, so set a new stop position
current_merged_stop = feature.stop
else:
# If feature.stop < current_merged_stop, it's completely
# within the previous feature. Nothing more to do.
continue
else:
# The start position is outside the merged feature, so we're
# done with the current merged feature. Prepare for output...
merged_feature = dict(
seqid=feature.chrom,
source=".",
featuretype=feature.featuretype,
start=current_merged_start,
end=current_merged_stop,
score=".",
strand=strand,
frame=".",
attributes="",
)
yield self._feature_returner(**merged_feature)
# and we start a new one, initializing with this feature's
# start and stop.
current_merged_start = feature.start
current_merged_stop = feature.stop
# need to yield the last one.
if len(features) == 1:
feature = features[0]
merged_feature = dict(
seqid=feature.chrom,
source=".",
featuretype=feature.featuretype,
start=current_merged_start,
end=current_merged_stop,
score=".",
strand=strand,
frame=".",
attributes="",
)
yield self._feature_returner(**merged_feature)
[docs]
def merge(
self,
features,
merge_criteria=(mc.seqid, mc.overlap_end_inclusive, mc.strand, mc.feature_type),
multiline=False,
):
"""
Merge features matching criteria together
Returned Features have a special property called 'children' that is
a list of the component features. This only exists for the lifetime of
the Feature instance.
Parameters
----------
features : iterable
Iterable of Feature instances to merge
merge_criteria : list
List of merge criteria callbacks. All must evaluate to True in
order for a feature to be merged. See notes below on callback
signature.
multiline : bool
True to emit multiple features with the same ID attribute, False
otherwise.
Returns
-------
Generator yielding merged Features
Notes
-----
See the `gffutils.merge_criteria` module (imported here as `mc`) for
existing callback functions. For writing custom callbacks, functions
must have the following signature::
callback(
acc: gffutils.Feature,
cur: gffutils.Feature,
components: [gffutils.Feature]
) -> bool
Where:
- `acc`: current accumulated feature
- `cur`: candidate feature to merge
- `components`: list of features that compose acc
The function should return True to merge `cur` into `acc`, False to set
`cur` to `acc` (that is, start a new merged feature).
If merge criteria allows different feature types then the merged
features' feature types should have their featuretype property
reassigned to a more specific ontology value.
"""
if not isinstance(merge_criteria, list):
try:
merge_criteria = list(merge_criteria)
except TypeError:
merge_criteria = [merge_criteria]
# To start, we create a merged feature of just the first feature.
features = iter(features)
last_id = None
current_merged = None
feature_children = []
for feature in features:
if current_merged is None:
if all(
criteria(feature, feature, feature_children)
for criteria in merge_criteria
):
current_merged = feature
feature_children = [feature]
else:
yield _finalize_merge(feature, no_children)
last_id = None
continue
if (
len(feature_children) == 0
): # current_merged is last feature and unchecked
if all(
criteria(current_merged, current_merged, feature_children)
for criteria in merge_criteria
):
feature_children.append(current_merged)
else:
yield _finalize_merge(current_merged, no_children)
current_merged = feature
last_id = None
continue
if all(
criteria(current_merged, feature, feature_children)
for criteria in merge_criteria
):
# Criteria satisfied, merge
# TODO Test multiline records and iron out the following code
# if multiline and (feature.start > current_merged.end + 1 or feature.end + 1 < current_merged.start):
# # Feature is possibly multiline (discontiguous), keep ID but start new record
# yield _finalize_merge(current_merged, feature_children)
# current_merged = feature
# feature_children = [feature]
if len(feature_children) == 1:
# Current merged is only child and merge is going to occur, make copy
current_merged = vars(current_merged).copy()
del current_merged["attributes"]
del current_merged["extra"]
del current_merged["dialect"]
del current_merged["keep_order"]
del current_merged["sort_attribute_values"]
current_merged = self._feature_returner(**current_merged)
if not last_id:
# Generate unique ID for new Feature
self._autoincrements[current_merged.featuretype] += 1
last_id = (
current_merged.featuretype
+ "_"
+ str(self._autoincrements[current_merged.featuretype])
)
current_merged["ID"] = last_id
current_merged.id = last_id
feature_children.append(feature)
# Set mismatched properties to ambiguous values
if feature.seqid not in current_merged.seqid.split(","):
current_merged.seqid += "," + feature.seqid
if feature.strand != current_merged.strand:
current_merged.strand = "."
if feature.frame != current_merged.frame:
current_merged.frame = "."
if feature.featuretype != current_merged.featuretype:
current_merged.featuretype = "sequence_feature"
if feature.start < current_merged.start:
# Extends prior, so set a new start position
current_merged.start = feature.start
if feature.end > current_merged.end:
# Extends further, so set a new stop position
current_merged.end = feature.end
else:
yield _finalize_merge(current_merged, feature_children)
current_merged = feature
feature_children = []
last_id = None
if current_merged:
yield _finalize_merge(current_merged, feature_children)
def merge_all(
self,
merge_order=("seqid", "featuretype", "strand", "start"),
merge_criteria=(mc.seqid, mc.overlap_end_inclusive, mc.strand, mc.feature_type),
featuretypes_groups=(None,),
exclude_components=False,
):
"""
Merge all features in database according to criteria.
Merged features will be assigned as children of the merged record.
The resulting records are added to the database.
Parameters
----------
merge_order : list
Ordered list of columns with which to group features before evaluating criteria
merge_criteria : list
List of merge criteria callbacks. See merge().
featuretypes_groups : list
iterable of sets of featuretypes to merge together
exclude_components : bool
True: child features will be discarded. False to keep them.
Returns
-------
list of merge features
"""
if not len(featuretypes_groups):
# Can't be empty
featuretypes_groups = (None,)
result_features = []
# Merge features per featuregroup
for featuregroup in featuretypes_groups:
for merged in self.merge(
self.all_features(featuretype=featuregroup, order_by=merge_order),
merge_criteria=merge_criteria,
):
# If feature is result of merge
if merged.children:
self._insert(merged, self.conn.cursor())
if exclude_components:
# Remove child features from DB
self.delete(merged.children)
else:
# Add child relations to DB
for child in merged.children:
self.add_relation(merged, child, 1, child_func=assign_child)
result_features.append(merged)
else:
pass # Do nothing, feature is already in DB
return result_features
[docs]
def children_bp(
self,
feature,
child_featuretype="exon",
merge=False,
merge_criteria=(mc.seqid, mc.overlap_end_inclusive, mc.strand, mc.feature_type),
**kwargs
):
"""
Total bp of all children of a featuretype.
Useful for getting the exonic bp of an mRNA.
Parameters
----------
feature : str or Feature instance
child_featuretype : str
Which featuretype to consider. For example, to get exonic bp of an
mRNA, use `child_featuretype='exon'`.
merge : bool
Whether or not to merge child features together before summing
them.
merge_criteria : list
List of merge criteria callbacks. All must evaluate to True in
order for a feature to be merged. Only used if merge=True. When
modifying this argument, you may want to use:
from gffutils import merge_criteria as mc
to access the available callbacks.
Returns
-------
Integer representing the total number of bp.
"""
if kwargs:
if "ignore_strand" in kwargs:
raise ValueError(
"'ignore_strand' has been deprecated; please use "
"merge_criteria to control how features should be merged. "
"E.g., leave out the 'mc.strand' criteria to ignore strand."
)
else:
raise TypeError(
"merge() got unexpected keyword arguments '{}'".format(
kwargs.keys()
)
)
children = self.children(
feature, featuretype=child_featuretype, order_by="start"
)
if merge:
children = self.merge(children, merge_criteria=merge_criteria)
total = 0
for child in children:
total += len(child)
return total
[docs]
def bed12(
self,
feature,
block_featuretype=["exon"],
thick_featuretype=["CDS"],
thin_featuretype=None,
name_field="ID",
color=None,
):
"""
Converts `feature` into a BED12 format.
GFF and GTF files do not necessarily define genes consistently, so this
method provides flexiblity in specifying what to call a "transcript".
Parameters
----------
feature : str or Feature instance
In most cases, this feature should be a transcript rather than
a gene.
block_featuretype : str or list
Which featuretype to use as the exons. These are represented as
blocks in the BED12 format. Typically 'exon'.
Use the `thick_featuretype` and `thin_featuretype` arguments to
control the display of CDS as thicker blocks and UTRs as thinner
blocks.
Note that the features for `thick` or `thin` are *not*
automatically included in the blocks; if you do want them included,
then those featuretypes should be added to this `block_features`
list.
If no child features of type `block_featuretype` are found, then
the full `feature` is returned in BED12 format as if it had
a single exon.
thick_featuretype : str or list
Child featuretype(s) to use in order to determine the boundaries of
the "thick" blocks. In BED12 format, these represent coding
sequences; typically this would be set to "CDS". This argument is
mutually exclusive with `thin_featuretype`.
Specifically, the BED12 thickStart will be the start coord of the
first `thick` item and the thickEnd will be the stop coord of the
last `thick` item.
thin_featuretype : str or list
Child featuretype(s) to use in order to determine the boundaries of
the "thin" blocks. In BED12 format, these represent untranslated
regions. Typically "utr" or ['three_prime_UTR', 'five_prime_UTR'].
Mutually exclusive with `thick_featuretype`.
Specifically, the BED12 thickStart field will be the stop coord of
the first `thin` item and the thickEnd field will be the start
coord of the last `thin` item.
name_field : str
Which attribute of `feature` to use as the feature's name. If this
field is not present, a "." placeholder will be used instead.
color : None or str
If None, then use black (0,0,0) as the RGB color; otherwise this
should be a comma-separated string of R,G,B values each of which
are integers in the range 0-255.
"""
if thick_featuretype and thin_featuretype:
raise ValueError(
"Can only specify one of `thick_featuertype` or " "`thin_featuretype`"
)
exons = list(
self.children(feature, featuretype=block_featuretype, order_by="start")
)
if len(exons) == 0:
exons = [feature]
feature = self[feature]
first = exons[0].start
last = exons[-1].stop
if first != feature.start:
raise ValueError(
"Start of first exon (%s) does not match start of feature (%s)"
% (first, feature.start)
)
if last != feature.stop:
raise ValueError(
"End of last exon (%s) does not match end of feature (%s)"
% (last, feature.stop)
)
if color is None:
color = "0,0,0"
color = color.replace(" ", "").strip()
# Use field names as defined at
# http://genome.ucsc.edu/FAQ/FAQformat.html#format1
chrom = feature.chrom
chromStart = feature.start - 1
chromEnd = feature.stop
orig = constants.always_return_list
constants.always_return_list = True
try:
name = feature[name_field][0]
except KeyError:
name = "."
constants.always_return_list = orig
score = feature.score
if score == ".":
score = "0"
strand = feature.strand
itemRgb = color
blockCount = len(exons)
blockSizes = [len(i) for i in exons]
blockStarts = [i.start - 1 - chromStart for i in exons]
if thick_featuretype:
thick = list(
self.children(feature, featuretype=thick_featuretype, order_by="start")
)
if len(thick) == 0:
thickStart = feature.start
thickEnd = feature.stop
else:
thickStart = thick[0].start - 1 # BED 0-based coords
thickEnd = thick[-1].stop
if thin_featuretype:
thin = list(
self.children(feature, featuretype=thin_featuretype, order_by="start")
)
if len(thin) == 0:
thickStart = feature.start
thickEnd = feature.stop
else:
thickStart = thin[0].stop
thickEnd = thin[-1].start - 1 # BED 0-based coords
tst = chromStart + blockStarts[-1] + blockSizes[-1]
assert tst == chromEnd, "tst=%s; chromEnd=%s" % (tst, chromEnd)
fields = [
chrom,
chromStart,
chromEnd,
name,
score,
strand,
thickStart,
thickEnd,
itemRgb,
blockCount,
",".join(map(str, blockSizes)),
",".join(map(str, blockStarts)),
]
return "\t".join(map(str, fields))
def seqids(self):
"""
Yield the unique sequence IDs (chromosomes, contigs) observed in the
database.
"""
c = self.conn.cursor()
c.execute(
"""
SELECT DISTINCT seqid from features
"""
)
for (i,) in c:
yield i
# Recycle the docs for _relation so they stay consistent between parents()
# and children()
children.__doc__ = children.__doc__.format(_relation_docstring=_relation.__doc__)
parents.__doc__ = parents.__doc__.format(_relation_docstring=_relation.__doc__)
# Add the docs for methods that call helpers.make_query()
for method in [parents, children, features_of_type, all_features]:
method.__doc__ = method.__doc__.format(_method_doc=_method_doc)