=========

Tutorial

=========

This tutorial is split in three parts.

Part I treats the Genomic Dataset

that are available through Janggu which can be directly consumed

by your keras model.

The tutorial illustates how to access genomics

data from various widely used file formats, including FASTA, BAM, BIGWIG, BED

and GFF for the purpose of using them as input to a deep learning application.

It illustrates a range of parameters to adapt the read out of genomics data

and it shows how the predictions or feature activities of a neural network

in numpy format can be converted to a genomic coverage representation

that can in turn be exported to BIGWIG file format

or visualized directly via a genome browser-like plot.

Part II treats utilities for defining a neural networks based on keras.

Part III illustrates Janggu's evaluation utilities.

Complementary to this tutorial, the janggu repository contains

a number of jupyter notebooks that illustrate for example with keras or

sklearn:

.. _notebook_tutorials:

+----------------------------------------------------+

| Example notebooks                                  |

+====================================================+

| `keras cnn example`_                               |

+----------------------------------------------------+

| `keras cnn example with JangguSequence`_           |

+----------------------------------------------------+

| `sklearn example`_                                 |

+----------------------------------------------------+

| `pytorch example`_                                 |

+----------------------------------------------------+

| `janggu example I`_                                |

+----------------------------------------------------+

| `janggu example II`_                               |

+----------------------------------------------------+

| `reusing datasets with view`_                      |

+----------------------------------------------------+

| `hyperparameter optimization`_                     |

+----------------------------------------------------+

| `randomizing HDF5 data`_                           |

+----------------------------------------------------+

| `variant effect prediction`_                       |

+----------------------------------------------------+

| `variant effect prediction - part II`_             |

+----------------------------------------------------+

| `plotting genome coverage`_                        |

+----------------------------------------------------+

.. _`keras cnn example`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/keras_convnet_example.ipynb

.. _`keras cnn example with JangguSequence`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/keras_convnet_example_w_sequence.ipynb

.. _`sklearn example`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/sklearn_example_with_kmers.ipynb

.. _`pytorch example`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/pytorch_convnet_example.ipynb

.. _`janggu example I`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/classify_nucleotide_sequences.ipynb

.. _`janggu example II`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/janggu_convnet_examples.ipynb

.. _`reusing datasets with view`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/janggu_convnet_examples_with_hdf5.ipynb

.. _`randomizing HDF5 data`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/janggu_convnet_examples_with_hdf5.ipynb

.. _`variant effect prediction`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/variant_effect_prediction.ipynb

.. _`variant effect prediction - part II`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/variant_effect_prediction-part2.ipynb

.. _`plotting genome coverage`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/plot_coverage.ipynb

.. _`hyperparameter optimization`: https://nbviewer.jupyter.org/github/BIMSBbioinfo/janggu/blob/master/src/examples/hyperparameter_optimization_w_janggu.ipynb

Furthermore, use cases

for predicting JunD binding, training and adapting published genomics models

as well as a regression model example are demonstrated in

the supplementary repository: `Janggu use cases`_

.. _`Janggu use cases`: https://github.com/wkopp/janggu_usecases

Part I) Introduction to Genomic Datasets

-----------------------------------------

.. sidebar:: Genomic Datasets

   Most of the parameters are consistent across

   :class:`Bioseq` and :class:`Cover`.

:mod:`janggu.data` provides Dataset classes

that can be used for

training and evaluating neural networks.

Of particular importance are the Genomics-specific dataset,

:class:`Bioseq` and :class:`Cover` which

allow easy access to genomics data,

including DNA sequences or coverage information.

Apart from accessing raw genomics data, :code:`Janggu`

also facilitates a method for converting an ordinary

numpy array (e.g. predictions obtained from a neural net)

to a :code:`Cover` object. This enables the user to export the predictions

as BIGWIG format or interactively plot genome browser tracks.

In this tutorial, we demonstrate some of the key functionality of

Janggu. Further details are available in :ref:`storage`

and :ref:`reference-label`.

Bioseq

^^^^^^^^^^

The :class:`Bioseq` can be used to load nucleotide

or protein sequence data from

