Gene Coordinate Annotation

Attaches genomic coordinates to a molecular phenotype matrix, and emits the region and gene lists that downstream steps split work by.

Overview¶

Molecular phenotype matrices arrive keyed by feature ID alone -- an ENSEMBL gene ID, a UniProt accession, a leafcutter intron-cluster label -- but every downstream QTL step needs a genomic position to define a cis window. This module attaches those coordinates, turning a plain matrix into a coordinate-sorted, bgzipped and tabix-indexed bed.gz file.

Coordinates come from a collapsed gene-model GTF, following the GTEx pipeline convention, so the annotation used here matches the one used to build the GTF in the first place. Feature IDs that the GTF does not contain are dropped rather than guessed at, which is why the row count of the output can be lower than that of the input.

Which route applies depends on how the features are named. Gene-level expression matrices match the GTF directly. Protein matrices carry gene_id|UniProt composite IDs and need an auxiliary mapping to reach a gene. Leafcutter intron clusters are not genes at all: they are genomic intervals that must first be mapped onto the genes they overlap, which also produces the phenotype-group file that grouped multiple-testing correction later consumes.

An older biomaRt route is retained for reference. It places a single-base TSS interval rather than the full gene body and is obsolete; new analyses should use the GTF-based route.

When to run it. After a molecular phenotype matrix has been assembled and normalised, and before any QTL association scan, since the association modules read the coordinates from this bed.gz.

Input¶

• --phenoFile (required): the molecular phenotype matrix to annotate. Feature IDs live in the column named by --phenotype-id-column; all remaining columns are samples.

• --coordinate-annotation (required): the collapsed gene-model GTF (or, for the biomaRt route, the reference the annotation is drawn from) supplying chromosome, start and end per feature.

• --phenotype-id-column (default gene_id): name of the feature-ID column in --phenoFile.

• --molecular-trait-type (default gene): gene for expression, protein for proteomics; selects how the feature ID is parsed before the GTF lookup.

• --auxiliary-id-mapping (default .): two-column map used for the protein trait type to reach a gene from a UniProt accession.

• --sample-participant-lookup (default .): optional two-column sample-to-participant map; supply it to rename matrix columns to participant IDs in the output.

• --ensembl-version: ENSEMBL release passed to biomaRt; only the annotate_coord_biomart workflow reads it, and it must be set explicitly there.

• --intron-count: leafcutter intron-count table, read by map_leafcutter_cluster_to_gene to enumerate the clusters that need a gene assignment.

• --map-stra / --cwd / --container / --entrypoint: intron-to-gene matching mode for cluster mapping (default site), output directory (default output), and the software environment.

Example inputs shipped with the protocol:

• tests/fixtures/gene_annotation/protocol_example.rnaseq.bed.gz -- toy gene-level bulk RNA-seq matrix.

• tests/fixtures/gene_annotation/protocol_example.rnaseq.gene_ID.tsv -- the same matrix with a gene_ID column, for the biomaRt example.

• tests/fixtures/gene_annotation/protocol_example.protein.no_coord.tsv -- toy protein matrix, gene_id|UniProt IDs, no coordinates.

• input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene -- the collapsed gene model.

First rows of the expression matrix:

#chr	start	end	ID	SAMPLE_001
chr22	10939387	10961337	ENSG00000283047	37.077
chr22	15528191	15529138	ENSG00000130538	0.0

Output¶

• <phenoFile stem>.bed.gz plus its .tbi index -- the annotated matrix, coordinate-sorted, with #chr start end ID prepended to the sample columns. Written by every annotation workflow.

• <phenoFile stem>.region_list.txt -- one row per retained feature giving its coordinates, ID and the path of the bed.gz it came from; QTL modules read this to enumerate cis regions.

• <phenoFile stem>.gene_list.tsv -- the feature IDs that were successfully matched to the GTF.

• *.leafcutter.clusters_to_genes.txt and *.exon_list -- from map_leafcutter_cluster_to_gene: the cluster-to-gene assignment and the exon table it was derived from.

