bio-genome-intervals-bed-file-basics

What this skill does

## Version Compatibility

Reference examples tested with: bedtools 2.31+, pybedtools 0.10+, pyranges 0.x (the `pyranges1` rewrite ships as a separate package), samtools 1.19+, pandas 2.2+, UCSC liftOver / CrossMap 0.7+ (the bare `CrossMap` entry point replaced `CrossMap.py` at 0.7.0).

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
- Python: `pip show <package>` then `help(module.function)` to check signatures

pyranges has a major-version API split: pyranges 0.x and the `pyranges1` rewrite differ in method names and DataFrame access; both keep `Chromosome/Start/End` columns and 0-based half-open coordinates. Check `import pyranges; pyranges.__version__` before chaining methods. Operations that need chromosome lengths (`slop`, `complement`, `shuffle`, `makewindows -g`, `-sorted` ordering) require a genome/chrom.sizes file. If code throws an error, introspect the installed tool and adapt rather than retrying.

# BED File Basics

**"Work with this interval file without shifting everything by one base"** -> Establish the coordinate convention from the format, read/create/validate/sort the intervals, and convert across format boundaries with the correct base shift.
- CLI: `bedtools sort -i in.bed`, `bedtools getfasta`, `bedtools makewindows -g genome.txt -w 10000`, `sort -k1,1 -k2,2n`
- Python: `pybedtools.BedTool('in.bed')`, `pr.read_bed('in.bed')` (pyranges), `pd.read_csv(sep='\t', comment='#')`

## The Single Most Important Modern Insight -- A Coordinate Is a Bare Integer With No Self-Describing Convention

A `start` column is just an `int`. Nothing in the file says whether it is 0-based or 1-based, so the convention lives in the analyst's head, keyed off the **format**, not the data. BED is 0-based half-open `[start, end)`; GTF/GFF, SAM, VCF, and wiggle are 1-based fully-closed `[start, end]`. Three load-bearing consequences:

1. **The conversion is `start - 1, end unchanged` -- and it throws no error if wrong.** A botched convention shift still parses, still runs, and silently shifts every answer by one base: gene bodies 1 bp short, boundary SNPs flipping in/out, exact-edge intersections toggling. The end is numerically identical between BED-half-open and GFF-closed because GFF's last *included* base and BED's first *excluded* position are the same boundary. The symmetry instinct (subtract 1 from both) is the classic bug. The reflex: convert `start_bed = start_1based - 1` (end unchanged) and **test the round-trip on a 1 bp feature** -- BED `chr1 5 6` == GFF `chr1 6 6`, both one base. Length is `end - start` in BED (no `+1`); `end - start + 1` in GTF.

2. **The two truly silent file failures.** (a) **Chrom-name mismatch** (`chr1` vs `1`, `chrM` vs `MT`): intersecting a `chr`-prefixed file against a bare-numeral one yields a perfectly valid **empty** result -- "no overlap" looks like biology, not a bug. Confirm shared naming (`cut -f1 a.bed | sort -u`) before any cross-file op. (b) **CRLF line endings** from Excel/Windows glue `\r` onto the last field (`end` becomes `"100\r"`); the tell is "works for some tools, breaks for others." `cat -A` shows `^M$`; fix with `dos2unix`. Never open a BED in Excel -- it date-mangles `SEPT9`->`9-Sep` and float-truncates large coordinates.

3. **bedtools `-sorted` assumes both inputs share the SAME chromosome order.** With a lexicographic (`chr1, chr10, chr2`) vs version (`chr1, chr2, chr10`) sort mismatch, modern bedtools (>=~2.25) detects the inconsistency and **errors out** (exit 1, `chromomsome sort ordering ... is inconsistent`); older versions silently swept past and **dropped chr10-chr22**. Either sort both files with the identical command, or pass `-g genome.txt` (derived from the same reference FASTA) to pin the expected order. For one-off work, omit `-sorted` (the in-memory path tolerates any order). The mismatch that stays SILENT on every version is a chromosome-NAME difference (`chr1` vs `1`), which returns an empty result with no error.

## Tool Taxonomy