fasta files or from a reference genome

along with a set of genomic coordinates defining the region of interest (ROI).

The class facilitates access the

*one-hot encoding* representation of the sequences.

Specifically,

the *one-hot encoding* is represented as a

4D array with dimensions corresponding

to :code:`(region, region_length, 1, alphabet_size)`.

The Bioseq offers a number of features:

1. Strand-specific sequence extraction (if DNA sequences are extracted from the reference genome)

2. Higher-order one-hot encoding, e.g. di-nucleotide based

Sequences can be loaded in two ways: using

:code:`Bioseq.create_from_seq` or

:code:`Bioseq.create_from_refgenome`.

The former constructor method can be used to load

DNA or protein sequences from fasta files directly

or from as list of :code:`Bio.SeqRecord.SeqRecord` entries.

An example is shown below:

.. code-block:: python

   from pkg_resources import resource_filename

   from janggu.data import Bioseq

   fasta_file = resource_filename('janggu',

                                  'resources/sample.fa')

   dna = Bioseq.create_from_seq(name='dna',

                                fastafile=fasta_file)

   # there are 3897 sequences in the in sample.fa

   len(dna)

   # Each sequence is 200 bp of length

   dna.shape  # is (3897, 200, 1, 4)

   # One-hot encoding for the first 10 bases of the first region

   dna[0][0, :10, 0, :]

   #array([[0, 1, 0, 0],

   #       [1, 0, 0, 0],

   #       [0, 1, 0, 0],

   #       [1, 0, 0, 0],

   #       [0, 0, 1, 0],

   #       [0, 1, 0, 0],

   #       [1, 0, 0, 0],

   #       [0, 0, 1, 0],

   #       [1, 0, 0, 0],

   #       [0, 0, 1, 0]], dtype=int8)

Furthermore, it is possible to trim variable sequence length using

the :code:`fixedlen` option. If specfied, all sequences will be truncated

or zero-padded to length `fixedlen`. For example,

.. code-block:: python

   dna = Bioseq.create_from_seq(name='dna',

                                fastafile=fasta_file,

                                fixedlen=205)

   # Each sequence is 205 bp of length

   dna.shape  # is (3897, 205, 1, 4)

   # the last 5 position were zero padded

   dna[0][0, -6:, 0, :]

   #array([[1, 0, 0, 0],

   #       [0, 0, 0, 0],

   #       [0, 0, 0, 0],

   #       [0, 0, 0, 0],

   #       [0, 0, 0, 0],

   #       [0, 0, 0, 0]], dtype=int8)

Alternatively, nucleotide sequences can be

obtained from a reference genome directly along with

a BED or GFF file that indicates the region of interest (ROI).

If each interval in the BED-file already corresponds

to a 'datapoint' that shall be consumed during training, like it

is the case for 'sample_equalsize.bed', the associated DNA sequences

can be loaded according to

.. code-block:: python

   roi = resource_filename('janggu',

                           'resources/sample_equalsize.bed')

   refgenome = resource_filename('janggu',

                                 'resources/sample_genome.fa')

   dna = Bioseq.create_from_refgenome(name='dna',

                                      refgenome=refgenome,

                                      roi=roi)

   dna.shape  # is (4, 200, 1, 4)

   # One-hot encoding of the first 10 nucleotides in region 0

   dna[0][0, :10, 0, :]

   #array([[0, 1, 0, 0],

   #       [1, 0, 0, 0],

   #       [0, 1, 0, 0],

   #       [1, 0, 0, 0],

   #       [0, 0, 1, 0],

   #       [0, 1, 0, 0],

   #       [1, 0, 0, 0],

   #       [0, 0, 1, 0],

   #       [1, 0, 0, 0],

   #       [0, 0, 1, 0]], dtype=int8)

Sometimes it is more convenient to provide the ROI

as a set of variable-sized broad intervals

(e.g. chr1:10000-50000 and chr3:4000-8000)

which should be divided into sub-intervals

of equal length (e.g. of length 200 bp).

This can be achieved

by explicitly specifying a desired :code:`binsize`

