### Install Basenji with Pip and Setuptools

Source: https://github.com/calico/basenji/blob/master/README.md

Installs Basenji and its dependencies using pip and setuptools. TensorFlow should be installed separately.

```bash
python setup.py develop
pip install tensorflow (or tensorflow-gpu)
```

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

### Setup Output Directory and Write Statistics

Source: https://github.com/calico/basenji/blob/master/jupyter/saluki_data.ipynb

Initializes the output directory by removing it if it exists and creating it anew. It then calculates and saves dataset statistics, including the number of folds and sequences per fold, to a JSON file.

```python
import os
import shutil
import json

# output directory

out_dir = 'data'
if os.path.isdir(out_dir):
    shutil.rmtree(out_dir)
os.mkdir(out_dir)

################################################################
# stats
stats_dict = {}
stats_dict['num_targets'] = 1
stats_dict['target_length'] = 1
stats_dict['length_t'] = full_len_t

stats_dict['num_folds'] = num_folds
for fi in range(num_folds):
    stats_dict['fold%d_seqs' % fi] = int(fold_masks[fi].sum())

with open('%s/statistics.json' % out_dir, 'w') as stats_json_out:
    json.dump(stats_dict, stats_json_out, indent=4)
```

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

### Install cooltools

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/README.md

Installs the cooltools Python package, a dependency for exploring Akita model predictions.

```bash
pip install cooltools
```

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

### Install Basenji with Anaconda

Source: https://github.com/calico/basenji/blob/master/README.md

Installs dependencies and Basenji using Anaconda. Ensure TensorFlow is installed separately.

```bash
conda env create -f environment.yml
conda install tensorflow (or tensorflow-gpu)
python setup.py develop --no-deps
```

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

### Verify Basenji Installation

Source: https://github.com/calico/basenji/blob/master/README.md

Verifies the Basenji installation by importing the library in a Python interpreter.

```python
import basenji
```

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

### Define Mouse Genome and Annotation File Paths (Commented Out)

Source: https://github.com/calico/basenji/blob/master/jupyter/saluki_data.ipynb

Provides commented-out example file paths for mouse genome assembly (FASTA) and gene annotation (GTF).

```python
# # mouse
# stability_file = 'mouse_stability.tsv'
# vikram_gtf_file = '/home/vagar/RNAbinding/forPaper/datasets/mouse/Mus_musculus.GRCm3.90.chosenTranscript.gtf.gz'
# genome_fasta_file = '%s/assembly/ensembl/Mus_musculus.GRCm38.dna.primary_assembly.fa' % os.environ['MM10']
```

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

### Install Basenji with Pre-specified Anaconda Environment

Source: https://github.com/calico/basenji/blob/master/README.md

Installs Basenji using a fully pre-specified Anaconda environment for guaranteed working dependency versions. TensorFlow must be installed separately.

```bash
conda env create -f prespecified.yml
conda install tensorflow (or tensorflow-gpu)
python setup.py develop --no-deps
```

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

### GPU Memory Requirement for Training

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Training Basenji models with default parameters requires a GPU with more than 12GB of RAM. Ensure your hardware meets this specification before starting.

```text
note that training with default parameters requires GPU with >12Gb RAM
```

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

### Set Basenji Environment Variables

Source: https://github.com/calico/basenji/blob/master/README.md

Sets essential environment variables for Basenji after installation. These variables define the Basenji directory and add its bin directory to the PATH and PYTHONPATH.

```bash
export BASENJIDIR=~/code/Basenji
export PATH=$BASENJIDIR/bin:$PATH
export PYTHONPATH=$BASENJIDIR/bin:$PYTHONPATH
```

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

### Find Indices for Example Sequences

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/explore_model.ipynb

Identifies the indices of specific genomic regions within the test dataset. This is useful for retrieving specific examples for visualization or analysis.

```python
fig2_examples = [   'chr12:115163136-116211712',
                    'chr11:75429888-76478464',
                    'chr15:63281152-64329728' ]
fig2_inds = []
for seq in fig2_examples:
    print(seq)
    chrm,start,stop = seq.split(':')[0], seq.split(':')[1].split('-')[0], seq.split(':')[1].split('-')[1]
    test_ind = np.where( (sequences_test['chr'].values== chrm) * 
                         (sequences_test['start'].values== int(start))* 
                         (sequences_test['stop'].values==  int(stop ))  )[0][0]
    fig2_inds.append(test_ind)
fig2_inds
```

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

### Example Gene and Transcript Data

Source: https://github.com/calico/basenji/blob/master/jupyter/saluki_data.ipynb

Represents processed gene and transcript data, including gene ID, transcript ID, saluki score, and genomic coordinates for exons.

```text
1 + gene_id "ENSG00000122484"; transcript_id "ENST00000610020"; saluki_score "1.0";  [92298965-92299146],[92300194-92300239],[92301476-92301590],[92303977-92304075],[92304284-92304349],[92307188-92307276],[92320599-92320634],[92323445-92324375],[92333391-92333473],[92336347-92336427],[92345846-92345914],[92380724-92380873],[92387011-92402056]
```

```text
3 + gene_id "ENSG00000168314"; transcript_id "ENST00000311042"; saluki_score "1.0";  [39467680-39467740],[39480040-39480123],[39502066-39502275],[39502463-39502876],[39513383-39516001]
```

