model_file = '../data/models/pretrained_model.th'
seqs_file = '../data/encode_roadmap.h5'

"""
Explanation: Once we've trained a model, we might want to better understand what sequence motifs the first convolutional layer has discovered and how it's using them. Basset offers two methods to help users explore these filters.
You'll want to double check that you have the Tomtom motif comparison tool for the MEME suite installed. Tomtom provides rigorous methods to compare the filters to a database of motifs. You can download it from here: http://meme-suite.org/doc/download.html
To run this tutorial, you'll need to either download the pre-trained model from https://www.dropbox.com/s/rguytuztemctkf8/pretrained_model.th.gz and preprocess the consortium data, or just substitute your own files here:
End of explanation
"""


import subprocess

cmd = 'basset_motifs.py -s 1000 -t -o motifs_out %s %s' % (model_file, seqs_file)
subprocess.call(cmd, shell=True)

"""
Explanation: First, we'll run basset_motifs.py, which will extract a bunch of basic information about the first layer filters. The script takes an HDF file (such as any preprocessed with preprocess_features.py) and samples sequences from the test set. It sends those sequences through the neural network and examines its hidden unit values to describe what they're doing.

-m specifies the Tomtom motif database: CIS-BP Homo sapiens database by default.
-s specifies the number of sequences to sample. 1000 is fast and sufficient.
-t asks the script to trim uninformative positions off the filter ends.
End of explanation
"""


# actual file is motifs_out/filter9_heat.pdf

from IPython.display import Image
Image(filename='motifs_eg/filter9_heat.png') 

"""
Explanation: Now there's plenty of information output in motifs_out. My favorite way to get started is to open the HTML file output by Tomtom's comparison of the motifs to a database. It displays all of the motifs and their database matches in a neat table.
Before we take a look though, let me describe where these position weight matrices came from. Inside the neural network, the filters are reprsented by real-valued matrices. Here's one:
End of explanation
"""


# actual file is motifs_out/filter9_logo.eps

Image(filename='motifs_eg/filter9_logo.png') 

"""
Explanation: Although it's matrix of values, this doesn't quite match up with the conventional notion of a position weight matrix that we typically use to represent sequences motifs in genome biology. To make that, basset_motifs.py pulls out the underlying sequences that activate the filters in the test sequences and passes that to weblogo.
End of explanation
"""


!sort -k6 -gr motifs_out/table.txt | head -n20

"""
Explanation: The Tomtom output will annotate the filters, but it won't tell you how the model is using them. So alongside it, let's print out a table of the filters with the greatest output standard deviation over this test set. Variance correlates strongly with how influential the filter is for making predictions.
The columns here are:
1. Index
2. Optimal sequence
3. Tomtom annotation
4. Information content
5. Activation mean
6. Activation standard deviaion
End of explanation
"""


!open motifs_out/tomtom/tomtom.html

"""
Explanation: As I discuss in the paper, unannotated low complexity filters tend to rise to the top here because low order sequence complexity influence accessibility.
The Tomtom output HTML is here:
End of explanation
"""


cmd = 'basset_motifs_infl.py'
cmd += ' -m motifs_out/table.txt'
cmd += ' -s 2000 -b 500'
cmd += ' -o infl_out'
cmd += ' --subset motifs_eg/primary_cells.txt'
cmd += ' -t motifs_eg/cell_activity.txt'
cmd += ' --width 7 --height 40 --font 0.5'
cmd += ' %s %s' % (model_file, seqs_file)

subprocess.call(cmd, shell=True)

"""
Explanation: The other primary tool that we have to understand the filters is to remove the filter from the model and assess the impact. Rather than truly remove it and re-train, we can just nullify it within the model by setting all output from the filter to its mean. This way the model around it isn't drastically affected, but there's no information flowing through.
This analysis requires considerably more compute time, so I separated it into a different script. To give it a sufficient number of sequences to obtain a good estimate influence, I typically run it overnight. If your computer is using too much memory, decrease the batch size. I'm going to run here with 1,000 sequences, but I used 20,000 for the paper.
To get really useful output, the script needs a few additional pieces of information:
* -m specifies the table created by basset_motifs.py above.
* -s samples 2000 sequences
* -b sets the batch_size to 500
* -o is the output directory
* --subset limits the cells displayed to those listed in the file.
* -t specifies a table where the second column is the target labels.
* --weight, --height, --font adjust the heatmap
End of explanation
"""