and :code:`stepsize` as shown below:

.. code-block:: python

   roi = resource_filename('janggu',

                           'resources/sample.bed')

   # loading non-overlapping intervals

   dna = Bioseq.create_from_refgenome(name='dna',

                                      refgenome=refgenome,

                                      roi=roi,

                                      binsize=200,

                                      stepsize=200)

   dna.shape  # is (100, 200, 1, 4)

   # loading mutually overlapping intervals

   dna = Bioseq.create_from_refgenome(name='dna',

                                      refgenome=refgenome,

                                      roi=roi,

                                      binsize=200,

                                      stepsize=50)

   dna.shape  # is (394, 200, 1, 4)

The argument :code:`flank` can be used to extend

the intervals up and downstream by a given length

.. code-block:: python

   dna = Bioseq.create_from_refgenome(name='dna',

                                      refgenome=refgenome,

                                      roi=roi,

                                      binsize=200,

                                      stepsize=200,

                                      flank=100)

   dna.shape  # is (100, 400, 1, 4)

Finally, sequences can be represented using **higher-order**

one-hot representation using the :code:`order` argument. An example

of a di-nucleotide-based one-hot representation is shown below

.. code-block:: python

   dna = Bioseq.create_from_refgenome(name='dna',

                                      refgenome=refgenome,

                                      roi=roi,

                                      binsize=200,

                                      stepsize=200,

                                      order=2)

   # is (100, 199, 1, 16)

   # that is the last dimension represents di-nucleotides

   dna.shape

   dna.conditions

   # ['AA', 'AC', 'AG', 'AT', 'CA', 'CC', 'CG', 'CT', 'GA', 'GC', 'GG', 'GT', 'TA', 'TC', 'TG', 'TT']

   dna[0][0, :5, 0, :]

   #array([[0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],

   #       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],

   #       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],

   #       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],

   #       [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0]], dtype=int8)

Cover

^^^^^

:class:`Cover` can be utilized to fetch different kinds of

coverage data from commonly used data formats,

including BAM, BIGWIG, BED and GFF.

Coverage data is stored as a 4D array with dimensions corresponding

to :code:`(region, region_length, strand, condition)`.

The following examples illustrate some use cases for :class:`Cover`,

including loading, normalizing coverage data.

Additional features are described in the :ref:`reference-label`.

**Loading read count coverage from BAM files** is supported for

single-end and paired-end alignments. For the single-end case

reads are counted on the 5'-end and and for paired-end alignments,

reads are optionally counted at the mid-points or 5' ends of the first mate.

The following example illustrate how to extract base-pair resolution coverage

with and without strandedness.

.. code:: python

   from janggu.data import Cover

   bam_file = resource_filename('janggu',

                                'resources/sample.bam')

   roi = resource_filename('janggu',

                           'resources/sample.bed')

   cover = Cover.create_from_bam('read_count_coverage',

                                 bamfiles=bam_file,

                                 binsize=200,

                                 stepsize=200,

                                 roi=roi)

   cover.shape  # is (100, 200, 2, 1)

   cover[0]  # coverage of the first region

   # Coverage regardless of read strandedness

   # sums reads from both strand.

   cover = Cover.create_from_bam('read_coverage',

                                 bamfiles=bam_file,

                                 binsize=200,

                                 stepsize=200,

                                 stranded=False,

                                 roi=roi)

   cover.shape  # is (100, 200, 1, 1)

Sometimes it is desirable to determine the read

count coverage in say 50 bp bins which can be

controlled by the :code:`resolution` argument.

Consequently, note that the second dimension amounts

to length 4 using `binsize=200` and `resolution=50` in the following example

.. code:: python

   # example with resolution=200 bp

   cover = Cover.create_from_bam('read_coverage',

                                 bamfiles=bam_file,

                                 binsize=200,

                                 resolution=50,

                                 roi=roi)

   cover.shape  # is (100, 4, 2, 1)

It might be desired to aggregate reads across entire interval

rather than binning the genome to equally sized bins of

length :code:`resolution`. An example application for this would

be to count reads per possibly variable-size regions (e.g. genes).

