step-2-intersect-genes

Metadata

Status	done
Assigned	`agent-15`
Agent identity	`3184716484e6f0ea08bb13539daf07686ee79d440505f1fdf2de0357707034c3`
Created	2026-03-31T21:03:36.016867686+00:00
Started	2026-03-31T21:08:30.389375549+00:00
Completed	2026-03-31T21:11:50.498248176+00:00
Tags	`impl`, `eval-scheduled`
Eval score	0.95
└ blocking impact	0.95
└ completeness	0.98
└ coordination overhead	0.93
└ correctness	0.95
└ downstream usability	0.96
└ efficiency	0.90
└ intent fidelity	0.88
└ style adherence	0.92

Description

Goal

Get the list of genes that fall within actual PHR intervals on CHM13.

Approach

Use bedtools intersect with the CHM13 gene annotation (RefSeq Liftoff or HPRC Liftoff) and the chm13.phrs.bed from Step 1.

zcat chm13v2.0_RefSeq_Liftoff_v5.2.gff3.gz \
  | awk '$3 == "gene"' \
  | bedtools intersect -a - -b chm13.phrs.bed -wa \
  > phrs.genes.gff3

Then extract gene names/IDs:

grep -oP 'Name=\K[^;]+' phrs.genes.gff3 > phrs.gene_names.txt
# Also extract Entrez/NCBI gene IDs if available for clusterProfiler

Context

Check the research task output for which annotation file to use
If the RefSeq Liftoff is insufficient, the research task should have identified the CHM13 file in Andrea's annotations at /moosefs/guarracino/HPRCv2/PHR_III/hprc_annotations/
Compare gene count to Angela's 327 unique genes — ours should be a subset since we're using tighter intervals

Validation

phrs.genes.gff3 exists with reasonable gene count
phrs.gene_names.txt has one gene per line, no empty lines
Log gene count, biotype breakdown if possible (protein-coding vs lncRNA vs pseudogene)
Briefly compare count to Angela's 327 unique genes

## Goal
Get the list of genes that fall within actual PHR intervals on CHM13.

## Approach
Use bedtools intersect with the CHM13 gene annotation (RefSeq Liftoff or HPRC Liftoff) and the `chm13.phrs.bed` from Step 1.

```bash
zcat chm13v2.0_RefSeq_Liftoff_v5.2.gff3.gz \
  | awk '$3 == "gene"' \
  | bedtools intersect -a - -b chm13.phrs.bed -wa \
  > phrs.genes.gff3
```

Then extract gene names/IDs:
```bash
grep -oP 'Name=\K[^;]+' phrs.genes.gff3 > phrs.gene_names.txt
# Also extract Entrez/NCBI gene IDs if available for clusterProfiler
```

## Context
- Check the research task output for which annotation file to use
- If the RefSeq Liftoff is insufficient, the research task should have identified the CHM13 file in Andrea's annotations at `/moosefs/guarracino/HPRCv2/PHR_III/hprc_annotations/`
- Compare gene count to Angela's 327 unique genes — ours should be a subset since we're using tighter intervals

## Validation
- `phrs.genes.gff3` exists with reasonable gene count
- `phrs.gene_names.txt` has one gene per line, no empty lines
- Log gene count, biotype breakdown if possible (protein-coding vs lncRNA vs pseudogene)
- Briefly compare count to Angela's 327 unique genes

Depends on

Required by

done step-3-run-go

Log

2026-03-31T21:03:36.015341845+00:00 Task paused
2026-03-31T21:04:09.339181799+00:00 Task published
2026-03-31T21:08:30.389379937+00:00 Spawned by coordinator --executor claude --model claude-sonnet-4-20250514
2026-03-31T21:09:04.845368561+00:00 Starting implementation - checking input files and performing gene intersection
2026-03-31T21:09:33.972613117+00:00 Found 61,312 total genes in annotation. Running bedtools intersect with PHR intervals
2026-03-31T21:10:47.460649218+00:00 Completed intersection: 412 gene entries (245 unique names), 70 unique Entrez IDs. Biotype breakdown: 204 pseudogenes, 108 lncRNAs, 51 miRNAs, 27 protein-coding, 21 transcribed_pseudogenes, 1 misc_RNA
2026-03-31T21:11:16.041415452+00:00 Validation completed: phrs.genes.gff3 exists (111KB, 412 entries), phrs.gene_names.txt has 412 lines with no empty lines. Found 245 unique gene names vs Angela's 327 - ours is smaller as expected with tighter intervals
2026-03-31T21:11:41.181487277+00:00 Committed: 7e2b9cd — pushed to remote
2026-03-31T21:11:50.498251773+00:00 Task marked as done