```text
3 + gene_id "ENSG00000185008"; transcript_id "ENST00000332191"; saluki_score "1.0";  [77040730-77040846],[77098014-77098340],[77477414-77477571],[77481099-77481219],[77493244-77493382],[77522775-77522902],[77546338-77546462],[77550818-77550989],[77557944-77558149],[77562651-77562732],[77563167-77563329],[77564954-77565120],[77568313-77568434],[77574499-77574730],[77577490-77577614],[77579947-77580118],[77588751-77588933],[77595142-77595184],[77596623-77596750],[77602210-77602491],[77607798-77607954],[77617513-77617773],[77622227-77622432],[77632474-77632656],[77634870-77635043],[77644704-77644904],[77646054-77647116]
```

```text
3 - gene_id "ENSG00000174963"; transcript_id "ENST00000525172"; saluki_score "1.0";  [147386048-147388859],[147390931-147391246],[147395852-147396469],[147402728-147402812],[147403964-147404284]
```

```text
4 + gene_id "ENSG00000126545"; transcript_id "ENST00000246891"; saluki_score "1.0";  [69931081-69931117],[69932544-69932606],[69934212-69934244],[69934690-69934710],[69935926-69935949],[69936456-69936479],[69936566-69936607],[69937121-69937144],[69937800-69937823],[69939176-69939208],[69940021-69940044],[69941019-69941060],[69942046-69942063],[69942536-69942577],[69944850-69945004],[69946196-69946574]
```

```text
4 - gene_id "ENSG00000164096"; transcript_id "ENST00000399075"; saluki_score "1.0";  [119296419-119298779],[119300336-119300733],[119304245-119304445]
```

```text
4 + gene_id "ENSG00000186352"; transcript_id "ENST00000335174"; saluki_score "1.0";  [185396484-185396950],[185397150-185397302],[185398937-185399028],[185399570-185399773],[185400017-185400241]
```

```text
7 - gene_id "ENSG00000048405"; transcript_id "ENST00000393313"; saluki_score "1.0";  [127370043-127371814],[127373342-127375034],[127377186-127377329],[127386060-127386155],[127391497-127391615],[127392158-127392691]
```

```text
X + gene_id "ENSG00000122824"; transcript_id "ENST00000376006"; saluki_score "1.0";  [51332231-51332367],[51332883-51333459],[51336239-51337525]
```

```text
X - gene_id "ENSG00000196368"; transcript_id "ENST00000375992"; saluki_score "1.0";  [51490011-51491749],[51495951-51496596]
```

```text
8 + gene_id "ENSG00000175445"; transcript_id "ENST00000311322"; saluki_score "1.0";  [19938971-19939528],[19948180-19948340],[19951769-19951948],[19953310-19953421],[19954120-19954353],[19955841-19956083],[19959260-19959380],[19960901-19961083],[19962115-19962219],[19965310-19967258]
```

```text
8 + gene_id "ENSG00000069206"; transcript_id "ENST00000380789"; saluki_score "1.0";  [24441084-24441160],[24442473-24442576],[24447186-24447262],[24463882-24463960],[24465699-24465775],[24466799-24466988],[24468767-24468820],[24476433-24476504],[24482142-24482311],[24485277-24485361],[24487187-24487317],[24489159-24489333],[24490799-24490888],[24491903-24492098],[24492495-24492597],[24493043-24493229],[24499236-24499316],[24500178-24500256],[24500790-24500895],[24501477-24501576],[24506091-24506156],[24507480-24507535],[24508546-24509548]
```

```text
8 + gene_id "ENSG00000104435"; transcript_id "ENST00000518111"; saluki_score "1.0";  [79611152-79611214],[79636802-79636897],[79641378-79641550],[79654871-79655062],[79663580-79663661],[79664815-79665007]
```

```text
9 - gene_id "ENSG00000122756"; transcript_id "ENST00000378980"; saluki_score "1.0";  [34551432-34552071],[34552161-34552329],[34552674-34552854],[34556255-34556418],[34557526-34557692],[34557867-34557984],[34564599-34564832],[34568897-34568981],[34581095-34581205],[34589555-34589737]
```

```text
10 - gene_id "ENSG00000203791"; transcript_id "ENST00000368836"; saluki_score "1.0";  [124759954-124760503],[124762300-124762558],[124765392-124765608],[124774675-124774782],[124789043-124789157],[124790273-124790338],[124791724-124791870]
```

```text
13 - gene_id "ENSG00000165480"; transcript_id "ENST00000314759"; saluki_score "1.0";  [21153595-21155150],[21155693-21155811],[21157922-21158125],[21159902-21159987],[21161790-21161875],[21167988-21168399],[21172339-21172504],[21172620-21172681],[21176375-21176602]
```

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

### Generate TFRecords for Basenji Training Data

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Use `basenji_data_write.py` to split genomic data into tfrecords for training, validation, and testing. Specify start and end coordinates, input FASTA, sequences BED, coverage BED, and output tfrecord path. For production, omit the `--sample` option.