This can be achived by :code:`resolution=None` which results

in the second dimension being collapsed to a length of one.

.. code:: python

   # example with resolution=None

   cover = Cover.create_from_bam('read_coverage',

                                 bamfiles=bam_file,

                                 binsize=200,

                                 resolution=None,

                                 roi=roi)

   cover.shape  # is (100, 1, 2, 1)

Similarly, if strandedness is not relevant we may use

.. code:: python

   # example with resolution=None without strandedness

   cover = Cover.create_from_bam('read_coverage',

                                 bamfiles=bam_file,

                                 binsize=200,

                                 resolution=None,

                                 stranded=False,

                                 roi=roi)

   cover.shape  # is (100, 1, 1, 1)

Finally, it is possible to normalize the coverage profile, e.g.

to account for differences in sequencing depth across experiments

using the :code:`normalizer` argument

.. code:: python

   # example with resolution=None without strandedness

   cover = Cover.create_from_bam('read_coverage',

                                 bamfiles=bam_file,

                                 binsize=200,

                                 resolution=None,

                                 stranded=False,

                                 normalizer='tpm',

                                 roi=roi)

   cover.shape  # is (100, 1, 1, 1)

More details on alternative normalization

options are discussed in :ref:`storage`.

**Loading signal coverage from BIGWIG files**

can be achieved analogously:

.. code-block:: python

   roi = resource_filename('janggu',

                           'resources/sample.bed')

   bw_file = resource_filename('janggu',

                               'resources/sample.bw')

   cover = Cover.create_from_bigwig('bigwig_coverage',

                                    bigwigfiles=bw_file,

                                    roi=roi,

                                    binsize=200,

                                    stepsize=200)

   cover.shape  # is (100, 200, 1, 1)

When applying signal aggregation using e.g :code:`resolution=50` or :code:`resolution=None`,

additionally, the aggregation method can be specified using

the :code:`collapser` argument.

For example, in order to represent the resolution sized

bin by its mean signal the following snippet may be used:

.. code-block:: python

   cover = Cover.create_from_bigwig('bigwig_coverage',

                                    bigwigfiles=bw_file,

                                    roi=roi,

                                    binsize=200,

                                    resolution=None,

                                    collapser='mean')

   cover.shape  # is (100, 1, 1, 1)

More details on alternative collapse

options are discussed in :ref:`storage`.

**Coverage from a BED files** is largely analogous to extracting coverage

information from BAM or BIGWIG files, but in addition the :code:`mode` option

enables to interpret BED-like files in various ways:

1. :code:`mode='binary'` Presence/Absence mode interprets the ROI as the union of positive and negative cases in a binary classification setting and regions contained in :code:`bedfiles` as positive examples.

2. :code:`mode='score'` reads out the real-valued score field value from the associated regions.

3. :code:`mode='score_category'` transforms integer-valued scores into a categorical one-hot representation.

4. :code:`mode='name_category'` transforms the name field into a categorical one-hot representation.

5. :code:`mode='bedgraph'` reads in the score from a file in bedgraph format.

Examples of loading data from a BED file are shown below

