A Python package for fast operations on 1-dimensional genomic signal tracks.
- Finds the mean, approx. mean, max, min, coverage, or standard deviation for a given interval in a bedGraph file
- Partly written in Cython for speed improvements
- Can look up exact statistics of 1 million regions in ~0.26 second on a conventional laptop
- An approximate mean for 10,000 regions can be computed in ~0.0012 second w/ error rate of less than 5 percent
- Uses memory to load files
- 16 bytes per line in bedGraph file
- 4 bytes per basePair in every chromosome loaded
- Loading the bedGraph file can take up to a minute or two
- Only works with sorted bedgraph files
Dependency requirements:
- Numpy >= v1.16.4
- pyBigWig >= v0.3.16 (For reading bigWig files)
- pyBigWig == 0.3.16 (For Benchmarking)
With pip:
pip3 install pyBedGraph
pip3 install pyBigWig # if using bigwig filesWith conda:
conda create -n test
conda activate test
conda install -c bioconda pyBedGraph
conda install -c bioconda pyBigWig # if using bigwig fileshttps://thejacksonlaboratory.ent.box.com/s/3jglutwf3d54pnomnp33ivo7a9546vhe
Test files are also included in this Github repository: test/test_files.
Enter the directory with the test files.
from pyBedGraph import BedGraph
# arg1 - chromosome sizes file
# arg2 - bedgraph file
# arg3 - (optional) chromosome_name
# Just load chromosome 'chr1' (uses less memory and takes less time)
bedGraph = BedGraph('myChrom.sizes', 'random_test.bedGraph', 'chr1')
# Load the whole bedGraph file
bedGraph = BedGraph('myChrom.sizes', 'random_test.bedGraph')
# Option to not ignore missing basePairs when calculating statistics
# Used the exact same way but produces slightly different results
inclusive_bedGraph = BedGraph('myChrom.sizes', 'random_test.bedGraph', ignore_missing_bp=False)bedGraph.load_chrom_data('chr1')
inclusive_bedGraph.load_chrom_data('chr1')For approx_mean:
- Smaller bin size -> more accurate but slower
- Larger bin size -> less accurate but faster
bedGraph.load_chrom_bins('chr1', 3)
inclusive_bedGraph.load_chrom_bins('chr1', 3)'mean''approx_mean'- an approximate mean is faster than exact mean, with < 5% error rate'max''min''coverage''std'- (population standard deviation)'sum'
import numpy as np
# Option 1
test_intervals = [
['chr1', 24, 26],
['chr1', 12, 15],
['chr1', 8, 12],
['chr1', 9, 10],
['chr1', 0, 5]
]
values = bedGraph.stats(intervals=test_intervals)
# [-1. 0.9 0.1 -1. 0.82 0.72222222]
print(values)
# Option 2
start_list = np.array([24, 12, 8, 9, 0], dtype=np.int32)
end_list = np.array([26, 15, 12, 10, 5], dtype=np.int32)
chrom_name = 'chr1'
# arg1 - (optional) stat (default is 'mean')
# arg2 - intervals
# arg3 - start_list
# arg4 - end_list
# arg5 - chrom_name
# must have either intervals or start_list, end_list, chrom_name
# returns a numpy array of values
result = bedGraph.stats(start_list=start_list, end_list=end_list, chrom_name=chrom_name)
# [-1. 0.9 0.1 -1. 0.82 0.72222222]
print(result)# arg1 - interval file
# arg2 - (optional) output_to_file (default is True and outputs to 'chr1_out.txt'
# arg3 - (optional) stat (default is 'mean')
# returns a dictionary; keys are chromosome names, values are numpy arrays
result = bedGraph.stats_from_file('test_intervals.txt', output_to_file=False, stat='mean')
# {'chr1': array([-1. , 0.9 , 0.1 , -1. , 0.82 ,
# 0.72222222])}
print(result)# [-1. 0.9 0.1 -1. 0.82 0.72222222]
bedGraph.stats('mean', test_intervals)
# [-1. 0.9 0.1 -1. 0.80769231 0.72222222]
bedGraph.stats('approx_mean', test_intervals)
# [0. 0.33333333 0.25 0. 1. 0.3 ]
bedGraph.stats('coverage', test_intervals)
# [-1. 0.9 0.1 -1. 0.7 0.1]
bedGraph.stats('min', test_intervals)
# [-1. 0.9 0.1 -1. 0.9 1. ]
bedGraph.stats('max', test_intervals)
# [-1. 0. 0. -1. 0.09797959 0.27799991]
bedGraph.stats('std', test_intervals)# [0. 0.3 0.025 0. 0.82 0.21666667]
inclusive_bedGraph.stats('mean', test_intervals)
# [0. 0.3 0.00833333 0. 0.7 0.21666667]
inclusive_bedGraph.stats('approx_mean', test_intervals)
# [0. 0.33333333 0.25 0. 1. 0.3 ]
inclusive_bedGraph.stats('coverage', test_intervals)
# [0. 0. 0.1 0. 0.7 0.1]
inclusive_bedGraph.stats('min', test_intervals)
# [0. 0.9 0.1 0. 0.9 1. ]
inclusive_bedGraph.stats('max', test_intervals)
# [0. 0.42426407 0.04330127 0. 0.09797959 0.36431061]
inclusive_bedGraph.stats('std', test_intervals)Actual values are found from the stats function in pyBigWig with the exact argument being True. The error for exact stats will be ~1e-8 due to rounding error of conversion of bigWig and bedGraph files.
Alternatively, one can make actual values be pyBedGraph's exact statistics.
Enter the graphs folder in the Github project repository.
from pyBedGraph import BedGraph
from Benchmark import Benchmark
# These files can be downloaded from the link given above
bedGraph = BedGraph('mm10.chrom.sizes', 'ENCFF376VCU.bedGraph', 'chr1')
# Alternatively using a bigwig file
# bedGraph = BedGraph('mm10.chrom.sizes', 'ENCFF376VCU.bigWig', 'chr1')
bedGraph.load_chrom_data('chr1')
bedGraph.load_chrom_bins('chr1', 100)
# arg1 - BedGraph object
# arg2 - bigwig file
bench = Benchmark(bedGraph, 'ENCFF376VCU.bigWig')
# arg1 - num_tests
# arg2 - interval_size
# arg3 - chrom_nam
# arg4 - bin_size
# arg5 - stats (optional) (Default is all stats)
# arg6 - just_runtime (optional) (Default is False)
# arg6 - bench_pyBigWig_approx (optional) (Default is True)
# arg6 - make_pyBigWig_baseline (optional) (Default is True)
# Test all statistics
result = bench.benchmark(10000, 5000, 'chr1', 100)
for key in result:
print(key, result[key])
# mean {'run_time': 0.008324861526489258, 'error': {'percent_error': 0.0, 'ms_error': 0.0, 'abs_error': 0.0, 'not_included': 0}}
# pyBigWig_mean {'approx_run_time': 1.4333949089050293, 'exact_run_time': 0.7698564529418945, 'error': {'percent_error': 0.06567272540694802, 'ms_error': 0.001222419386871348, 'abs_error': 0.023540340949669364, 'not_included': 79}}
# approx_mean {'run_time': 0.002111673355102539, 'error': {'percent_error': 0.006529644707171326, 'ms_error': 7.858080037556034e-06, 'abs_error': 0.001824641073039555, 'not_included': 4}}
# max {'run_time': 0.005040645599365234, 'error': {'percent_error': 0.0, 'ms_error': 0.0, 'abs_error': 0.0, 'not_included': 0}}
# pyBigWig_max {'approx_run_time': 1.2673799991607666, 'exact_run_time': 0.7933700084686279, 'error': {'percent_error': 0.10220448242023446, 'ms_error': 1.2678718593032368, 'abs_error': 0.25865022624731066, 'not_included': 79}}
# min {'run_time': 0.005083560943603516, 'error': {'percent_error': 0.0, 'ms_error': 0.0, 'abs_error': 0.0, 'not_included': 0}}
# pyBigWig_min {'approx_run_time': 1.2120039463043213, 'exact_run_time': 0.7468140125274658, 'error': {'percent_error': 0.0001, 'ms_error': 7.109862619931795e-07, 'abs_error': 8.432000130414962e-06, 'not_included': 0}}
# coverage {'run_time': 0.0063626766204833984, 'error': {'percent_error': 0.0, 'ms_error': 0.0, 'abs_error': 0.0, 'not_included': 0}}
# pyBigWig_coverage {'approx_run_time': 1.2101118564605713, 'exact_run_time': 0.7483360767364502, 'error': {'percent_error': 0.0, 'ms_error': 0.0, 'abs_error': 0.0, 'not_included': 0}}
# std {'run_time': 0.0422673225402832, 'error': {'percent_error': 9.690484548456011e-05, 'ms_error': 4.764358150024449e-09, 'abs_error': 6.25265457158463e-05, 'not_included': 0}}
# pyBigWig_std {'approx_run_time': 1.219078540802002, 'exact_run_time': 0.7484426498413086, 'error': {'percent_error': 0.04560011737269686, 'ms_error': 0.005008324729263816, 'abs_error': 0.02569405301725115, 'not_included': 79}}Some tests are provided in test/test.py. Additional bedgraph and bigwig files for ENCFF376VCU are needed to run extensive_test.py. Build badge comes from a forked repository, https://github.com/c0ver/pyBedGraph, that has the same version as this repository.
pyBedGraph: a Python package for fast operations on 1-dimensional genomic signal tracks, Zhang et al., Bioinformatics, 2020
pyBedGraph/include_missing_bp.pyx in pyBedGraph.include_missing_bp.get_exact_means()
pyBedGraph/include_missing_bp.pyx in pyBedGraph.include_missing_bp.get_exact_means()
IndexError: Out of bounds on buffer access (axis 0)
This is frequently caused by giving an interval that is outside the chromosome size.
Note that bedgraph files need to be sorted. To do so, try
sort -k1,1 -k2,2 example.bedgraph > example.sorted.bedgraph
To report bugs, contact Henry (henrybzhang.99@gmail.com) and Minji (minji.kim@jax.org) or visit the Issues page.