```bash
basenji_data_write.py -s 0 -e 128 ./data/hg38.ml.fa ./data/1m/sequences.bed ./data/1m/seqs_cov ./data/1m/tfrecords/train-0.tfr
basenji_data_write.py -s 128 -e 256 ./data/hg38.ml.fa ./data/1m/sequences.bed ./data/1m/seqs_cov ./data/1m/tfrecords/train-1.tfr
basenji_data_write.py -s 256 -e 350 ./data/hg38.ml.fa ./data/1m/sequences.bed ./data/1m/seqs_cov ./data/1m/tfrecords/train-2.tfr
basenji_data_write.py -s 350 -e 478 ./data/hg38.ml.fa ./data/1m/sequences.bed ./data/1m/seqs_cov ./data/1m/tfrecords/valid-0.tfr
basenji_data_write.py -s 478 -e 606 ./data/hg38.ml.fa ./data/1m/sequences.bed ./data/1m/seqs_cov ./data/1m/tfrecords/valid-1.tfr
basenji_data_write.py -s 606 -e 667 ./data/hg38.ml.fa ./data/1m/sequences.bed ./data/1m/seqs_cov ./data/1m/tfrecords/valid-2.tfr
basenji_data_write.py -s 667 -e 795 ./data/hg38.ml.fa ./data/1m/sequences.bed ./data/1m/seqs_cov ./data/1m/tfrecords/test-0.tfr
basenji_data_write.py -s 795 -e 923 ./data/hg38.ml.fa ./data/1m/sequences.bed ./data/1m/seqs_cov ./data/1m/tfrecords/test-1.tfr
basenji_data_write.py -s 923 -e 981 ./data/hg38.ml.fa ./data/1m/sequences.bed ./data/1m/seqs_cov ./data/1m/tfrecords/test-2.tfr
```

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

### Create Samples Table

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Generates a tab-separated file named 'microc_cools.txt' in the 'data' directory. This file lists the cooler files and their associated metadata, including index, identifier, file path, clip value, sum statistic, and description.

```python
lines = [['index','identifier','file','clip','sum_stat','description']]
lines.append(['0', 'HFF', './data/coolers/HFF_hg38_4DNFIP5EUOFX.mapq_30.2048.cool', '2', 'sum', 'HFF'])
lines.append(['1', 'H1hESC', './data/coolers/H1hESC_hg38_4DNFI1O6IL1Q.mapq_30.2048.cool', '2', 'sum', 'H1hESC'])

samples_out = open('data/microc_cools.txt', 'w')
for line in lines:
    print('\t'.join(line), file=samples_out)
samples_out.close()
```

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

### Download and Extract Tutorial Data

Source: https://github.com/calico/basenji/blob/master/tutorials/train_test.ipynb

This snippet downloads and extracts a small HDF5 dataset for tutorial purposes if it doesn't already exist. It uses curl and tar for data handling.

```python
if len(glob.glob('data/heart_l131k/tfrecords/*.tfr')) == 0:
    subprocess.call('curl -o data/heart_l131k.tgz https://storage.googleapis.com/basenji_tutorial_data/heart_l131k.tgz', shell=True)
    subprocess.call('tar -xzvf data/heart_l131k.tgz', shell=True)
```

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

### Download Micro-C Cooler Files

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Creates a directory for cooler files if it doesn't exist and downloads two specific Micro-C cooler files (HFF and H1hESC) if they are not already present.

```python
if not os.path.exists('./data/coolers'):
    os.mkdir('./data/coolers')
if not os.path.isfile('./data/coolers/HFF_hg38_4DNFIP5EUOFX.mapq_30.2048.cool'):
    subprocess.call('curl -o ./data/coolers/HFF_hg38_4DNFIP5EUOFX.mapq_30.2048.cool'+ 
            ' https://storage.googleapis.com/basenji_hic/tutorials/coolers/HFF_hg38_4DNFIP5EUOFX.mapq_30.2048.cool', shell=True)
    subprocess.call('curl -o ./data/coolers/H1hESC_hg38_4DNFI1O6IL1Q.mapq_30.2048.cool'+ 
            ' https://storage.googleapis.com/basenji_hic/tutorials/coolers/H1hESC_hg38_4DNFI1O6IL1Q.mapq_30.2048.cool', shell=True)
```

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

### List Cooler Files

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Lists the contents of the data/coolers directory to verify the downloaded cooler files.

```bash
ls ./data/coolers/
```

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

### Import necessary libraries

Source: https://github.com/calico/basenji/blob/master/tutorials/sad.ipynb

Imports the required Python libraries for the tutorial.

```python
import h5py
import os
import subprocess
```

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

### Get Data Dimensions

Source: https://github.com/calico/basenji/blob/master/jupyter/basenji_sat_vcf_viz.ipynb

Retrieves the number of sequences and targets from the loaded HDF5 file and targets DataFrame.

```python
num_seqs, seq_len, _ = sat_h5['seqs'].shape
num_targets = targets_df.shape[0]
```

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

### Train Akita Model

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Use this command to train an Akita model. Ensure you have sufficient GPU RAM (>12Gb) for training with default parameters. Specify output directory and parameter file.

```bash
! akita_train.py -k -o ./data/1m/train_out/  ./data/1m/params_tutorial.json ./data/1m/
```

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

### Load Data Parameters and Targets

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/explore_model.ipynb

Loads target identifiers, file mappings, and data statistics including sequence length, target length, and cropping information. Assumes data is in './data/'.