| Tool | Role | Mechanism | When |
|------|------|-----------|------|
| bedtools | CLI interval algebra (Quinlan 2010 *Bioinformatics* 26:841) | streaming, sorted-input reference implementation | shell pipelines, large files, reproducible one-liners |
| pybedtools | Python wrapper over bedtools (Dale 2011 *Bioinformatics* 27:3423) | shells out per op; BedTool objects + iterators | inside a Python analysis; chaining with pandas |
| pyranges | pure-Python interval engine (Stovner 2020 *Bioinformatics* 36:918) | vectorized PyRanges/pandas, no bedtools binary | large in-memory joins, dataframe-native workflows |
| pandas | flat tabular read | `read_csv(sep='\t')`; knows NO coordinate semantics | quick filter/inspect; the analyst enforces 0-based + sort manually |
| UCSC bedToBigBed / tabix | indexed/compressed BED for random access | requires `sort -k1,1 -k2,2n`, no track lines | browser tracks, region queries on huge files |

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Quick create/sort/filter on the command line | bedtools + coreutils `sort -k1,1 -k2,2n` | no Python overhead; reproducible |
| Inside a pandas/Python pipeline | pybedtools or pyranges | stays in-process; pyranges if no bedtools binary |
| Convert VCF/GTF/SAM positions to BED | subtract 1 from start, end unchanged | the convention boundary; test on a 1 bp feature |
| Empty intersect / "no overlap found" | check chrom naming (`chr1` vs `1`) FIRST | the most common silent null result |
| Using `-sorted` for speed/RAM | sort both files identically, or pass `-g genome.txt` | modern bedtools errors on a lexicographic-vs-version mismatch; old versions dropped chroms silently |
| Need chromosome lengths (slop/complement/windows) | generate genome.txt from the SAME FASTA | a stale/generic chrom.sizes rots slop/complement |
| Set operations on these intervals | -> interval-arithmetic | this skill is the format/coordinate substrate |
| Parse a GTF/GFF gene model | -> gtf-gff-handling | 1-based, parent/child hierarchy, not a flat BED |
| Peaks not yet called | -> chip-seq/peak-calling or atac-seq/atac-peak-calling | this category operates on existing intervals |
| Convert between assemblies (hg19<->hg38) | liftOver/CrossMap, report unmapped | a different problem from convention shifts |

## BED Columns (BED3 -> BED12)

The first 3 fields are required; the rest are optional but **positional** (cannot supply field N without 1..N-1), and the field count must be identical on every line.

```
BED3   chrom  start  end
BED4   + name
BED5   + score (int 0-1000; '.' allowed)
BED6   + strand (+/-/.)            # the common stranded-interval form
BED12  + thickStart thickEnd itemRgb blockCount blockSizes blockStarts   # transcript/exon models
```

narrowPeak is **BED6+4** (`signalValue pValue qValue peak`); broadPeak is **BED6+3** (drops `peak`). The `peak` column is a **0-based offset from chromStart** (absolute summit = `chromStart + peak`), `-1` if none; `pValue`/`qValue` are `-log10` scaled with `-1` meaning "not assigned", NOT p=0.1.

## Create and Read BED Files

```python
import pybedtools
import pandas as pd

intervals = [('chr1', 100, 200, 'peak1', 100, '+'), ('chr1', 300, 400, 'peak2', 200, '-')]
bed = pybedtools.BedTool(intervals)                                       # from list of tuples
bed = pybedtools.BedTool.from_dataframe(pd.read_csv('peaks.tsv', sep='\t'))  # from a DataFrame
bed.saveas('peaks.bed')

for iv in pybedtools.BedTool('peaks.bed'):
    print(iv.chrom, iv.start, iv.end, len(iv))   # start/end are ints; len(iv) == end - start
df = pybedtools.BedTool('peaks.bed').to_dataframe(names=['chrom', 'start', 'end', 'name', 'score', 'strand'])
```

pandas reads BED as a flat table but knows nothing about coordinates: `pd.read_csv('in.bed', sep='\t', header=None, comment='#')` -- the analyst enforces 0-based and sorts manually.

## Generate the G