<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Preparing a dataset for Basenji training involves a series of design choices.
The input you bring to the pipeline is
Step1:  Next, let's grab a few CAGE datasets from FANTOM5 related to heart biology.
These data were processed by
1. Aligning with Bowtie2 with very sensitive alignment parameters.
2. Distributing multi-mapping reads and estimating genomic coverage with bam_cov.py
Step2:  Then we'll write out these BigWig files and labels to a samples table.
Step3:  Next, we want to choose genomic sequences to form batches for stochastic gradient descent, divide them into training/validation/test sets, and construct TFRecords to provide to downstream programs.
The script basenji_data.py implements this procedure.
The most relevant options here are
Step4:  Now, data/heart_l131k contains relevant data for training.
contigs.bed contains the original large contiguous regions from which training sequences were taken (possibly strided).
sequences.bed
contains the train/valid/test sequences.
<ASSISTANT_TASK:>
Python Code:
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)                
Explanation: Preparing a dataset for Basenji training involves a series of design choices.
The input you bring to the pipeline is:
* BigWig coverage tracks
* Genome FASTA file
First, make sure you have an hg19 FASTA file visible. If you have it already, put a symbolic link into the data directory. Otherwise, I have a machine learning friendly simplified version you can download in the next cell.
End of explanation
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)
Explanation: Next, let's grab a few CAGE datasets from FANTOM5 related to heart biology.
These data were processed by
1. Aligning with Bowtie2 with very sensitive alignment parameters.
2. Distributing multi-mapping reads and estimating genomic coverage with bam_cov.py
End of explanation
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()
Explanation: Then we'll write out these BigWig files and labels to a samples table.
End of explanation
! basenji_data.py -d .1 -g data/unmap_macro.bed -l 131072 --local -o data/heart_l131k -p 8 -t .1 -v .1 -w 128 data/hg19.ml.fa data/heart_wigs.txt
Explanation: Next, we want to choose genomic sequences to form batches for stochastic gradient descent, divide them into training/validation/test sets, and construct TFRecords to provide to downstream programs.
The script basenji_data.py implements this procedure.
The most relevant options here are:
| Option/Argument | Value | Note |
|:---|:---|:---|
| -d | 0.1 | Down-sample the genome to 10% to speed things up here. |
| -g | data/unmap_macro.bed | Dodge large-scale unmappable regions like assembly gaps. |
| -l | 131072 | Sequence length. |
| --local | True | Run locally, as opposed to on my SLURM scheduler. |
| -o | data/heart_l131k | Output directory |
| -p | 8 | Uses multiple concourrent processes to read/write. |
| -t | .1 | Hold out 10% sequences for testing. |
| -v | .1 | Hold out 10% sequences for validation. |
| -w | 128 | Pool the nucleotide-resolution values to 128 bp bins. |
| fasta_file| data/hg19.ml.fa | FASTA file to extract sequences from. |
| targets_file | data/heart_wigs.txt | Target samples table with BigWig paths. |
End of explanation
! cut -f4 data/heart_l131k/sequences.bed | sort | uniq -c
! head -n3 data/heart_l131k/sequences.bed
! grep valid data/heart_l131k/sequences.bed | head -n3
! grep test data/heart_l131k/sequences.bed | head -n3
! ls -l data/heart_l131k/tfrecords/*.tfr
Explanation: Now, data/heart_l131k contains relevant data for training.
contigs.bed contains the original large contiguous regions from which training sequences were taken (possibly strided).
sequences.bed
contains the train/valid/test sequences.
End of explanation
<END_TASK>