```python
data_dir =   './data/'

hic_targets = pd.read_csv(data_dir+'/targets.txt',sep='	')
hic_file_dict_num = dict(zip(hic_targets['index'].values, hic_targets['file'].values) )
hic_file_dict     = dict(zip(hic_targets['identifier'].values, hic_targets['file'].values) )
hic_num_to_name_dict = dict(zip(hic_targets['index'].values, hic_targets['identifier'].values) )

# read data parameters
data_stats_file = '%s/statistics.json' % data_dir
with open(data_stats_file) as data_stats_open:
    data_stats = json.load(data_stats_open)
seq_length = data_stats['seq_length']
target_length = data_stats['target_length']
hic_diags =  data_stats['diagonal_offset']
target_crop = data_stats['crop_bp'] // data_stats['pool_width']
target_length1 = data_stats['seq_length'] // data_stats['pool_width']
```

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

### Import Necessary Libraries

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/explore_model.ipynb

Imports core libraries for data manipulation, model building, and visualization. Sets the environment to use the CPU for execution.

```python
import os
import json
import subprocess
os.environ["CUDA_VISIBLE_DEVICES"] = '-1' ### run on CPU

from cooltools.lib.numutils import set_diag
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pysam
import tensorflow as tf

from basenji import dataset, dna_io, seqnn
```

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

### Write Data to TFRecord Files

Source: https://github.com/calico/basenji/blob/master/jupyter/saluki_data.ipynb

Writes gene sequences and their associated targets into TFRecord files, organized by fold. This function handles slicing data, creating features, and writing TFRecord examples with ZLIB compression. It uses helper functions like `tfr_feature_dna` and `tfr_feature_targets`.

```python
def write_fold(fold_index, out_dir, seqs_per_tfr=128):
    tf_opts = tf.io.TFRecordOptions(compression_type='ZLIB')
    
    # slice fold's genes
    fold_df = stability_df[fold_masks[fold_index]]
    
    gi = 0
    num_tfr = int(np.ceil(fold_df.shape[0] / seqs_per_tfr))
    for fi in range(num_tfr):
        tfr_file = '%s/fold%d-%d.tfr' % (out_dir,fold_index,fi)
        with tf.io.TFRecordWriter(tfr_file, tf_opts) as writer:
            for ti in range(seqs_per_tfr):
                if gi < fold_df.shape[0]:
                    sgene = fold_df.iloc[gi]
                    gene_id = sgene.Gene
                    
                    # trim sequence
                    seq1 = gene_seq1[gene_id][:full_len_t]
                    seq_len = seq1.shape[0]
                    
                    # split
                    sequence = seq1[:,:4]
                    coding = seq1[:,4]
                    splice = seq1[:,5]
                        
                    # targets
                    # targets = np.log2(sgene.Stability + pseudocount)
                    targets = sgene.Stability

                    # hash to bytes
                    features_dict = {
                        'lengths': tfr_feature_ints([seq_len]),
                        'sequence': tfr_feature_dna(sequence),
                        'coding': tfr_feature_track(coding),
                        'splice': tfr_feature_track(splice),
                        'targets': tfr_feature_targets(targets)
                    }

                    # write example
                    example = tf.train.Example(features=tf.train.Features(feature=features_dict))
                    writer.write(example.SerializeToString())
                
                    # next gene
                    gi += 1
```

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

### Display First Three Lines of Sequences BED File

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Use `head -n3` to display the first three lines of the `sequences.bed` file, showing the format and content of sequence entries.

```bash
head -n3 data/1m/sequences.bed
```

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

### Display Accuracy Metrics

Source: https://github.com/calico/basenji/blob/master/tutorials/train_test.ipynb

Use this command to display the accuracy metrics (Pearson correlation and R2) for each dataset after training. The output is typically found in an 'acc.txt' file.

```bash
! cat output/heart_test/acc.txt
```

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

### Load Model Parameters and Initialize SeqNN

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/explore_model.ipynb

Loads model parameters from a JSON file and initializes a Basenji SeqNN model. Ensure 'params.json' and 'model_best.h5' are in the specified model directory.

```python
model_dir = './'
params_file = model_dir+'params.json'
model_file  = model_dir+'model_best.h5'
with open(params_file) as params_open:
    params = json.load(params_open)
    params_model = params['model']
    params_train = params['train']

seqnn_model = seqnn.SeqNN(params_model)
```

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

### Train Basset-style Model

Source: https://github.com/calico/basenji/blob/master/manuscripts/basset/README.md

Command to train a Basset-style model using specified parameters and dataset. Ensure the data_basset directory contains the generated annotation dataset.

```bash
basenji_train.py -k -o train_basset params_basset.json data_basset
```

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

### Download and Prepare Genome FASTA

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/explore_model.ipynb

Downloads the hg38 genome FASTA file if it doesn't exist and prepares it for sequence fetching using pysam. Requires 'os', 'subprocess', and 'pysam' libraries.

```python
import os
import subprocess
import pysam

if not os.path.isfile('./data/hg38.ml.fa'):
    print('downloading hg38.ml.fa')
    subprocess.call('curl -o ./data/hg38.ml.fa.gz https://storage.googleapis.com/basenji_barnyard2/hg38.ml.fa.gz', shell=True)
    subprocess.call('gunzip ./data/hg38.ml.fa.gz', shell=True)

fasta_open = pysam.Fastafile('./data/hg38.ml.fa')
```

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

