Nimrod

This is a NN designed as a binary classifier of putative Transcription Factor binding sites (TFbs). It takes both ATACseq signal and genomic sequence as input.

To use it, you will need:

• A computer with a GPU that is compatible with tensorflow. The details of how to enable all the libraries that the GPU and tensorflow expect (CUDA drivers, CUDA toolkit) are system dependent so we cannot guide you more here. Nimrod should run without a GPU, but it'll be too slow for any meaningful analysis.
• Python3 and some python libraries (such as TensorFlow, pysam, pybedtools [see below])
• ATACseq and CHIPseq (for training) data in the same biological context (same dev. stage, or same tissue, or same cell line etc)
• Position Weight Matrix hits on your organism's genome.
• All work has been done in linux but the underlying tools (python, TensorFlow) should run in other OS as well, no guarantee

[[TOC]]

Installation

Before setting up the python installation, you will need to have the two following tools installed in your system. We don't use them directly but they are used by the python libraries that we use: * bedtools * samtools

Make sure you have the header files and static libraries for python-dev. This needs to be installed in the system and you'll need root access.

# in debia/ubuntu -like linux systems:
>> sudo apt-get install python3-dev

Your python3 installation will need the following libraries. We suggest creating a python virtual environment to keep things tidy.

# to make a new venv:
# (Be advised, you cannot move this folder after you've made it or things will break)
>> python3 -m 'venv' ./aNewVirtualEnvironment
# to activate it (you will need to do this everytime you want to work with this virtual environment)
>> source ./aNewVirtualEnvironment/bin/activate 
# now when you say 'python' in your command line, it refers to the python of this virtualenv. try this:
>> which python
# it should point to the venv folder
# to 'exit' the virtual env, just enter the following command
>> deactivate 
#but don't do this yet, we need to install the libraries

With your virtual environment activated, or on the global python (if you're sure you know what you're doing):

pip install --upgrade pip
# you can in theory run all these in one command but I find it goes smoother one by one
pip install numpy
pip install scipy
pip install pandas
pip install seaborn
pip install pysam
pip install pyBigWig
pip install pybedtools
pip install SharedArray
pip install tensorflow-gpu
pip install scikit-learn

After that, you can download this repository and install nimrod with pip

>> git clone https://gitlab.com/panosfirbas/nimrod.git
>> pip install ./nimrod

Nimrod should now be available to you in the command line, try:

>> nimrod --help

on Training

Good and bad TFbs

As a binary classifier, Nimrod will tell us yes or no for any putative TFbs that we present to it. The first step with a NN is training it by showing it what "good" and "bad" TFbs look like. This is a pretty vague task, as there are multiple ways to define "good" or "bad" sites, so this could change depending on what exactly you're trying to do.

In our work, we defined "good" sites as PWM hits that overlap a CHIP-seq peak and an ATACseq peak (as defined by a peak caller) in a specific biological context and "bad" as those sites that do not overlap a CHIP peak (we only test hits inside ATAC peaks, everything outside those was automatically dismissed). This means that the same PWM hit can be considered as good in one tissue but bad in another. This is possible thanks to the inclusion of ATACseq signal.

