Metadata-Version: 2.4
Name: prokapeak
Version: 0.1.0
Summary: ProkaPeak: assign MACS3 peaks to target genes through TSS and gene-start windows
Author: Alexandre Dentler, Julian Trouillon
License-Expression: MIT
Project-URL: Repository, https://github.com/trouillon-lab/prokapeak
Keywords: ChIP-seq,DAP-seq,MACS3,transcription factor,operon,bacteria
Classifier: Programming Language :: Python :: 3
Classifier: Operating System :: OS Independent
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Bio-Informatics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.21
Requires-Dist: pandas>=1.4
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Dynamic: license-file

<img align="left" width="170" src="https://raw.githubusercontent.com/trouillon-lab/prokapeak/v0.1.0/images/ProkaPeak_logo.png" />

### A Python package to assign **Proka**ryotic **Peak**s to target genes.
#### ➡️ [_ProkaPeak web app_](https://prokapeak.trouillonlab.org/)
&nbsp;

&nbsp;
<br clear="left"/>

# ProkaPeak

ProkaPeak assigns transcription factor (TF) binding peaks from ChIP-seq, DAP-seq, or similar assays to the genes
they likely regulate. A peak is assigned to a gene when its summit lies in a window around that gene's start or 
transcription start site (TSS), depending on available annotation. 

ProkaPeak can run with only a peak list and gene annotation, but also supports 
both TSS and operon information for more precise predictions. If no operon information is provided, ProkaPeak can predict operons first.
ProkaPeak uses default window sizes and peak assignment rules chosen to maximize target gene recovery, as defined based on experimental
binding sites. All parameters are optionally modifiable.