### Initialize Gene and Genome Data Structures

Source: https://github.com/calico/basenji/blob/master/jupyter/saluki_data.ipynb

Initializes data structures for gene sequences and opens a GTF file and a genome FASTA file for reading. Requires genes_gtf_file and genome_fasta_file to be defined.

```python
########################################################
# read gene data

# open genes GTF
genes_gtf = pygene.GTF(genes_gtf_file)

# open genome FASTA
genome_open = pysam.Fastafile(genome_fasta_file)

# initialize dicts
gene_seq1 = {}
```

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

### Display First 3 Lines of BED File

Source: https://github.com/calico/basenji/blob/master/tutorials/preprocess.ipynb

Shows the first three lines of the sequences BED file. Helps in understanding the file format and content.

```python
! head -n3 data/heart_l131k/sequences.bed
```

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

### Create Samples Table for BigWig Files

Source: https://github.com/calico/basenji/blob/master/tutorials/preprocess.ipynb

Generates a tab-separated samples table ('heart_wigs.txt') from provided BigWig file paths and metadata. This table is used by `basenji_data.py` to identify target data.

```python
lines = [['index','identifier','file','clip','sum_stat','description']]
lines.append(['0', 'CNhs11760', 'data/CNhs11760.bw', '384', 'sum', 'aorta'])
lines.append(['1', 'CNhs12843', 'data/CNhs12843.bw', '384', 'sum', 'artery'])
lines.append(['2', 'CNhs12856', 'data/CNhs12856.bw', '384', 'sum', 'pulmonic_valve'])

samples_out = open('data/heart_wigs.txt', 'w')
for line in lines:
    print('\t'.join(line), file=samples_out)
samples_out.close()
```

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

### Download Trained Model

Source: https://github.com/calico/basenji/blob/master/tutorials/sat_mut.ipynb

Downloads a pre-trained Basenji model ('model_best.h5') and creates the necessary directory if they do not already exist. This model is used for predictions.

```python
if not os.path.isdir('models/heart'):
    os.mkdir('models/heart')
if not os.path.isfile('models/heart/model_best.h5'):
    subprocess.call('curl -o models/heart/model_best.h5 https://storage.googleapis.com/basenji_tutorial_data/model_best.h5', shell=True)
```

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

### Load Test Dataset

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/explore_model.ipynb

Loads the test dataset, including sequences and their corresponding targets and inputs. Requires the 'get_data.sh' script to be run first. The test data is structured into sequences, inputs (1-hot encoded), and targets (log(obs/exp) data).

```python
# note: run %%bash get_data.sh 
# if you have not already downloaded the data

sequences = pd.read_csv(data_dir+'sequences.bed', sep='	', names=['chr','start','stop','type'])
sequences_test = sequences.iloc[sequences['type'].values=='test']
sequences_test.reset_index(inplace=True, drop=True)

test_data = dataset.SeqDataset(data_dir, 'test', batch_size=8)

# test_targets is a float array with shape 
# [#regions, #pixels, target #target datasets]
# representing log(obs/exp)data, where #pixels 
# corresponds to the number of entries in the flattened
# upper-triangular representation of the matrix

# test_inputs are 1-hot encoded arrays with shape
# [#regions, 2^20 bp, 4 nucleotides datasets]

test_inputs, test_targets = test_data.numpy(return_inputs=True, return_outputs=True)
```

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

### Load HDF5 File and List Keys

Source: https://github.com/calico/basenji/blob/master/tutorials/sad.ipynb

Opens an HDF5 file containing SNP activity difference predictions and lists all the keys (datasets) within it. This is useful for understanding the structure of the output file.

```python
sad_h5 = h5py.File('output/rfx6_sad/sad.h5', 'r')
list(sad_h5.keys())
```

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

### Download HG19 FASTA and Index

Source: https://github.com/calico/basenji/blob/master/tutorials/preprocess.ipynb

Downloads a simplified HG19 FASTA file and its index if they do not already exist. This is a prerequisite for further data processing.

```python
import os, subprocess

if not os.path.isfile('data/hg19.ml.fa'):
    subprocess.call('curl -o data/hg19.ml.fa https://storage.googleapis.com/basenji_tutorial_data/hg19.ml.fa', shell=True)
    subprocess.call('curl -o data/hg19.ml.fa.fai https://storage.googleapis.com/basenji_tutorial_data/hg19.ml.fa.fai', shell=True)
```

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

### Read and Bin Target Data

Source: https://github.com/calico/basenji/blob/master/jupyter/basenji_sad_viz.ipynb

Reads target data for the determined region from a W5 file, reshapes it into bins, and clips the values. Displays the shape of the binned target data.

```python
# read targets

with h5py.File(targets_df.iloc[ti].file,'r') as target_w5:
    region_targets = target_w5[chrm][start:end]
    
# sum in bins
region_targets = np.reshape(region_targets, (-1,128)).sum(axis=-1)

# clip
region_targets = np.clip(region_targets, 0, 384)
region_targets.shape
```

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

### Run Basenji Training