As an example, say we have data from a cell line. We scanned the genome for hits of a PWM (such as M1957 from here), with relaxed criteria. Out of those, we keep all the hits that overlap an ATACseq peak of our cell line and discard the rest. We check these PWM hits against CHIPseq peaks of our cell line, and mark them as "overlapping" (they'll later have an isTrue(iT) value of 1) or "non-overlapping" (iT = 0).

The model learns by becoming better at predicting if a putative site isTrue or not. The genomic and ATACseq data at each putative site are fed through the network and at the output the network gives us two values, predicting how well the site fits label 0, or label 1. We will then test the model on the 'test' subset of putative sites that we kept aside. This test gives us an expected error rate of the model. We can now apply the model on datasets for which we do not have available iT values (in other words, we don't have a CHIPseq experiment of the factor available), and expect our results to be as good as they were on the test set.

Data subsets

We split our full set of data into three subsets:

• Test set
  This part of the data is only used after training to evaluate how well the model is trained, the model will never encounter this dataset during training. We set aside 10% of our data for this.

• Validation set
  This part of the data is used during training to evaluate the progress of the training. During training, we test our model with this dataset to determine when the model has stopped improving. We set another 10% of our total data for this.

• Training set
  This is the data given to our model for training, after setting aside 10% each for test and validation, 80% of the datapoints were left for training

These datasets tend to be very unbalanced, in our cases the "bad" sites outnumber the good ones with ratios from 1:3 to 1:10 or worse. To mitigate the effect of this imbalance, we make sure that the ratio of good/bad sites is maintained (as much as possible) in all subsets.

Preparing your data

Let us assume that you are looking to train your model with data from two contexts where you have CHIPseq and use it to infer good TFbs in a third context where CHIP data is not available.

Prepare the PWM hits for training

You will need to prepare a tsv (tab-separated) file that looks like this (including the header with exactly the following names):

nimrod prep dataset (1/4)

This subcommand will help you prepare this dataset. It assumes you have a set of PWM hits on the whole genome, from those, it keeps only those that overlap an ATAC seq peak and based on the chipseq peaks, it marks the retained sites as true (iT=1) if they overlap a chipseq peak or false (iT=0) if they don't. Furthermore, it divides the dataset into train/validation/test subsets. You can choose to only make a test dataset by passing the --onlytest option.

As an example, here's how I used the subcommand to prepare a dataset on xenopus tropicalis to apply a pretrained model on it. We want to get a set of high-confidense PWM hits from nimrod, by using our ATACseq data from a developmental stage. We decided to test any PWM hits that fall inside an ATACseq peak (would you trust a putative hit that is not inside an atac? maybe so for other TFs but in our experience, CTCF sites tend to be found inside atac peaks).

We need to make a simple json file that points nimrod to our data. We only have one biological context (one developmental stage) in this dataset, and no chipseq data, and we will only be 'testing' this data, not training with it. Therefore, our json looks like this:

{
"chip_peaks" : {},                  
"atac_peaks" : {"stage1": "/path/to/atac_seq_peaks.bed"},
"genomic_pwm" : "/path/to/CTCF_M1957_genomichits_onXenTro9_p0.0005.bed"
}

I had previously created the genomic hits with moods:

moods_dna.py -p 0.0005 -S CTCF_M1957.pfm -o CTCF_M1957_onXenTro9_p0.0005.csv -s xenTro9.fa 

I converted the csv moods output to a bedfile with a small python script:

csvfp = "CTCF_M1957_onXenTro9_p0.0005.csv"
outfp = "CTCF_M1957_genomichits_onXenTro9_p0.0005.bed"
tfname = 'CTCF_M1957'
tflen = 19 #the number of bases in the m1957 pfm

thecsv = pd.read_csv(csvfp,sep=',', header=None,comment='#')
thecsv['end'] = thecsv[2] + tflen
thecsv['name'] = tfname
thecsv[[0,2,'end','name',4,3]].to_csv(outfp, sep='\t', header=False, index=False)

You can find this function in nimrod's auxiliary functions: "nimrod aux moods --help"

So to make the dataset:

>> nimrod prep dataset --json ./theonewejustmade.json --out ./xentro9_CTCF_M19_fornimrod.tsv --nochip

The --nochip option is important, it tells nimrod that it doesn't need to work with chipseq peaks and to mark everything for 'test'.

If this subcommand doesnn't fit your needs, you'll need to manually create this tsv file

chrom, start, end, score, strand

These are the typical fields of a bed file, the chromosome, start end and score of a PWM
The start and end refer to the start and end of the PWM itself, the model will extend this to 1000bp automatically. In actuallity, the PWM score field is not needed but it's usefull to keep it there for later analyses

stage

This is the name of your biological context

iT

stands for "is True", you will have to mark the 'good' sites with 1 and the bad with 0. As mentioned above, we marked the sites overlapping a chipseq peak with 1 and the rest with 0

set

This allows the model to split the data into the three subsets that were discussed above

Further data preparation

Nimrod comes with a number of subcommands to help you prepare the data and run the training and testing. Three steps are necessary :

nimrod prep genome (2/4)

This creates a pickled version of a typical fasta genome that is much faster for nimrod to load from the disk, you only need to do this once for any genome you might use.

# example
>> nimrod prep genome --fa hg19.fa --out hg19.pickle

nimrod prep bigwig (3/4)

Besides a genome, nimrod also requires the ATAC signal from each condition. This command constructs a series of big wig files and an index from your .bam files. You should run this once for each collection of bam files that you intend to use with the model.

You will need:

• a json file matching your biological contexts to their bam files.
  This is a simple text file, you can have nimrod print an example of such a file with

>> nimrod prep bigwig --printjson > nimrod_seed_example.json

Here is what the final json could look like:

{
"cell_line_a": ["/path/to/cell_line_a.bam"],

"cell_line_b":    ["/path/to/cell_line_b_1.bam",
                        "/path/to/cell_line_b_2.bam",
                        "/path/to/cell_line_b_3.bam"]
}

Make sure the name of the contexts here ("cell_line_a", etc) exactly matches the names used in the PWM tsv table.

• A chrominfo file
  This is a standard bioinformatics text file containing the names of the chromosomes in your orgnanisms's genome and their respective size. Make sure the chromosome names here are the same as the chromosome names in your bam file. Here are the top three lines of our hg19 chrominfo file:

chr1    249250621
chr2    243199373
chr3    198022430

Once you have the chrominfo file and the json file, you can have nimrod prepare the BW folder:

# example:
>> nimrod prep bigwig --procs 16 --json human_lines.json --chrominfo hg19.txt --out human_lines_nimrod_bw

The folder will eventually contain a number of .bw files as well as a new json file that works as an index. If you move the folder around and the paths of the files change, nimrod will not be able to find them, but you can just edit the index file manually to fix the faulty paths.

nimrod prep ndata (4/4)

This extracts the ATACseq data from around the PWM sites that you are about to analyze and packs it in a very fast loading format.

#example:
>> nimrod prep ndata --genome hg19.pickle --bwf human_lines_nimrod_bw --df CTCF_sites_human_lines.tsv --out nimrod_preproc_hg_lines

Having completed the final step (nimrod prep ndata) you can now train the model.

nimrod train

This runs the training of the network and creates a trained model that you can later use to evaluate unknown data.

Necessary arguments

• genome
  Path to a genome file (.fa or .pickle)

• df
  path to a table with the motif hits, stage, isReal, set

• ndata
  path to folder created with the 'makechunks' command

• out
  path to output folder

Optional arguments

• epochs
  The model will go over the entire training set more than one times. Each such time is called an epoch. This option defines a maximum number of epochs to run (default: 300)

• valfre
  Validation Frequency: Every X number of batches, the model tests itself against the validation subset. This argument defines this validation frequency (default: 5000).

• patience
  If we don't improve after this many validation tests, stop the training (default: 10)

• slr
  Starting Learning Rate. (default: 1e-05)

• dexe
  The learning rate decays to (learning rate)0.9 every dexe*epochs. If you set it to something very high, like (10000), the model will run with practically unchanged learning rate. (default: 5)

• bsm
  Batch size multiplier, sets the size of each batch that the model examines. Try values from 2-50+ (default: 2). The 'nimrod aux bsm_analysis' can help you make a decision. Run the analysis and choose the smallest bsm value that gives a 'ratio' value (the fourth column of the analysis output) as close to 0 as possible.

Example training command:

# example:
>> nimrod train --epochs 10 --dexe 1000 --bsm 10 --ndata nimrod_preproc_1 --genome hg19.pickle --df CTCF_sites_human_lines.tsv --out nimrod_training_1

training output

After a few seconds of loading, and some warnings from tensorflow, nimrod will start printing its own output:

The most important of these columns is probably the 7th one, where the AUC score of the validation set is printed. You expect to see the value of this column to start at 0.5 or higher and continue rising. If after 3-4 lines of output the auc value is still around 0.5, the model is not learning, stop it and try again with different parameters (learning rate, bsm).

The fourth column outputs 1 in the cases where the model did better than before and will thus be saved.

In the folder created by nimrod for training, you will find the following. These are quite redundant so exercise your judgment to remove those you don't need.

• arguments_overview.txt An overview of the training arguments you used, because let's face it, you already forgot what values you used.

• df_test.tsv df_train.tsv df_valid.tsv The tsv file that was used as input, split into the three subsets, and with the regions expanded to 1kb. Note that some hits might have been removed.

• full_motif_dataframe.tsv The dataframe you used as input, just to make sure you don't lose it.

• the_model.py The exact TF model that was used in this run

• training.log A copy of the output that was printed during the training

• test_results/ A csv file for every validation test that was run during training

• model_states/ The saved, trained models of nimrod from this run. The last two best states according to AUC test, the last two according to loss and the very final state. Your best model will typically be the lat auc-based model but YMMV.

nimrod continue-train

This allows you to continue training on a previously created model.

nimrod test

This allows you to use a previously trained model on unknown data. We would use this to test the model on the test subset of data at the end of training.

Acknowledgment

The project's avatar is a detail of a work from Assaf Hinden, on the statue of Nimrod from Yitzhak Danziger.