### Web app
In addition to the Python package, we provide the [ProkaPeak web app](https://prokapeak.trouillonlab.org/), which allows for flexible peak assignment with visualisation of summary metrics. The web app allows to assign multiple peak files simultaneously, provides pre-made TSS and operon maps for selected bacterial species, and runs entirely locally in the user's browser without the need of data upload.

## Installation

```bash
pip install prokapeak
```

Python ≥ 3.9; installs numpy ≥ 1.21, pandas ≥ 1.4 and the commands `prokapeak` and `prokapeak-operons`. The
development version installs with `pip install git+https://github.com/trouillon-lab/prokapeak`.

## Usage

```bash
prokapeak --gff genome.gff --peaks peaks/*.narrowPeak --operons operons.tsv --tss tss.bed --out out/
```

Without installing, `python assign_targets.py` in a clone takes the same options (`python predict_operons.py` for
`prokapeak-operons`). Only `--gff`, `--peaks` and `--out` are required. `--predict-operons` replaces `--operons`
with operons predicted from gene spacing (see Predicted operons). Each peak file is processed in turn and gives one
output file.

`assign_targets` can also be imported: `read_genome` and `read_peaks` read the input files, `build_units` and `assign`
take the resulting tables; `main()` is the command line.

## Example

`example/` in the [repository](https://github.com/trouillon-lab/prokapeak) holds an artificial genome to show the file
formats: 100 genes on two chromosomes (`chrom1`, `plasmid1`), 45 operons, 28 TSSs and peak files for four made-up
TFs. `example/expected/` is the output of

```bash
cd example
prokapeak --gff genome.gff --operons operons.tsv --tss tss.bed --peaks peaks/*.narrowPeak --out out
```

which should match `expected/` file for file. `expected/predicted_operons.tsv` and `predicted_operon_pairs.tsv` are
the operons predicted from the same genome by

```bash
prokapeak-operons --gff genome.gff --out expected/predicted_operons.tsv --pairs expected/predicted_operon_pairs.tsv
```

(25 operons instead of the 45 the example was built with: a genome of 100 genes is far too small for the method,
see Predicted operons).

## Input files

All coordinates refer to the same genome; the chromosome name (column 1, or the GFF seqid) must be written
identically in every file. Genomes with several chromosomes or plasmids are supported. Files are plain text, named
by format: `.gff` or `.gff3`, `.tsv`, `.bed`, `.narrowPeak`. Lines starting with `#` are skipped and extra columns are
ignored; a missing file or a file of the wrong shape (too few columns, header lines, non-integer coordinates, missing
summits) stops the run with a message naming it.
Rows that cannot be matched to the annotation (loci with rows on several chromosomes or strands, peaks on unknown
chromosomes, operon genes not in the GFF, operons with genes on several chromosomes or strands, linked
TSSs whose gene is not found) are counted in warnings on stderr.

### `--gff` genome annotation (required)

GFF3. Rows of the listed feature types are read (`--gff-feature`, default `gene,CDS`); a locus with several rows,
of one type or of several, spans all of them. Every row needs the attribute holding the gene identifier
(`--gff-attribute`, default `locus_tag`), which is used in all other files and in the output.

| column | used |
|---|---|
| 1 seqid | chromosome |
| 3 type | one of `--gff-feature` (default `gene,CDS`) |
| 4 start, 5 end | 1-based, inclusive |
| 7 strand | `+` or `-`; the start codon is `start` on `+` and `end` on `-` |
| 9 attributes | `<--gff-attribute>=...` |

### `--peaks` peak files (required, one or more)

One file per TF, either MACS3 narrowPeak (`.narrowPeak`, 10 tab-separated columns, no header) or a BED file of
summits (`.bed`, at least 3 columns, every row a 1-bp interval at the summit, as in MACS3 `NAME_summits.bed`). Both
give the same assignments; a summit BED has no signal, p- or q-value, so those output columns are `.` (see Output).

| column | content |
|---|---|
| 1 chrom | chromosome |
| 2 start | 0-based peak start |
| 3 end | peak end |
| 4 name | peak id (copied to the output) |
| 5 score | copied to the output |
| 6 strand | ignored |
| 7 signal | fold enrichment, copied to the output |
| 8 pValue, 9 qValue | −log10, copied to the output |
| 10 peak | summit offset from column 2 |

The summit is column 2 + column 10 (0-based), or column 3 of a summit BED. The TF name is the file name without
extension, so file names must be distinct. An empty peak file gives an output file with only the header line.

### `--operons` operon map (optional)

Tab-separated, no header, two columns: operon id, locus tags of its genes separated by commas. One row per operon;
rows with the same operon id are merged, so one gene per row works too. The order of the genes in the list does
not matter: the genes of an operon are ordered 5'→3' from their GFF coordinates. Genes not listed are single-gene
units. Listed genes that are not in the GFF, and operons whose genes lie on several chromosomes or strands, are
skipped with a warning.

```
op1	b0002,b0003,b0004
op2	b0005
```

### `--tss` TSS map (optional)

BED6, tab-separated, no header, one row per TSS.

| column | content |
|---|---|
| 1 chrom | chromosome |
| 2 start, 3 end | the TSS as a 1-bp interval (0-based start; TSS position = end, 1-based) |
| 4 name | locus tag of the gene the TSS starts; `.` = not linked to a gene |
| 5 score | ignored |
| 6 strand | `+` or `-` |

A TSS is an anchor for the unit its gene starts: the operon or internal unit beginning with that gene, or the gene
itself. TSSs that are unlinked, or whose gene is not in the GFF or lies on another strand or chromosome, are skipped.

Example: `NC_000913.3  147  148  b0001  0  +`

## Assignment

**Units.** With `--operons`, every operon is a unit (operon start), and every gene of an operon other than the
first starts an internal unit: that gene and the genes downstream of it. Genes in no operon are single-gene
units. Without `--operons`, every gene is a unit.

**Windows** (bp, along the unit's strand, negative = upstream; summit position minus anchor position):

| unit | anchor | default | option |
|---|---|---|---|
| any unit with a linked TSS | the TSS | −200 to +100 | `--tss-window` |
| operon start or single gene without a TSS | start codon of the first gene | −400 to +25 | `--start-window` |
| internal start without a TSS | start codon of that gene | −250 to +50 | `--internal-window` |

Pass windows with `=`, as in `--tss-window=-150,100`, so the leading minus is not read as an option.

**Overlapping windows** (function `resolve` in the script).
1. A unit with several TSSs is assigned once, through the TSS nearest the summit.
2. Of several units on one strand, only the one with the nearest anchor is assigned (`--same-strand`).
3. A summit in windows on both strands (between divergent genes) is assigned only on the strand with the nearest
   anchor (`--opp-strand`), compared over the units kept by rule 2.
   Units or strands tied in distance are all kept.
4. An assigned unit assigns all its genes; a gene assigned through several units is listed once, through the
   nearest anchor.

## Predicted operons

Without an operon map, operons can be predicted from gene spacing using the method from UniOP (Su, Zhang and
Söding, bioRxiv 2024, doi:10.1101/2024.11.11.623000). Pairs at or above the cutoff (default 0.5) are joined; operons
of two or more genes are kept, in the operon map format.

The simplest is to predict them as part of the run, with `--predict-operons` in place of `--operons`:

```bash
prokapeak --gff genome.gff --peaks peaks/*.narrowPeak --predict-operons --out out/
```

`predicted_operons.tsv` and `predicted_operon_pairs.tsv` are written into `--out`; `--operon-cutoff` sets the cutoff.

Alternatively, `prokapeak-operons` predicts them on their own, to inspect them or pass them later to `--operons`:

```bash
prokapeak-operons --gff genome.gff --out operons.tsv [--pairs pairs.tsv] [--cutoff 0.5]
```

`--gff-feature` and `--gff-attribute` work as in `prokapeak`. The pair table has one row per pair of consecutive genes:

| column | content |
|---|---|
| chrom | chromosome |
| a, b | locus tags, in genome order |
| gap | start of b − end of a + 1 |
| class | `same`, `convergent` or `divergent` |
| prob | probability that a and b are in one operon; empty for opposite-strand pairs |

The probabilities equal those of the UniOP script (`tests/`), except that pairs are formed within a chromosome
only, the prior is clipped to 0.01–0.99, and a genome without convergent or without divergent pairs uses the
other class alone. The method needs a genome-sized gene set: below about 1,000 same-strand pairs the gaps stop
separating the pairs and nearly all same-strand neighbours are joined.

## Options

| option | effect |
|---|---|
| `--gff-feature` | GFF feature types to read, comma-separated (default `gene,CDS`); a locus spans all its rows |
| `--gff-attribute` | GFF attribute holding the gene id (default `locus_tag`) |
| `--predict-operons` | operons predicted from gene spacing instead of `--operons` |
| `--operon-cutoff` | pair probability that joins two genes for `--predict-operons` (default 0.5) |
| `--no-gene-start` | no start-codon windows; units without a TSS are not assigned |
| `--no-internal` | operon starts only, no internal units |
| `--same-strand closest\|both\|farthest` | a summit in windows of several units on one strand: assign the unit with the nearest anchor (default), all of them, or the one with the farthest anchor |
| `--opp-strand closest\|both\|farthest` | a summit in windows on both strands: assign the strand with the nearest anchor (default), both strands, or only the other strand |
| `--tss-window`, `--start-window`, `--internal-window` | window `up,down` in bp |

## Output

One `<TF>_assigned_peaks.bed` per peak file in `--out`: tab-separated, a `#` header line naming the columns, then
one row per peak and assigned gene, sorted by position. A peak assigned to an operon of four genes gives four rows;
a peak in no window is not listed. With `--predict-operons`, `--out` also holds `predicted_operons.tsv` and
`predicted_operon_pairs.tsv` (see Predicted operons).

| column | content |
|---|---|
| 1 chrom | chromosome of the peak |
| 2 start, 3 end | summit as a 1-bp interval (0-based start) |
| 4 target | locus tag of the assigned gene |
| 5 peak_score | column 5 of the peak file; `.` for a 3-column summit BED |
| 6 strand | strand of the target gene |
| 7 signal, 8 pscore, 9 qscore | narrowPeak columns 7–9; `.` for a summit BED |
| 10 target_chrom | chromosome of the target gene |
| 11 target_start, 12 target_end | target gene coordinates (0-based start) |
| 13 peak_id | column 4 of the peak file; `.` for a 3-column summit BED |
| 14 operon | operon id, `.` for a gene in no operon |
| 15 anchor | `TSS` or `start_codon` |
| 16 distance_to_anchor | summit minus anchor, along the target's strand (negative = upstream) |

`prokapeak` prints one line per TF: peaks assigned and number of target genes.

## License

ProkaPeak is released under the [MIT License](https://github.com/trouillon-lab/prokapeak/blob/main/LICENSE),
copyright (c) 2026 Trouillon lab.