Source: https://github.com/calico/basenji/blob/master/tutorials/train_test.ipynb

Execute the basenji_train.py script to train a model. This command requires specifying output directories, parameter files, and the data directory. Training can take several hours depending on hardware.

```python
! basenji_train.py -o models/heart models/params_small.json data/heart_l131k
```

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

### Run Saturation Mutagenesis with BED Input

Source: https://github.com/calico/basenji/blob/master/tutorials/sat_mut.ipynb

Executes the `basenji_sat_bed.py` script to perform saturation mutagenesis. This command requires the genome FASTA, output directory, target WIG file list, model parameters, trained model, and a BED file of regions.

```bash
basenji_sat_bed.py -f data/hg19.ml.fa -l 200 -o output/gata4_sat --rc -t data/heart_wigs.txt models/params_small.json models/heart/model_best.h5 data/gata4.bed
```

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

### List Output Files

Source: https://github.com/calico/basenji/blob/master/tutorials/sat_mut.ipynb

Use this command to list the generated plot files in the specified output directory. Ensure the directory exists and contains the plot images.

```python
! ls output/gata4_sat/plots
```

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

### Akita Model Architecture Summary

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

This output details the layers, output shapes, and parameter counts of the Akita model. It is generated after the training command is executed.

```text
Model: "model_1"
__________________________________________________________________________________________________
 Layer (type)                Output Shape                 Param #   Connected to                  
==================================================================================================
 sequence (InputLayer)       [(None, 1048576, 4)]         0         []                            
                                                                                                  
 stochastic_reverse_complem  ((None, 1048576, 4),         0         ['sequence[0][0]']            
 ent (StochasticReverseComp   ())
 lement)
                                                                                                  
 stochastic_shift (Stochast  (None, 1048576, 4)           0         ['stochastic_reverse_complemen
 icShift)                                                           t[0][0]']                     
                                                                                                  
 re_lu (ReLU)                (None, 1048576, 4)           0         ['stochastic_shift[0][0]']    
                                                                                                  
 conv1d (Conv1D)             (None, 1048576, 96)          4224      ['re_lu[0][0]']               
                                                                                                  
 batch_normalization (Batch  (None, 1048576, 96)          384       ['conv1d[0][0]']              
 Normalization)
                                                                                                  
 max_pooling1d (MaxPooling1  (None, 524288, 96)           0         ['batch_normalization[0][0]'] 
 D)
                                                                                               
                                                                                                  
 re_lu_1 (ReLU)              (None, 524288, 96)           0         ['max_pooling1d[0][0]']       
                                                                                                  
 conv1d_1 (Conv1D)           (None, 524288, 96)           46080     ['re_lu_1[0][0]']             
                                                                                                  
 batch_normalization_1 (Bat  (None, 524288, 96)           384       ['conv1d_1[0][0]']            
 chNormalization)
                                                                                                                
 max_pooling1d_1 (MaxPoolin  (None, 262144, 96)           0         ['batch_normalization_1[0][0]'
 g1D)                                                               ]                             
                                                                                                  
 re_lu_2 (ReLU)              (None, 262144, 96)           0         ['max_pooling1d_1[0][0]']     
                                                                                                  
 conv1d_2 (Conv1D)           (None, 262144, 96)           46080     ['re_lu_2[0][0]']             
                                                                                                  
 batch_normalization_2 (Bat  (None, 262144, 96)           384       ['conv1d_2[0][0]']            
 chNormalization)
                                                                                                                
 max_pooling1d_2 (MaxPoolin  (None, 131072, 96)           0         ['batch_normalization_2[0][0]'
 g1D)                                                               ]                             
                                                                                                  
 re_lu_3 (ReLU)              (None, 131072, 96)           0         ['max_pooling1d_2[0][0]']     
                                                                                                  
 conv1d_3 (Conv1D)           (None, 131072, 96)           46080     ['re_lu_3[0][0]']             
                                                                                                  
 batch_normalization_3 (Bat  (None, 131072, 96)           384       ['conv1d_3[0][0]']            
 chNormalization)

```

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

### Log Device Placement Output

Source: https://github.com/calico/basenji/blob/master/docs/train.md

This output shows how TensorFlow assigns operations to devices. In this specific case, training occurs on the CPU, indicating that the GPU might be unsupported or not utilized.

```text
Device mapping:
...
2017-07-23 12:31:25.796354: I tensorflow/core/common_runtime/simple_placer.cc:847] cnn1/BatchNorm/Const: (Const)/job:localhost/replica:0/task:0/cpu:0
cnn0/BatchNorm/Const_1: (Const): /job:localhost/replica:0/task:0/cpu:0
2017-07-23 12:31:25.796361: I tensorflow/core/common_runtime/simple_placer.cc:847] cnn0/BatchNorm/Const_1: (Const)/job:localhost/replica:0/task:0/cpu:0
cnn0/BatchNorm/Const: (Const): /job:localhost/replica:0/task:0/cpu:0
2017-07-23 12:31:25.796368: I tensorflow/core/common_runtime/simple_placer.cc:847] cnn0/BatchNorm/Const: (Const)/job:localhost/replica:0/task:0/cpu:0
Initialization time 15.614956
```

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

### List TFRecord Files

Source: https://github.com/calico/basenji/blob/master/tutorials/preprocess.ipynb

