seqdataloader

Sequence data label generation and ingestion into deep learning models

Installation

pip install seqdataloader

If you plan to modify the code, you can install it in development mode: pip install -e seqdataloader

Please note: to use the "dbingest" functionality in seqdataloader, python>=3.7 is needed.

Quick Start

labelgen

The input for the labelgen submodule is a 4 column tab-delimited file with the following fields:

• "task" -- required. User-specified task name
• "narrowPeak" -- Path to narrowPeak file. (Optional if "bigwig" is specified.)
• "bigwig" -- Path to bigwig file (optional if "narrowPeak" is specified.)
• "ambig" -- bed file containing user-specified regions to label as ambiguous (optional)

genomewide_labels --task_list tasks.tsv \
		  --outf classificationlabels.SummitWithin200bpCenter.tsv.gz \
		  --output_type gzip \ # (one of gzip, bz2, hdf5, pkl) 
		  --chrom_sizes hg38.chrom.sizes \
		  --bin_stride 50 \
		  --left_flank 400 \
		  --right_flank 400 \
		  --bin_size 200 \
		  --task_hreads 10 \
		  --chrom_threads 4 \
		  --allow_ambiguous \
		  --labeling_approach peak_summit_in_bin_classification 

And for regression:

genomewide_labels --task_list tasks.tsv \
       --outf regressionlabels.allbins.hg38.hdf5 \
       --output_type hdf5 \
       --chrom_sizes hg38.chrom.sizes \
       --bin_stride 50 \
       --left_flank 400 \
       --right_flank 400 \
       --chrom_threads 24 \
       --task_threads 2 \
       --label_transformer asinh \ one of None, asinh, log10, log, default is asinh 
       --labeling_approach all_genome_bins_regression

labeling_approach can be one of:

"peak_summit_in_bin_classification"

"peak_percent_overlap_with_bin_classification"

"peak_summit_in_bin_regression"

"peak_percent_overlap_with_bin_regression"

"all_genome_bins_regression"

How to run

Sample datasets are included in the folder examples/peak_files_from_encode_for_label_comparison and examples/bigwig_files_from_encode_for_label_comparison

Executing seqdataloader as a script:

Execute the script:

examples/genomewide_labels.sh for examples on how to generate classification and regression labels on sample datasets. The script generates binary classification labels (1,0,-1 for ambiguous) or continuous regression labels reflective of bigWig coverage in a bin in bed file format:

http://mitra.stanford.edu/kundaje/seqdataloader/classificationlabels.50PercentOverlap.tsv.gz

http://mitra.stanford.edu/kundaje/seqdataloader/classificationlabels.SummitWithin200bpCenter.tsv.gz

http://mitra.stanford.edu/kundaje/seqdataloader/regressionlabels.50PercentOverlap.tsv.gz

http://mitra.stanford.edu/kundaje/seqdataloader/regressionlabels.SummitWithin200bpCenter.tsv.gz

Corresponding WashU Browser Tracks with optimal narrowPeak and associated bin labels are here: http://epigenomegateway.wustl.edu/legacy/?genome=hg38&session=GDB2BTMGnB&statusId=1154897038

calling seqdataloader as a Python function:

from seqdataloader import *
classification_params={
    'task_list':"tasks.tsv",
    'outf':"classificationlabels.SummitWithin200bpCenter.tsv.gz",
    'output_type':'gzip',
    'chrom_sizes':'hg38.chrom.sizes',
    'chroms_to_keep':['chr21'],
    "store_positives_only":True,
    'bin_stride':50,
    'left_flank':400,
    'right_flank':400,
    'bin_size':200,
    'chrom_threads':10,
    'task_threads':4,
    'allow_ambiguous':True,
    'labeling_approach':'peak_summit_in_bin_classification'
    }
genomewide_labels(classification_params)

regression_params={
    'task_list':"tasks.tsv",
    'outf':"regressionlabels.all_genome_bins_regression.hdf5",
    'output_type':'hdf5',
    'chrom_sizes':'hg38.chrom.sizes',
    'store_values_above_thresh': 0,
    'chroms_to_keep':['chr21'],
    'bin_stride':50,
    'left_flank':400,
    'right_flank':400,
    'bin_size':200,
    'chrom_threads':10,
    'task_threads':4,
    'labeling_approach':'all_genome_bins_regression',
    'label_transformer':'log10',
    'label_transfomer_pseudocount':0.001
    }
genomewide_labels(regression_params)

Regression label transformations

In regression mode ( "peak_summit_in_bin_regression", "peak_percent_overlap_with_bin_regression", "all_genome_bins_regression"), the generated labels can be transformed in one of several ways. You can use the arguments label_transformer and label_transformer_pseudocount to specify the desired tranformation. Allowed values are:

• asinh -- numpy.arcsinh(values) will be computed (this is the default)
• None -- no label transformation will be performed
• log10 -- numpy.log10(values + pseudocount) will be computed using a pseudocount specified by label_transformer_pseudocount argument. If this argument is not provided,a default pseudocout of 0.001 is used.
• log -- numpy.log(values + pseudocount) will be computed using a pseudcount as above.

A note on file outputs

The code supports several output types: hdf5, gzip, pkl, bz2. Specify your desired output type with the flag --output_type. The default setting for this flag is gzip Please note that the large bottleneck in the code is writing the files to disk. hdf5 has negligible overhead, but using gzip or bz2 may increase runtime. Timining benchmarks are provided in examples/genomewide_labels.sh

You may speed up i/o by writing chromosome outputs to separate files in parallel. This is currently only supported for the gzip and bz2 output types, as i/o is less of a bottleneck for hdf5 and pkl output formats. Use the flag --split_output_by_chrom to invoke this parallelized saving of chromosomes.

dbingest

The input tsv file must have a subset of the columns corresponding to the supported configurations:

 * encode_config 
        ** dataset
	** fc_bigwig 
	** pval_bigwig
	** count_bigwig_plus_5p
	** count_bigwig_minus_5p
	** count_bigwig_unstranded_5p
	** idr_peak
	** overlap_peak 
	** ambig_peak 
	
* generic_bigwig 
	** bigwig_track 
	

Dependencies

Please make sure the following dependencies are installed on your system to use SeqDataLoader:

• pybedtools
• pyBigWig
• pandas
• numpy
• multiprocessing

Documentation and benchmarks

Testing, benchmarks, and documentation can be found in the docs folder