.. code-block:: python

   roi = resource_filename('janggu',

                           'resources/sample.bed')

   score_file = resource_filename('janggu',

                                  'resources/scored_sample.bed')

   # binary mode (default)

   cover = Cover.create_from_bed('binary_coverage',

                                 bedfiles=score_file,

                                 roi=roi,

                                 binsize=200,

                                 stepsize=200,

                                 collapser='max',

                                 resolution=None)

   cover.shape

   # (100, 1, 1, 1)

   cover[4]

   # array([[[[1.]]]])

   # score mode

   cover = Cover.create_from_bed('score_coverage',

                                 bedfiles=score_file,

                                 roi=roi,

                                 binsize=200,

                                 stepsize=200,

                                 resolution=None,

                                 collapser='max',

                                 mode='score')

   cover.shape

   # (100, 1, 1, 1)

   cover[4]

   # array([[[[5.]]]])

   # scoreclass (or categorical) mode

   # Interprets the integer-valued score as class-label,

   # which will then be one-hot encoded.

   cover = Cover.create_from_bed('cat_coverage',

                                 bedfiles=score_file,

                                 roi=roi,

                                 binsize=200,

                                 stepsize=200,

                                 resolution=None,

                                 collapser='max',

                                 mode='score_category')

   # there are 4 categories

   cover.shape

   # (100, 1, 1, 4)

   # category names

   cover.conditions

   # ['1', '2', '4', '5']

   cover[4]

   # array([[[[0., 0., 0., 1.]]]])

   # nameclass mode

   # Interprets the name field as class-label.

   # The class labels will be one hot encoded.

   cover = Cover.create_from_bed('cat_coverage',

                                 bedfiles=score_file,

                                 roi=roi,

                                 binsize=200,

                                 stepsize=200,

                                 resolution=None,

                                 collapser='max',

                                 mode='name_category')

   cover.shape

   # (100, 1, 1, 2)

   # category names

   cover.conditions

   # ['state1', 'state2']

   cover[4]

   # array([[[[0., 1.]]]])

   # bedgraph-format mode

   bedgraph_file = resource_filename('janggu',

                                     'resources/sample.bedgraph')

   cover = Cover.create_from_bed('bedgraph_coverage',

                                 bedfiles=bedgraph_file,

                                 roi=roi,

                                 binsize=200,

                                 stepsize=200,

                                 resolution=None,

                                 collapser='max',

                                 mode='bedgraph')

   cover.shape

   # (100, 1, 1, 1)

   cover[4]

   # array([[[[0.5]]]])

Dataset wrappers

^^^^^^^^^^^^^^^^^

In addition to the core datset :code:`Bioseq` and :code:`Cover`, Janggu offers convenience wrappers

to transform them in various ways.

For instance, :code:`ReduceDim` can be used to convert a 4D coverage dataset into 2D table like object.

That is it may be used to transform the dimensions

:code:`(region, region_length, strand, condition)` to :code:`(region, condition)` by

aggregating over the middle two dimensions.

.. code:: python

   from janggu.data import ReduceDim

   cover.shape

   # (100, 1, 1, 1)

   data = ReduceDim(cover, aggregator='sum')

   data.shape

   # (100, 1)

Other dataset wrappers can be used in order to perform data augmentation, including

:code:`RandomSignalScale` and :code:`RandomOrientation` which can be used

to randomly alter the signal intensity during model fitting and randomly flipping

the 5' to 3' orientations of the coverage signal.

For more specialized cases, these wrappers might also be a good starting point

to derive or adapt from.

Using the Genomic Datasets with keras or sklearn

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The above mentioned datasets :code:`Bioseq` and :code:`Cover`

are directly compatible with keras and sklearn models. An illustration of a

simple convolutional neural network with keras is shown in

`keras cnn example`_.

Moreover, an example of a logistic regression model from sklearn used with Janggu

is shown in

`sklearn example`_.

Converting a Numpy array to :code:`Cover`

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

After having trained and performed predictions with a model, the data

is represented as numpy array. A convenient way to reassociate the

predictions with the genomic coordinates they correspond to is achieved

using :code:`create_from_array`.

.. code:: python

   import numpy as np

   # True labels may be obtained from a BED file

   cover = Cover.create_from_bigwig('cov',

                                     bigwigfiles=bw_file,

                                     roi=roi,

                                     binsize=200,

                                     resolution=50)

   # Let's pretend to have derived predictions from a NN

   # of the same shape

   predictions = np.random.randn(*cover.shape)*.1 + cover[:]

   # We can reassociate the predictions with the genomic coordinates

   # of a :code:`GenomicIndexer` (in this case, cover.gindexer).

   predictions = Cover.create_from_array('predictions',

                                         predictions, cover.gindexer)

Exporting and visualizing :code:`Cover`

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

After having converted the predictions or feature activities of a neural

network to a :code:`Cover` object, it is possible to export the results

as BIGWIG format for further investigation in a genome browser of your choice

.. code:: python

   # writes the predictions to a specified folder

   predictions.export_to_bigwig(output_dir = './')