Lists all TFRecord files within the specified directory. Useful for verifying data availability and organization.

```bash
! ls -l data/heart_l131k/tfrecords/*.tfr
```

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

### Load Data and Targets

Source: https://github.com/calico/basenji/blob/master/jupyter/basenji_sad_viz.ipynb

Opens the SAD HDF5 file and reads target information into a pandas DataFrame. Determines the number of targets.

```python
# open scores h5
sad_h5 = h5py.File(sad_h5_file, 'r')

# read targets
targets_df = pd.read_table(targets_file, index_col=0, usecols=['index', 'identifier', 'description', 'file'])
num_targets = targets_df.shape[0]
```

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

### Run Basenji Test Script

Source: https://github.com/calico/basenji/blob/master/tutorials/train_test.ipynb

Executes the basenji_test.py script with specified parameters to benchmark model accuracy. This includes options for averaging predictions and specifying model/data files.

```bash
! basenji_test.py --ai 0,1,2 -o output/heart_test --rc --shifts "1,0,-1" models/params_small.json models/heart/model_best.h5 data/heart_l131k
```

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

### Display First 3 Test Sequences

Source: https://github.com/calico/basenji/blob/master/tutorials/preprocess.ipynb

Filters the sequences BED file for 'test' entries and displays the first three. Useful for inspecting the 'test' split data.

```python
! grep test data/heart_l131k/sequences.bed | head -n3
```

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

### Import Python Libraries

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Imports essential Python libraries for file operations, JSON handling, and subprocess execution.

```python
import json
import os
import shutil
import subprocess
```

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

### Download and Unzip FASTA File

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Checks if the simplified human genome FASTA file (hg38.ml.fa) exists. If not, it downloads and unzips the gzipped version.

```python
if not os.path.isfile('./data/hg38.ml.fa'):
    print('downloading hg38.ml.fa')
    subprocess.call('curl -o ./data/hg38.ml.fa.gz https://storage.googleapis.com/basenji_barnyard2/hg38.ml.fa.gz', shell=True)
    subprocess.call('gunzip ./data/hg38.ml.fa.gz', shell=True)
```

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

### Import Libraries

Source: https://github.com/calico/basenji/blob/master/tutorials/sat_mut.ipynb

Imports necessary Python libraries for file operations, subprocess execution, and displaying interactive elements.

```python
import os
import subprocess
from IPython.display import IFrame
```

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

### Model Configuration Parameters

Source: https://github.com/calico/basenji/blob/master/tutorials/sad.ipynb

These parameters define the strides, lengths, and crops used for the model's target processing. They are essential for configuring the input data and model's receptive field.

```text
model_strides [128]
target_lengths [1024]
target_crops [0]
```

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

### Run Akita Data Processing Script

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/tutorial.ipynb

Executes the `akita_data.py` script with specified parameters to process genomic and Hi-C data, creating TFRecords for training, validation, and testing. This includes down-sampling, specifying gap regions, sequence length, output directory, and train/validation/test splits.

```bash
! akita_data.py --sample 0.05 -g ./data/hg38_gaps_binsize2048_numconseq10.bed -l 1048576 --crop 65536 --local -o ./data/1m --as_obsexp -p 8 -t .1 -v .1 -w 2048 --snap 2048 --stride_train 262144 --stride_test 32768 ./data/hg38.ml.fa ./data/microc_cools.txt
```

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

### Restore Model

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/explore_model.ipynb

Restores a pre-trained model from a file. Ensure the model file is downloaded before running.

```python
# note: run %%bash get_model.sh 
# if you have not already downloaded the model
seqnn_model.restore(model_file)
print('successfully loaded')
```

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

### Download Genome FASTA

Source: https://github.com/calico/basenji/blob/master/tutorials/sat_mut.ipynb

Downloads the human genome FASTA file (hg19.ml.fa) and its index if they do not already exist. This is required for sequence extraction.

```python
if not os.path.isfile('data/hg19.ml.fa'):
    subprocess.call('curl -o data/hg19.ml.fa https://storage.googleapis.com/basenji_tutorial_data/hg19.ml.fa', shell=True)
    subprocess.call('curl -o data/hg19.ml.fa.fai https://storage.googleapis.com/basenji_tutorial_data/hg19.ml.fa.fai', shell=True)
```

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

### Check for Near Equality with Different Tolerance

Source: https://github.com/calico/basenji/blob/master/tests/test_sad_debug.ipynb

Another check using np.allclose with different tolerance values. This helps in understanding the sensitivity of the comparison to tolerance settings.

```python
np.allclose(saved_sad, this_sad, atol=0.2, rtol=0.2)
```

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

### Define File Paths

Source: https://github.com/calico/basenji/blob/master/jupyter/basenji_sad_viz.ipynb

Specifies the paths to the SAD HDF5 file and the targets file.

```python
sad_h5_file = 'ponies_sad/sad.h5'
targets_file = '/home/drk/farm/data/basenji/data/v4/targets_h.txt'
```

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

### Import Libraries for SAD Debugging

Source: https://github.com/calico/basenji/blob/master/tests/test_sad_debug.ipynb

Imports necessary libraries including h5py for HDF5 file handling, numpy for numerical operations, sklearn for metrics, and matplotlib/seaborn for plotting.