• *.phenotype_group.txt -- from the leafcutter isoform workflow: the isoform-to-gene grouping used by grouped multiple-testing correction.

tests/fixtures/gene_annotation/expected/protocol_example.rnaseq.bed.bed.gz:

#chr	start	end	ID	strand	SAMPLE_001
chr22	10939387	10961337	ENSG00000283047	-	37.0768403309204
chr22	15528191	15529138	ENSG00000130538	+	0.0

tests/fixtures/gene_annotation/expected/protocol_example.rnaseq.bed.region_list.txt (the path column holds the absolute location of the bed.gz, abbreviated here):

#chr	start	end	ID	strand	path
chr22	10939387	10961337	ENSG00000283047	-	.../output/gene_annotation/protocol_example.rnaseq.bed.bed.gz
chr22	15528191	15529138	ENSG00000130538	+	.../output/gene_annotation/protocol_example.rnaseq.bed.bed.gz

Minimal Working Example¶

Each block below is a self-contained route: pick the one that matches how the features in the matrix are named. All of them read from input/ and write to output/gene_annotation/, and none depends on the others having been run.

Gene Expression¶

Annotate a gene-level expression matrix. Each feature ID is an ENSEMBL gene ID that is matched directly against the GTF to obtain its chromosome, start and end coordinates.

Timing: ~1 min (on toy dataset)

The output protocol_example.rnaseq.bed.gz is written to output/gene_annotation/, with genomic coordinates (#chr start end ID) prepended to the expression matrix.

Proteins¶

Annotate a protein matrix. Here each feature ID has the form gene_id\|UniProt; the gene ID portion is matched against the GTF for coordinates and the UniProt portion is retained in the output feature ID. Use --molecular-trait-type protein to select this behaviour.

Timing: ~1 min (on toy dataset)

Leafcutter clusters to genes¶

For leafcutter splicing data, intron clusters first need to be assigned to genes. The input is a leafcutter intron-count table whose feature IDs have the form chr:start:end:clu_N_strand; --coordinate-annotation is a full gene-model GTF (with exon records). The default --map-stra site maps introns to genes by matching splice-site coordinates to exon boundaries.

Timing: ~2 min (on toy dataset)

Leafcutter isoforms¶

This example turns the raw leafcutter intron-excision output into a phenotype matrix that TensorQTL can use. Building on the cluster-to-gene mapping from the previous example, it attaches genomic coordinates and a gene assignment to every intron, then writes a coordinate-annotated phenotype BED together with a phenotype-group file that links each intron back to its gene. The same --intron-count table and --coordinate-annotation GTF used above are supplied here.

Along the way the step produces three intermediate files, shown below with a few example rows from the toy data.

1. Exon list — exon records parsed from the annotation GTF, used to match intron splice sites to genes (<intron_count>.exon_list):

2. Clusters-to-genes — each leafcutter cluster mapped to the gene it belongs to (<intron_count>.leafcutter.clusters_to_genes.txt):

3. Phenotype group — links each annotated intron feature to its gene, so introns from the same gene are grouped together (<phenotype>.phenotype_group.txt):

Timing: ~2 min (on toy dataset)

Alternative coordinate annotation with biomaRt¶

Instead of a local GTF, coordinates can be fetched from the Ensembl biomaRt web service. The phenotype matrix must have a gene_ID column of ENSEMBL gene IDs. --ensembl-version selects the Ensembl release; use a release whose server is reachable from your network (e.g. a current release). This step requires internet access to Ensembl.

Timing: ~2 min (on toy dataset)

phenotype_group

Command Interface¶

Workflow implementation¶

Implementation using pyqtl¶

Implementation based on GTEx pipeline.

Following step serves to annotate coord for gene expression file.

Implementation using biomaRt¶

This workflow adds the annotations of chr pos(TSS where start = end -1) and gene_ID to the bed file. This workflow is obsolete.

Annotation of leafcutter isoform¶

The gtf used here should be the collapsed gtf, i.e. the final output of reference_data gtf processing and the one used to called rnaseq.