Is it possible to create a db only for one chromosome?

[lldelisle]

Hi,
I am using gffutils to plot genes in a given region. This is very useful but quite slow because I am doing gffutils.create_db(file_path, ':memory:') which creates a database for the full genome. Is it possible to create a database only for a region?
Thank you very much,

Lucille

[daler]

Sorry for the late reply . . .

GFF and GTF files don't have a defined order, so in the general case you'd need to read the whole file even if you only want a subset. Even if you know your files are sorted, gffutils doesn't have a good mechanism of extracting the regions you want without first entering them into the database.

If your plotting requires a gffutils database but only of a small genomic region, I would use bedtools or pybedtools to do the subsetting first, and then only make a database on the subset.

Here's some code to get you started:

import pybedtools
import gffutils

gff = pybedtools.BedTool('full_genome_annotation.gff')

def db_for_region(region):

    # only the GFF features that intersect regions.
    # May want u=True to keep entire GFF features,
    # or keep as-is to truncate them to the region
    subset = gff.intersect(region)

    # create database on subset
    return gffutils.create_db(subset.fn, ':memory:')

region = [pybedtools.Interval('chr1', 1, 100)]
# or
region = pybedtools.BedTool('chr1 1 100', from_string=True)

db = db_for_region(region)

[lldelisle]

Many thanks