```python
import h5py
import numpy as np
from sklearn.metrics import explained_variance_score

import matplotlib.pyplot as plt
import seaborn as sns
```

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

### Import Libraries

Source: https://github.com/calico/basenji/blob/master/jupyter/basenji_sat_vcf_viz.ipynb

Imports necessary libraries for data manipulation, plotting, and scientific computing.

```python
import h5py
import logomaker
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns
```

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

### Download CAGE BigWig Files

Source: https://github.com/calico/basenji/blob/master/tutorials/preprocess.ipynb

Downloads three CAGE BigWig files related to heart biology if they do not already exist. These files represent genomic coverage data.

```python
if not os.path.isfile('data/CNhs11760.bw'):
    subprocess.call('curl -o data/CNhs11760.bw https://storage.googleapis.com/basenji_tutorial_data/CNhs11760.bw', shell=True)
    subprocess.call('curl -o data/CNhs12843.bw https://storage.googleapis.com/basenji_tutorial_data/CNhs12843.bw', shell=True)
    subprocess.call('curl -o data/CNhs12856.bw https://storage.googleapis.com/basenji_tutorial_data/CNhs12856.bw', shell=True)
```

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

### Basenji Model Architecture Summary

Source: https://github.com/calico/basenji/blob/master/tutorials/train_test.ipynb

Displays a summary of the Basenji model architecture, including layer types, output shapes, and parameter counts. This is useful for understanding the model's structure.

```text
Model: "model_1"
__________________________________________________________________________________________________
 Layer (type)                Output Shape                 Param #   Connected to                  
==================================================================================================
 sequence (InputLayer)       [(None, 131072, 4)]          0         []                            
                                                                                                  
 stochastic_reverse_complem  ((None, 131072, 4),          0         ['sequence[0][0]']            
 ent (StochasticReverseComp   ())
 lement)
                                                                                                  
 stochastic_shift (Stochast  (None, 131072, 4)            0         ['stochastic_reverse_complemen
 icShift)
                                                                                                 
                                                                                                  
 tf.nn.gelu (TFOpLambda)     (None, 131072, 4)            0         ['stochastic_shift[0][0]']    
                                                                                                  
 conv1d (Conv1D)             (None, 131072, 64)           3840      ['tf.nn.gelu[0][0]']          
                                                                                                  
 
```

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

### View Gene Prediction and Measurement Table

Source: https://github.com/calico/basenji/blob/master/tutorials/archive/genes.ipynb

Unzip and display the first few lines of the gene prediction and experimental measurement table using 'gunzip -c' and 'head'. This file contains specific gene predictions and their corresponding experimental values.

```bash
! gunzip -c output/gencode_chr9_test/gene_table.txt.gz | head
```

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

### Import Libraries

Source: https://github.com/calico/basenji/blob/master/jupyter/basenji_sad_viz.ipynb

Imports necessary libraries for data manipulation, file handling, and plotting.

```python
import h5py
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns
```

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

### Prepare DNA Sequence for Prediction

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/explore_model.ipynb

Fetches a DNA sequence from a FASTA file for a given region, converts it to 1-hot encoding, and checks if its length matches the expected sequence length. Requires 'pysam' and 'dna_io'.

```python
# this example uses the sequence for the test set region
# with the corresponding test_index, but
# predictions can be made for any DNA sequence of length = seq_length = 2^20

chrm, seq_start, seq_end = sequences_test.iloc[test_index][0:3] 
seq = fasta_open.fetch( chrm, seq_start, seq_end ).upper()
if len(seq) != seq_length: raise ValueError('len(seq) != seq_length')

# seq_1hot is a np.array with shape [2^20 bp, 4 nucleotides]
# representing 1-hot encoded DNA sequence
seq_1hot = dna_io.dna_1hot(seq)
```

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

### TensorFlow CUDA Warnings

Source: https://github.com/calico/basenji/blob/master/manuscripts/akita/explore_model.ipynb

These are informational warnings from TensorFlow indicating that CUDA dynamic libraries could not be loaded and that the driver might not be running. The binary is optimized for specific CPU instructions.

```text
2022-06-10 15:02:50.723562: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcuda.so.1'; dlerror: libcuda.so.1: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /home/drk/code/torch_1-2018/install/lib:/usr/local/cuda-11.2/targets/x86_64-linux/lib:/home/drk/code/torch_1-2018/install/lib:/usr/local/cuda-11.2/targets/x86_64-linux/lib:
2022-06-10 15:02:50.723690: W tensorflow/stream_executor/cuda/cuda_driver.cc:269] failed call to cuInit: UNKNOWN ERROR (303)
2022-06-10 15:02:50.723735: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:156] kernel driver does not appear to be running on this host (cb2-drk): /proc/driver/nvidia/version does not exist
2022-06-10 15:02:50.724429: I tensorflow/core/platform/cpu_feature_guard.cc:151] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.
```

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

### Model Training Progress

Source: https://github.com/calico/basenji/blob/master/tutorials/train_test.ipynb

Displays the training progress of the Basenji model, showing loss and evaluation metrics per epoch. This output is typical during the training phase.

```text
45/45 [==============================] - 10s 184ms/step - loss: 0.3362 - pearsonr: 0.5476 - r2: 0.2953
```