genomelake

Efficient random access to genomic data for deep learning models.

Supports the following types of input data:

bigwig
DNA sequence

genomelake extracts signal from genomic inputs in provided BED intervals.

Requirements

python 2.7, 3.5, or 3.6
tiledb
bcolz
cython
numpy
pybedtools
pysam

Installation

Clone the repository and run:

python setup.py install

Getting started: training a protein-DNA binding model

Extract genome-wide sequence data into a genomelake data source:

from genomelake.backend import extract_fasta_to_file

genome_fasta = "/mnt/data/annotations/by_release/hg19.GRCh37/hg19.genome.fa"
genome_data_directory = "./hg19_data_directory"
extract_fasta_to_file(genome_fasta, genome_data_directory)

Using a BED intervals file with labels, a genome data source, and genomelake's ArrayExtractor, generate input DNA sequences and labels:

import pybedtools
from genomelake.extractors import ArrayExtractor
import numpy as np

def batch_iter(iterable, batch_size):
    it = iter(iterable)
    try:
        while True:
            values = []
            for n in range(batch_size):
                values += (next(it),)
            yield values
    except StopIteration:
        yield values

def generate_inputs_and_labels(intervals_file, data_source, batch_size=128):
    bt = pybedtools.BedTool(intervals_file)
    extractor = ArrayExtractor(data_source)
    intervals_generator = batch_iter(bt, batch_size)
    for intervals_batch in intervals_generator:
        inputs = extractor(intervals_batch)
        labels = []
        for interval in intervals_batch:
            labels.append(float(interval.name))
        labels = np.array(labels)
        yield inputs, labels

Train a keras model of JUND binding to DNA using 101 base pair intervals and labels in ./examples/JUND.HepG2.chr22.101bp_intervals.tsv.gz:

from keras.models import Sequential
from keras.layers import Conv1D, Flatten, Dense

intervals_file = "./examples/JUND.HepG2.chr22.101bp_intervals.tsv.gz"
inputs_labels_generator = generate_inputs_and_labels(intervals_file, genome_data_directory)

model = Sequential()
model.add(Conv1D(15, 25, input_shape=(101, 4)))
model.add(Flatten())
model.add(Dense(1, activation='sigmoid'))

model.compile(loss='binary_crossentropy', optimizer='adam', metrics=['accuracy'])
model.fit_generator(inputs_labels_generator, steps_per_epoch=100)

Here is the expected result:

100/100 [==============================] - 7s - loss: 0.0584 - acc: 0.9905

License

genomelake is released under the BSD-3 license. See LICENSE for details.