which should result in a file 'predictions.Cond_0.bigwig'.

Furthermore, it is possible to visualize the tracks interactively

.. code:: python

   from janggu.data import LineTrack

   from janggu.data import plotGenomeTrack

   fig = plotGenomeTrack([LineTrack(cover), LineTrack(predictions)], 'chr1', 16000, 18000).figsave('coverage.png')

.. image:: coverage.png

   :width: 70%

   :alt: Coverage tracks

   :align: center

Part II) Building a neural network with Janggu

-----------------------------------------------

While the Genomic Dataset may be used directly with keras,

this part of the tutorial discusses the :class:`Janggu` wrapper class

for a keras model.

It offers the following features:

1. Building models using automatic input and output layer shape inference

2. Built-in logging functionality

3. Automatic evaluation through the attachment of Scorer callbacks

A list of examples can be found in the :ref:`Table <notebook_tutorials>` at the beginning.

.. sidebar:: Datasets are named

   Dataset names must match with the Input and Output layers of the neural

   network.

A neural network can be created by

instantiating a :class:`Janggu` object.

There are two ways of achieving this:

1. Similar as with `keras.models.Model`, a :class:`Janggu` object can be created from a set of native keras Input and Output layers, respectively.

2. Janggu offers a `Janggu.create` constructor method which helps to reduce redundant code when defining many rather similar models.

Example 1: Instantiate Janggu similar to keras.models.Model

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

.. sidebar:: **Model name**

   Model results,

   e.g. trained parameters, are automatically stored with the associated model name. To simplify the determination of a unique name for the model, Janggu automatically derives the model name based on a md5-hash of the network configuration. However, you can also specify a name yourself.

.. code-block:: python

  from keras.layers import Input

  from keras.layers import Dense

  from janggu import Janggu

  # Define neural network layers using keras

  in_ = Input(shape=(10,), name='ip')

  layer = Dense(3)(in_)

  output = Dense(1, activation='sigmoid',

                 name='out')(layer)

  # Instantiate model name.

  model = Janggu(inputs=in_, outputs=output)

  model.summary()

Example 2: Specify a model using a model template function

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

As an alternative to the above stated variant, it is also possible to specify

a network via a python function as in the following example

.. code-block:: python

   def model_template(inputs, inp, oup, params):

       inputs = Input(shape=(10,), name='ip')

       layer = Dense(params)(inputs)

       output = Dense(1, activation='sigmoid',

                      name='out')(layer)

       return inputs, output

   # Defines the same model by invoking the definition function

   # and the create constructor.

   model = Janggu.create(template=model_template,

                         modelparams=3)

The model template function must adhere to the

signature :code:`template(inputs, inp, oup, params)`.

Notice, that :code:`modelparams=3` gets passed on to :code:`params`

upon model creation. This allows to parametrize the network

and reduces code redundancy.

From the model template it is also possible to obtain

a keras model directly, rather than the Janggu model wrapper if this is preferred

.. code-block:: python

   from janggu import create_model

   # This will construct a keras model directly

   model = create_model(template=model_template,

                        modelparams=3)

Example 3: Automatic Input and Output layer extension

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

A second benefit to invoke :code:`Janggu.create` is that it can automatically

determine and append appropriate Input and Output layers to the network.

This means, only the network body remains to be defined.

.. code-block:: python

    import numpy as np

    from janggu import inputlayer, outputdense

    from janggu.data import Array

    # Some random data

    DATA = Array('ip', np.random.random((1000, 10)))

    LABELS = Array('out', np.random.randint(2, size=(1000, 1)))

    # inputlayer and outputdense automatically

    # extract dataset shapes and extend the

    # Input and Output layers appropriately.

    # That is, only the model body needs to be specified.

    @inputlayer

    @outputdense('sigmoid')

    def model_body_template(inputs, inp, oup, params):

        with inputs.use('ip') as layer:

            # the with block allows

            # for easy access of a specific named input.

            output = Dense(params)(layer)

        return inputs, output

    # create the model.

    model = Janggu.create(template=model_body_template,

                          modelparams=3,

                          inputs=DATA, outputs=LABELS)

    model.summary()

