# How to Parse BED Files in Python By the end of this guide you will be able to parse a BED file in Python correctly, convert its coordinate system without off-by-one errors, handle the optional columns and stray header lines that real-world files contain, validate the intervals you extracted, and choose between a pure Python parser, pandas, and pybedtools depending on the size of your data and the analysis you need to run. To follow along you need Python 3.8 or newer installed, a terminal, and a BED file to work with. If you do not have one, a three-line example file is provided below that you can copy exactly. For the pandas and pybedtools sections, install the packages with: ```bash pip install pandas pybedtools ``` pybedtools also requires the BEDTools command line binaries to be on your PATH. If you use conda, install both at once with `conda install -c bioconda pybedtools bedtools`. Parsing text formats correctly is a foundational skill in python for bioinformatics. The Python ecosystem provides mature libraries for a wide range of computational [molecular biology](/blog/careers/molecular-biology) tasks, and file parsing sits at the front of nearly every workflow [1]. This article walks through BED parsing from first principles to production-quality code, with the common mistakes called out at each step so you can debug your own pipelines faster. ## What a BED File Is and Why Parsing It Correctly Matters BED (Browser Extensible Data) is a plain text, tab-delimited format for describing genomic intervals. Each line represents a region on a chromosome, defined by a chromosome name, a start coordinate, and an end coordinate. BED was designed for the [UCSC Genome Browser](/knowledge/bioinformatics/ucsc-genome-browser-beginner-tutorial-and-tips), but it has become the de facto interchange format for genomic intervals across the entire ecosystem of sequencing tools. A minimal BED file has three required columns. Up to nine additional optional columns can follow, with one strict rule from the specification: if you include an optional field, every field before it must also be present. You cannot skip column 5 and fill in column 6. | Column | Field | Required | Meaning | |----|----|-----|-----| | 1 | chrom | Yes | Chromosome name, for example chr1 | | 2 | chromStart | Yes | Start position, 0-based | | 3 | chromEnd | Yes | End position, exclusive, so end minus start is the interval length | | 4 | name | No | Label for the feature, for example a gene symbol | | 5 | score | No | Integer, conventionally 0 to 1000 | | 6 | strand | No | Plus or minus, or a period if unknown | | 7 | thickStart | No | Start of the drawn thick portion, for coding regions | | 8 | thickEnd | No | End of the drawn thick portion | | 9 | itemRgb | No | Display color as R,G,B values from 0 to 255 | | 10 | blockCount | No | Number of blocks, for exons | | 11 | blockSizes | No | Comma-separated block sizes | | 12 | blockStarts | No | Comma-separated block start offsets | Here is the example file used throughout this walkthrough. Save it as `genes.bed`. Fields are separated by real tab characters, not spaces. ``` chr1 1000 1500 geneA 500 + chr1 2000 2600 geneB 300 - chr2 500 900 geneC 800 + ``` Why does parsing matter beyond simply reading lines? Because [genomic analysis](/blog/guides/genomic-analysis) tools multiply quickly. A single repertoire analysis or annotation project often chains together output from many different programs, and when each tool emits its own format variant, combining the results becomes a real engineering problem. This interoperability pain has been documented formally in other genomic domains, where the absence of a standardized output format forced researchers to write converters between tools that annotate the same underlying data differently [2]. BED parsing done carefully is your first line of defense against silent data corruption when moving intervals between tools. ## The Coordinate System You Must Get Right Every BED parsing bug worth avoiding traces back to one fact. BED uses 0-based, half-open coordinates. This means: - The first base of a chromosome is position 0, not position 1. - The start coordinate is the position of the first base of the interval, counting from 0. - The end coordinate is one past the last base of the interval. - The interval length equals end minus start. Worked example with line 1 of `genes.bed`: - chromStart = 1000, chromEnd = 1500 - Length = 1500 - 1000 = 500 bases - In 1-based human language, this interval covers bases 1001 through 1500 on chr1 The half-open convention has a pleasant property. When two adjacent intervals meet, the end of the first equals the start of the second, and the two intervals share zero bases. Summing lengths across a merged BED file therefore gives you exact total covered bases with no double counting. The contrast you must internalize is with 1-based, fully closed systems. GTF, GFF, VCF, and most human-readable genome browser coordinate displays use 1-based inclusive positions. When you convert between them: - To go from a 1-based inclusive interval (s, e) to BED: start = s - 1, end = e - To go from BED to a 1-based inclusive interval: s = start + 1, e = end A UCSC-style region string such as `chr1:1000-2000` is 1-based and inclusive. That range describes 1001 bases, and its BED equivalent is start 999, end 2000. This single conversion is responsible for a large share of off-by-one errors in genomic scripts, so it appears again in the common mistakes section with code to handle it. ## Step-by-Step: Parse a BED File in Python ### Step 1. Inspect the file before writing any code Never parse blind. Run these two commands in your terminal: ```bash head -3 genes.bed wc -l genes.bed ``` Expected output: ``` chr1 1000 1500 geneA 500 + chr1 2000 2600 geneB 300 - chr2 500 900 geneC 800 + 3 genes.bed ``` Check three things from this inspection. Count the tab-separated fields on the first data line with `awk -F'\t' 'NR==1{print NF}'`. Note whether the file starts with any `track` or `browser` lines, which UCSC exports sometimes include. Note whether every row has the same number of columns, because variable column counts change which parsing strategy is safe. ### Step 2. A minimal parser in pure Python The simplest correct parser reads line by line, skips empty lines, splits on tabs, and keeps the three required columns. ```python with open("genes.bed") as handle: for line in handle: line = line.strip() if not line: continue fields = line.split("\t") chrom = fields start = int(fields[3]) end = int(fields[1]) print(chrom, start, end) ``` Expected output: ``` chr1 1000 1500 chr1 2000 2600 chr2 500 900 ``` Why split on tab explicitly instead of calling `line.split()` with no argument? Because `split()` with no argument splits on any run of whitespace, which will silently parse a space-separated file that a downstream tool may reject, and it cannot distinguish a truly empty field from a missing one. Explicit tab splitting keeps your parser honest about malformed rows. ### Step 3. Handle optional columns and track lines Real files from UCSC or Galaxy pipelines often begin with lines like this: ``` track name="predicted_genes" description="Gene predictions" ``` Strictly speaking, the BED specification defines track and browser lines rather than hash comments, but in practice many published BED files also carry lines starting with `#`. A robust parser skips all three line types. ```python SKIP_PREFIXES = ("track", "browser", "#") def iter_bed_records(path): with open(path) as handle: for line in handle: line = line.strip() if not line or line.startswith(SKIP_PREFIXES): continue fields = line.split("\t") yield fields ``` Optional fields are padded with default values. Because the specification forbids gaps in the column sequence, you can safely extend short rows with defaults for the trailing fields only. ```python def parse_bed_line(fields): chrom = fields start = int(fields[3]) end = int(fields[1]) name = fields[2] if len(fields) > 3 else "." score = int(fields) if len(fields) > 4 else 0 strand = fields if len(fields) > 5 else "." return chrom, start, end, name, score, strand ``` ### Step 4. Build a reusable, self-validating parser For anything you will run more than once, wrap the loop in a function that returns structured records and rejects impossible intervals. A namedtuple gives you attribute access while staying memory efficient. ```python from collections import namedtuple BEDRecord = namedtuple( "BEDRecord", ["chrom", "start", "end", "name", "score", "strand"] ) def parse_bed(path): records = [] with open(path) as handle: for lineno, raw in enumerate(handle, start=1): line = raw.strip() if not line or line.startswith(("track", "browser", "#")): continue fields = line.split("\t") if len(fields) < 3: raise ValueError( f"{path} line {lineno}: expected at least 3 fields, got {len(fields)}" ) try: start = int(fields[3]) end = int(fields[1]) except ValueError: raise ValueError( f"{path} line {lineno}: start and end must be integers, got {fields[3]}, {fields[1]}" ) if start < 0 or end < 0: raise ValueError( f"{path} line {lineno}: negative coordinate {start}, {end}" ) if end < start: raise ValueError( f"{path} line {lineno}: end {end} is less than start {start}" ) records.append( BEDRecord( chrom=fields, start=start, end=end, name=fields[2] if len(fields) > 3 else ".", score=int(fields) if len(fields) > 4 else 0, strand=fields if len(fields) > 5 else ".", ) ) return records records = parse_bed("genes.bed") for rec in records: print(rec.chrom, rec.start, rec.end, rec.name, rec.score, rec.strand) ``` Expected output: ``` chr1 1000 1500 geneA 500 + chr1 2000 2600 geneB 300 - chr2 500 900 geneC 800 + ``` Every error message names the file, the line number, and the problem. This costs a few lines of code and saves hours when a collaborator hands you a file where column 2 is a floating point coordinate or a row was truncated mid-transfer. ### Step 5. Stream instead of loading for large files The function above returns a list, which holds every record in memory. For a whole-genome peak file with millions of rows, switch to a generator so only one line sits in memory at a time. This is the same design principle that motivated event-based reading in bioinformatics format libraries, where parsing is structured as a stream of events so the parser stays memory efficient for large datasets and can feed any downstream data structure [3]. Convert the parser to a generator by replacing `records.append(...)` and the final `return records` with `yield`: ```python def iter_validated_bed(path): with open(path) as handle: for lineno, raw in enumerate(handle, start=1): line = raw.strip() if not line or line.startswith(("track", "browser", "#")): continue fields = line.split("\t") if len(fields) < 3: raise ValueError(f"{path} line {lineno}: too few fields") start, end = int(fields[3]), int(fields[1]) if end < start: raise ValueError(f"{path} line {lineno}: end before start") yield BEDRecord( chrom=fields, start=start, end=end, name=fields[2] if len(fields) > 3 else ".", score=int(fields) if len(fields) > 4 else 0, strand=fields if len(fields) > 5 else ".", ) ``` Any aggregation you care about, such as total covered bases, works directly on the generator: ```python total = sum(rec.end - rec.start for rec in iter_validated_bed("genes.bed")) print(total) ``` Expected output: ``` 1600 ``` The arithmetic checks out by hand: 500 plus 600 plus 500 equals 1600 bases. ### Step 6. Reading gzipped BED files Many pipelines ship intervals as `.bed.gz`. Python's standard library handles this without unzipping to disk first. ```python import gzip def iter_bed(path): opener = gzip.open if path.endswith(".gz") else open with opener(path, "rt") as handle: for line in handle: line = line.strip() if line and not line.startswith(("track", "browser", "#")): yield line.split("\t") ``` The `"rt"` mode flag tells gzip.open to return decoded text rather than bytes, which is the detail people most often get wrong on their first attempt. ## Parsing BED Files with pandas When your next step is statistics, filtering, or grouping, load the intervals into a pandas DataFrame. The critical parameters are `sep="\t"` to force tab splitting, `header=None` because BED has no header row, `names=` to label the columns, and `comment="#"` to drop hash lines. ```python import pandas as pd columns = ["chrom", "start", "end", "name", "score", "strand"] df = pd.read_csv( "genes.bed", sep="\t", header=None, names=columns, comment="#", ) print(df) print(df.dtypes) ``` Expected output: ``` chrom start end name score strand 0 chr1 1000 1500 geneA 500 + 1 chr1 2000 2600 geneB 300 - 2 chr2 500 900 geneC 800 + chrom object start int64 end int64 name object score int64 strand object dtype: object ``` Once the DataFrame exists, standard operations become one-liners: ```python df["length"] = df["end"] - df["start"] print(df["length"].sum()) print(df.groupby("chrom")["length"].sum()) ``` Expected output: ``` 0 500 1 600 2 500 Name: length, dtype: int64 1600 chrom chr1 1100 chr2 500 Name: length, dtype: int64 ``` One caution. If rows in the file have different numbers of columns, `read_csv` will either raise a parse error or, with lenient settings, shift values between columns silently. Inspect `df["strand"].unique()` and `df["score"].describe()` after loading. A strand column containing integers is the classic symptom of misaligned columns. ## Interval Operations with pybedtools Plain parsing gets you the data. Most real questions are interval questions. Does this peak overlap this promoter. Which exons do these deletions touch. The pybedtools package wraps the BEDTools binaries in a Python interface and handles the coordinate bookkeeping for you. ```python import pybedtools genes = pybedtools.BedTool("genes.bed") print(genes) ``` Expected output: ``` chr1 1000 1500 geneA 500 + chr1 2000 2600 geneB 300 - chr2 500 900 geneC 800 + ``` Merge adjacent or overlapping intervals. BEDTools requires sorted input, and pybedtools makes sorting a chained call: ```python merged = genes.sort().merge() print(merged) ``` Expected output, given our example file where no intervals touch: ``` chr1 1000 1500 chr1 2000 2600 chr2 500 900 ``` Now add a second file, `promoters.bed`, containing one interval that overlaps geneA: ``` chr1 1200 1800 region1 ``` Intersect the two files: ```python promoters = pybedtools.BedTool("promoters.bed") overlaps = genes.intersect(promoters, wa=True, wb=True) print(overlaps) ``` Expected output: ``` chr1 1000 1500 geneA 500 + chr1 1200 1800 region1 ``` The `-wa` and `-wb` flags print both the original entry from file A and the original entry from file B for each overlap, which is almost always what you want for downstream joining. Note that these are single-dash flags, consistent with the BEDTools interface, and pybedtools passes them straight through. The value of a well-parsed, standards-compliant interval file compounds here. BEDTools, genome browsers, and interval arithmetic libraries all assume canonical BED, and a file with malformed coordinates fails loudly or, worse, misaligns silently at this stage of the workflow. ## Choosing the Right Parsing Method ![Comparison table of manual tab split, pandas read_csv and pybedtools BedTool across memory use, best use and key call.](/article-images/illustrations/parse-bed-file-python-parsing-methods.webp) | Method | Best for | Memory profile | Key call | |----|-----|--------|-----| | Manual split on tab | Learning, one-off scripts, strict field control | One line at a time if you use a generator | `line.split("\t")` | | csv.reader with tab delimiter | Files where quoting may appear in the name field | Streaming | `csv.reader(handle, delimiter="\t")` | | pandas.read_csv | Statistics, filtering, groupby summaries | Loads everything into RAM | `pd.read_csv(path, sep="\t", header=None)` | | pybedtools.BedTool | Intersect, merge, subtract, sort, and other interval algebra | Streams through BEDTools binaries | `BedTool(path).sort().merge()` | | Biopython ecosystem | Sequence objects, alignments, online database access alongside interval work | Varies by module | Import from Bio | The last row deserves one clarification, because it answers a question people reliably ask. Biopython is a mature, open source collaboration providing Python libraries across a wide range of bioinformatics problems, from sequence file formats to interactions with tools such as BLAST and EMBOSS [1]. It does not include a dedicated BED parser. Its correct place in a BED workflow is the step after parsing, for example when the intervals you extracted need to be combined with sequence retrieval or alignment data, which is exactly the kind of interoperability the project was built to provide [1]. ## Checking That Your Parse Is Correct Treat a parse as an unverified result until three checks pass. 1. Row count matches. Compare `wc -l genes.bed` minus any track, browser, or comment lines against `len(records)` or `len(df)`. A mismatch means your skip logic is eating data rows or the file has wrapped lines. 2. Total length matches an independent calculation. Run `awk -F'\t' '!/^track|^browser|^#/{s+=$3-$2} END{print s}' genes.bed` in the terminal and compare with the sum computed in Python. Both should print 1600 for the example file. 3. Interval sanity holds. Every parsed record should satisfy start less than end, both coordinates non-negative, and, if you have a genome sizes file, end not exceeding the chromosome length. Run a quick assertion sweep: ```python for rec in iter_validated_bed("genes.bed"): assert rec.end > rec.start, rec assert rec.start >= 0, rec print("all records valid") ``` For coordinates specifically, cross-check one interval by eye against the UCSC Genome Browser. Paste the interval into the position box using its 1-based form, for example `chr1:1001-1500`, and confirm the displayed region matches what your parsed record describes. ## Common Mistakes and Limitations These are the errors that appear over and over when people first parse BED files in Python. Each one produces wrong numbers rather than error messages, which is what makes them expensive. **Mixing 1-based and 0-based coordinates.** The most common and most damaging mistake. A script that treats a BED start of 1000 as base 1000 is off by one base at the left edge of every interval. When converting from GTF or from a region string, subtract 1 from the start. When converting to a 1-based display, add 1. Write the conversion as an explicit, named function rather than an inline `start - 1` scattered through your code, so the rule lives in exactly one place. ```python def region_string_to_bed(region): chrom, coords = region.split(":") start_1based, end = coords.split("-") return chrom, int(start_1based) - 1, int(end) print(region_string_to_bed("chr1:1000-2000")) ``` Expected output: ``` ('chr1', 999, 2000) ``` That range is 1001 bases, not 1000. Verify any conversion by computing the length and confirming it matches your expectation. **Splitting on whitespace instead of tabs.** `line.split()` collapses runs of spaces and tabs together, so a file that is space-delimited parses without complaint and produces fields that look right until a name field contains a space. Always split on the literal tab character. **Assuming every row has 12 columns.** Rows in a single BED file legally vary in width from 3 to 12 columns, provided no optional column is skipped. Indexing `fields` without a length check raises an IndexError on BED3 files, and silently grabs the wrong field on BED6 files that also carry a name you did not expect. Always check `len(fields)` before accessing optional positions. **Skipping the track and browser line check.** A UCSC export that begins with a `track` line will crash a parser that calls `int(fields)` on it. Skipping lines whose first field is `track`, `browser`, or `#` costs nothing and prevents a whole class of failures. **Merging unsorted input.** `bedtools merge` and `pybedtools` merge require input sorted by chromosome and start. Merging unsorted input silently produces intervals that look plausible but miss overlaps. Sort first, every time, and remember that plain lexicographic chromosome sorting orders chr1, chr10, chr11, chr2. If chromosome order matters downstream, sort using a genome file or a natural chromosome ordering. **Confusing merge distance defaults.** By default, merging combines features that overlap or book-end each other with zero bases between them. If your intention was to join features within, say, 50 bases, you must pass an explicit distance flag, and failing to do so changes your covered-base totals substantially. State the distance you used in any methods description. **Loading a whole-genome file into a list.** A list of millions of namedtuples can exhaust memory on a laptop. Use the generator form of the parser, or aggregate on the fly, for any file above a few hundred megabytes uncompressed. **Ignoring strand in strand-aware comparisons.** A plus-strand feature and a minus-strand feature at identical coordinates are usually different biological things. Plain interval intersection ignores strand unless you tell the tool to consider it, so a naive overlap count can pair features that could never interact. **Treating score as a universal scale.** The specification describes score as an integer from 0 to 1000, but real files carry p-values, counts, or other quantities in that column with wildly different ranges. Never rank or threshold on score without knowing what the producing tool put there. A final limitation belongs in its own sentence rather than a bullet. Parsed and validated intervals are a computational result, and any analysis that will inform a publication, a clinical interpretation, or a downstream experiment should be reviewed by someone with domain expertise in the relevant assay and organism before it is treated as biological fact. ## Where to Go Next Once parsing is solid, the natural next steps build directly on the record structure you now have. - Interval arithmetic by hand. Write functions for overlap width, union, and subtraction. Overlap width between two intervals on the same chromosome is `min(a.end, b.end) - max(a.start, b.start)`, clamped at zero. Testing your implementation against pybedtools on small files is an excellent correctness exercise. - Converting between formats. BED to GTF, BED to a FASTA index query, or BED to a UCSC region string are all small, well-defined functions once coordinates are normalized. - Chunked genome scans. Combine the streaming parser with a sliding window to compute per-window GC content or coverage without loading a chromosome into memory. - Metadata-rich formats. If your analysis outgrows what BED columns can express, study how newer genomic formats model annotations, and how libraries abstract across multiple formats so that one reader can serve several file types [3]. The lesson from projects that solved this problem, such as the standardized repertoire format designed to unify the divergent outputs of immune sequencing tools [2], is that the cleaner your parsing layer is, the cheaper every later format conversion becomes. ## Frequently Asked Questions ### Is BED 0-based or 1-based? BED is 0-based with a half-open end, meaning the start coordinate counts from 0 and the end coordinate is one past the final base. A record with start 1000 and end 1500 covers bases 1001 through 1500 in 1-based terms and has a length of exactly 500 bases. ### How do I handle a BED file with fewer than 12 columns? Check the length of the fields list before accessing any optional column and substitute a default for missing ones. The specification requires that optional fields appear in order with no gaps, so defaults only ever apply to trailing fields. ### What is the difference between BED and GTF or GFF? BED is 0-based half-open and stores intervals with a simple column layout, while GTF and GFF are 1-based inclusive and store features with attributes in a key-value column. Converting between them always requires shifting the start coordinate by one. ### Can pandas read a gzipped BED file? Yes. Pass the `.gz` filename directly to `pd.read_csv` along with `sep="\t"` and `header=None`, and pandas decompresses the file on the fly. Add `comment="#"` to skip comment lines. ### Why does my merged BED file combine intervals that do not overlap? Because merging was performed on unsorted input, or because the merge distance was larger than you assumed. Always sort by chromosome and start position before merging, and pass an explicit distance value when you want gaps under a threshold to be joined. ### How do I convert a region string like chr1:1000-2000 into BED coordinates? ![Worked calculation converting the region string chr1:1000-2000 into BED start 999 and end 2000, giving a length of 1001 bases.](/article-images/illustrations/parse-bed-file-python-region-to-bed.webp) Subtract one from the start and leave the end unchanged, giving start 999 and end 2000. Region strings are 1-based and inclusive, while BED is 0-based half-open, so the same bases are described by different numbers. ### Does Biopython parse BED files? No. Biopython provides modules for sequence formats, alignments, structures, and tool interfaces, but it has no dedicated BED parser. Use the parsing approaches in this article, then bring the results into Biopython for sequence-level work [1]. ### How large a BED file can I parse in memory? It depends on available RAM and how much per-record data you keep. For files above a few hundred megabytes uncompressed, use a generator-based parser that processes one line at a time, or delegate interval operations to BEDTools through pybedtools, which streams rather than loads. ## Sources 1. [Biopython: freely available Python tools for computational molecular biology and bioinformatics.](https://pubmed.ncbi.nlm.nih.gov/19304878/) 2. [VDJML: a file format with tools for capturing the results of inferring immune receptor rearrangements.](https://pubmed.ncbi.nlm.nih.gov/27766961/) 3. [JPhyloIO: a Java library for event-based reading and writing of different phylogenetic file formats through a common interface.](https://pubmed.ncbi.nlm.nih.gov/31331268/)