Sequence data label generation and ingestion into deep learning models
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.
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"
Sample datasets are included in the folder examples/peak_files_from_encode_for_label_comparison
and examples/bigwig_files_from_encode_for_label_comparison
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
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)
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.
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.
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
Please make sure the following dependencies are installed on your system to use SeqDataLoader:
- pybedtools
- pyBigWig
- pandas
- numpy
- multiprocessing
Testing, benchmarks, and documentation can be found in the docs
folder