As is illustrated by the example, automatic Input and Output layer determination

can be achieved by using the decorators :code:`inputlayer` and/or

:code:`outputdense` which extract the layer dimensions from the

provided input and output Datasets in the create constructor.

Fit a neural network on DNA sequences

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In the previous sections, we learned how to acquire data and

how to instantiate neural networks. Now let's

create and fit a simple convolutional neural network that learns

to discriminate between two classes of sequences. In the following example

the sample sequences are of length 200 bp each. `sample.fa` contains Oct4 CHip-seq

peaks and sample2.fa contains Mafk CHip-seq peaks. We shall use a simple

convolutional neural network with 30 filters of length 21 bp to learn

the sequence features that discriminate the two sets of sequences.

The example makes use of two more janggu utilities: First,

:code:`DnaConv2D` constitutes a keras layer wrapper that facilitates scanning

of both DNA strands with the same kernels. That is it simulataneously applies

a convolution and a cross-correlation and aggregates the resulting activities.

Second, the example illustrates the dataset wrapper :code:`ReduceDim` which

allows to collapse 4D the signal contained in the Cover object

across the sequence length and strand dimension. The result is yields a 2D

table-like dataset which is used in the subsequent model fitting example.

.. code:: python

   from keras.layers import Conv2D

   from keras.layers import GlobalAveragePooling2D

   from janggu import inputlayer

   from janggu import outputconv

   from janggu import DnaConv2D

   from janggu.data import ReduceDim

   # load the dataset which consists of

   # 1) a reference genome

   REFGENOME = resource_filename('janggu', 'resources/pseudo_genome.fa')

   # 2) ROI contains regions spanning positive and negative examples

   ROI_FILE = resource_filename('janggu', 'resources/roi_train.bed')

   # 3) PEAK_FILE only contains positive examples

   PEAK_FILE = resource_filename('janggu', 'resources/scores.bed')

   # DNA sequences are loaded directly from the reference genome

   DNA = Bioseq.create_from_refgenome('dna', refgenome=REFGENOME,

                                      roi=ROI_FILE,

                                      binsize=200)

   # Classification labels over the same regions are loaded

   # into the Coverage dataset.

   # It is important that both DNA and LABELS load with the same

   # binsize, stepsize to ensure

   # the correct correspondence between both datasets.

   # Finally, the ReduceDim dataset wrapper transforms the 4D Coverage

   # object into a 2D table like object (regions by conditions)

   LABELS = ReduceDim(Cover.create_from_bed('peaks', roi=ROI_FILE,

                                  bedfiles=PEAK_FILE,

                                  binsize=200,

                                  resolution=None), aggregator='mean')

   # 2. define a simple conv net with 30 filters of length 15 bp

   # and relu activation.

   # outputconv as opposed to outputdense will put a conv layer as output

   @inputlayer

   @outputdense('sigmoid')

   def double_stranded_model(inputs, inp, oup, params):

       with inputs.use('dna') as layer:

           # The DnaConv2D wrapper can be used with Conv2D

           # to scan both DNA strands with the weight matrices.

           layer = DnaConv2D(Conv2D(params[0], (params[1], 1),

                                    activation=params[2]))(layer)

       output = GlobalAveragePooling2D(name='motif')(layer)

       return inputs, output

   # 3. instantiate and compile the model

   model = Janggu.create(template=double_stranded_model,

                         modelparams=(30, 15, 'relu'),

                         inputs=DNA, outputs=LABELS)

   model.compile(optimizer='adadelta', loss='binary_crossentropy',

                 metrics=['acc'])

   # 4. fit the model

   model.fit(DNA,LABELS,epochs=100)

An illustration of the network architecture is depicted below.

Upon creation of the model a network depiction is

automatically produced in :code:`<results_root>/models` which is illustrated

below

.. image:: dna_peak.png

   :width: 70%

   :alt: Prediction from DNA to peaks

   :align: center

After the model has been trained, the model parameters and the

illustration of the architecture are stored in :code:`<results_root>/models`.

Furthermore, information about the model fitting, model and dataset dimensions

