annotshift

annotshift shifts genomic data between whole-genome coordinates and gene-defined zone coordinates.

The to-zones command defines genomic zones from one or more GFF3 annotations, extracts the corresponding FASTA sequences, shifts each GFF3 to zone-relative coordinates, and optionally extracts and shifts an indexed BAM.

The to-genome command shifts a zone-relative GFF3 back to whole-genome coordinates.

Regional sequences are never reverse-complemented. Strand and CDS phase are therefore preserved unchanged.

Input GFF table

The source table is tab-separated and must contain the following columns:

gff	use_for_zone_definition
path/to/annot1.gff	true
path/to/annot2.gff	true
path/to/annot3.gff	false

Relative GFF paths are resolved relative to the TSV file.

Zone definition

For every GFF with use_for_zone_definition=true:

1. Every gene feature is retained.
2. A zone starts at gene.start - flank and ends at gene.end + flank, clipped to chromosome boundaries.
3. Defining genes are processed by genomic start. A next gene joins the current zone only when next_gene.start <= current_zone.end.
4. Therefore, with flank=1000, a next gene starting at previous_gene.end + 1000 joins the zone, whereas one starting at previous_gene.end + 1001 does not.

For every GFF with use_for_zone_definition=false:

1. The coding span of each gene is defined as min(CDS.start)..max(CDS.end) across all CDS features and transcripts belonging to that gene.
2. Genes without CDS features are ignored.
3. A gene is retained when its coding span overlaps the coding span of at least one defining gene.
4. A retained false gene expands its associated zone enough to contain the complete gene feature, without adding additional flank around that gene.
5. A retained false gene can merge defining groups when its complete gene interval reaches a defining gene from another group. Overlapping selected false-gene intervals can also bridge defining groups.
6. Defining zones that remain separate under the explicit next_gene.start <= current_zone.end rule are not merged merely because their flanking regions overlap. Consequently, distinct output zones can share genomic sequence.

Outputs

For an input FASTA named genome.fasta:

output/
├── genome.zones.fasta
├── genome.zones.fasta.fai
├── zones.tsv
├── gff/
│   ├── annot1.gff
│   ├── annot2.gff
│   └── annot3.gff
├── <input-bam-stem>.zones.bam      # only with --bam
└── <input-bam-stem>.zones.bam.bai

The input FASTA extension is preserved. For example:

genome.fa    -> genome.zones.fa
genome.fasta -> genome.zones.fasta

Zone reference names are encoded as:

<source_seqid>_<source_start>_<source_end>

For example:

Chr1A_1500000_1600000

zones.tsv contains the same mapping explicitly.

Output GFF files contain only feature lines. GFF comments and directives beginning with # are omitted.

BAM behavior

The input BAM must be coordinate-sorted and indexed.

Extraction is performed with indexed regional fetch() calls, one final zone at a time. Only alignments whose complete reference-consuming span is contained within the zone are retained.

For paired-end data, retained alignments are buffered only for the current zone:

• if both primary R1 and R2 alignments are retained in the same zone, mate coordinates and template length are preserved after coordinate shifting;
• otherwise, the retained alignment remains marked as originating from a paired template, but proper pair is cleared and RNEXT, PNEXT, and TLEN are made unavailable;
• secondary and supplementary alignments are retained when fully contained within the zone;
• coordinate-dependent SA, CC, and CP tags are removed because their whole-genome coordinates would otherwise become stale.

Zones are written in output-reference order, so the output BAM is produced directly in coordinate order before indexing.

Because the zone-definition rule can intentionally produce overlapping regional sequences, an input alignment that is fully contained within more than one output zone is projected once into each relevant zone. This is useful for independent regional visualization, but the resulting BAM should not be used to sum read counts across all zones without accounting for duplicated projections.

Requirements

The input FASTA must already have a .fai index:

samtools faidx genome.fa

The optional BAM must also already have a BAI or CSI index.

Install the package with:

python3 -m pip install --user .

For development and testing:

python3 -m pip install --user -e '.[test]'

Installing only the dependencies with:

python3 -m pip install --user -r requirements.txt

does not install the annotshift package itself.

Shift whole-genome data to zone coordinates

annotshift to-zones \
    --fasta genome.fa \
    --gffs annotations.tsv \
    --flank 10000 \
    --bam expression.sorted.bam \
    --output-dir output \
    --threads 8

The BAM input is optional:

annotshift to-zones \
    --fasta genome.fa \
    --gffs annotations.tsv \
    --flank 10000 \
    --output-dir output

The module entry point is equivalent:

python3 -m annotshift to-zones ...

Shift a zone-relative GFF to whole-genome coordinates

The inverse operation uses the coordinates encoded in each zone sequence name and does not require zones.tsv.

annotshift to-genome \
    --in-gff output/gff/annot1.gff \
    --out-gff annot1.genome.gff

For a feature on Chr1A_1500000_1600000, whole-genome coordinates are calculated as:

genome_start = 1500000 + zone_start - 1
genome_end   = 1500000 + zone_end - 1

The module entry point is equivalent:

python3 -m annotshift to-genome ...

Tests

Install the test dependencies and run:

python3 -m pip install --user -e '.[test]'
pytest -q

The tests cover:

• zone boundary rules;
• intentionally overlapping flanking regions;
• CDS-envelope selection;
• false-gene zone bridging;
• GFF hierarchy and strand preservation;
• shifting GFF coordinates back to the whole genome;
• spliced BAM alignments;
• boundary-crossing alignments;
• complete paired-end reads;
• orphaned paired reads.