are written to :code:`<results_root>/logs`.

Note that in the example above the output dimensionality of the network is 4D.

However, it might be more convenient at times to remove the single dimensional

elements of the array.

This can be achieved by wrapping the LABELS dataset using :code:`ReduceDim`.

In this case the example becomes

.. code:: python

   @inputlayer

   @outputdense('sigmoid')

   def double_stranded_model(inputs, inp, oup, params):

       with inputs.use('dna') as layer:

           # The DnaConv2D wrapper can be used with Conv2D

           # to scan both DNA strands with the weight matrices.

           layer = DnaConv2D(Conv2D(params[0], (params[1], 1),

                                    activation=params[2]))(layer)

       output = GlobalAveragePooling2D(name='motif')(layer)

       return inputs, output

   # 3. instantiate and compile the model

   model = Janggu.create(template=double_stranded_model,

                         modelparams=(30, 15, 'relu'),

                         inputs=DNA, outputs=LABELS)

   model.compile(optimizer='adadelta', loss='binary_crossentropy',

                 metrics=['acc'])

   # 4. fit the model

   model.fit(DNA, LABELS, epochs=100)

Part III) Evaluation and interpretation of the model

-----------------------------------------------------

Janggu supports various methods to evaluate and interprete a trained model,

including evaluating summary scores, inspecting the results in

the built-in genome browser (see Part I), evaluating the integrated gradients

which allows to visualized input feature importance and by

offering support for variant effect predictions.

In this last part we will illustrate these aspects.

Evaluation of summary scores

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

After the model has been trained, the quality of the predictions

is usually summarized by its agreement with the ground truth, e.g. by

evaluating the area under the ROC curve in a binary classification application

or by computing the correlation between predictions and targets in a regression

setting.

For some commonly used evaluation criteria, the evaluate method directly allows

to determine and export the given metric results.

For example, for a classification task the following line

evaluates the ROC and PRC and exports a figure and a tsv file, respectively,

for each measure.

.. code-block:: python

   model.evaluate(DNA_TEST, LABELS_TEST, callbacks=['roc', 'prc', 'auprc', 'auroc'])

The results are stored in :code:`<results_root>/evaluation/{roc,prc}.png`

as well as :code:`<results_root>/evaluation/{auroc,auprc}.tsv`.

Furthermore, for a regression setting it is possible to invoke

.. code-block:: python

   model.evaluate(DNA_TEST, LABELS_TEST, callbacks=['cor', 'mae', 'mse', 'var_explained'])

which evaluates the Pearson's correlation, the mean absolute error, the mean squared error

and the explained variance into tsv files.

It is also possible to customize the scoring callbacks by instantiating a

:code:`Scorer` objects which can be passed to

:code:`model.evaluate` and :code:`model.predict`. Further details about

customizing the scoring callbacks are given in :doc:`custom_scorer`.

Input feature importance

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

In order to inspect what the model has learned,

it is possible to identify the most important features in the input space

using the integrated gradients method.

This is illustrated on a toy example for discriminating Oct4 and Mafk binding sites (see

`variant effect prediction`_).

Variant effect prediction

^^^^^^^^^^^^^^^^^^^^^^^^^^^

In order to measure the effect of single nucleotide variant on the predict

network output can be tested via the :code:`Janggu.predict_variant_effect`

based on a Bioseq object and single nucleotide variants in VCF format.

This method evaluates the network for each variant (using its sequence context)

as well as its respective reference sequence.

As a result, an hdf5 file and a bed file will be produced which

contain the network predictions for each variant and the associated genomic

loci.

An illustration of the variant effect prediction in the notebook (see

`variant effect prediction`_).

Browse through the results

^^^^^^^^^^^^^^^^^^^^^^^^^^

Finally, after you have fitted and evaluated your results

you can browse through the results in an web browser of your choice.

To this end, first start the web application server

::

   janggu -path <results-root>

Then you can inspect the outputs in a browser of your choice

(default: localhost:8050)

.. image:: janggu_example.png

   :width: 70%

   :alt: Prediction from DNA to peaks

   :align: center