Gene Coordinate Annotation

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. Psichomics PSI output is already grouped by gene symbol, so it only needs that symbol translated to an ENSEMBL ID.

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:

input/rnaseq/protocol_example.rnaseq.bed.gz – toy gene-level bulk RNA-seq matrix.
input/rnaseq/protocol_example.rnaseq.gene_ID.tsv – the same matrix with a gene_ID column, for the biomaRt example.
input/proteomics/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 and psichomics isoform workflows: the isoform-to-gene grouping used by grouped multiple-testing correction.

output/gene_annotation/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

output/gene_annotation/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)

sos run pipeline/gene_annotation.ipynb annotate_coord \
    --cwd output/gene_annotation \
    --phenoFile input/rnaseq/protocol_example.rnaseq.bed.gz \
    --coordinate-annotation input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf \
    --phenotype-id-column gene_id

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)

sos run pipeline/gene_annotation.ipynb annotate_coord \
    --cwd output/gene_annotation \
    --phenoFile input/proteomics/protocol_example.protein.no_coord.tsv \
    --coordinate-annotation input/reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf \
    --phenotype-id-column gene_id \
    --molecular-trait-type protein

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)

sos run pipeline/gene_annotation.ipynb map_leafcutter_cluster_to_gene \
    --cwd output/gene_annotation \
    --phenoFile input/rnaseq/protocol_example.leafcutter.phenotype.bed.gz \
    --intron-count input/rnaseq/protocol_example.leafcutter.intron_count.tsv \
    --coordinate-annotation input/reference_data/Homo_sapiens.GRCh38.103.chr.gtf \
    --map-stra site

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):

chr	start	end	strand	gene_id	gene_name
1	111869	112227	+	ENSG00000223972	DDX11L1
1	112613	112721	+	ENSG00000223972	DDX11L1

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

clu	genes
22:clu_10_+	ENSG00000236052
22:clu_11_-	ENSG00000185264

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

ID	gene
22:19149095:19149663:clu_14_-:ENSG00000063515	ENSG00000063515
22:19149916:19150025:clu_14_-:ENSG00000063515	ENSG00000063515

Timing: ~2 min (on toy dataset)

sos run pipeline/gene_annotation.ipynb annotate_leafcutter_isoforms \
    --cwd output/gene_annotation \
    --phenoFile input/rnaseq/protocol_example.leafcutter.phenotype.bed.gz \
    --intron-count input/rnaseq/protocol_example.leafcutter.intron_count.tsv \
    --coordinate-annotation input/reference_data/Homo_sapiens.GRCh38.103.chr.gtf \
    --map-stra site

Psichomics isoforms#

For psichomics splicing quantifications, each event ID ends in _<gene_id>; the gene ID is matched against the GTF to obtain coordinates. The output is a coordinate-annotated phenotype BED plus a phenotype-group file.

Timing: ~2 min (on toy dataset)

sos run pipeline/gene_annotation.ipynb annotate_psichomics_isoforms \
    --cwd output/gene_annotation \
    --phenoFile input/rnaseq/protocol_example.psichomics.phenotype.tsv \
    --coordinate-annotation input/reference_data/Homo_sapiens.GRCh38.103.chr.gtf

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)

sos run pipeline/gene_annotation.ipynb annotate_coord_biomart \
    --cwd output/gene_annotation \
    --phenoFile input/rnaseq/protocol_example.rnaseq.gene_ID.tsv \
    --ensembl-version 115

chr	start	end	strand	gene_id	gene_name
chr1	29554	30039	+	ENSG00000243485	MIR1302-2HG
chr1	30564	30667	+	ENSG00000243485	MIR1302-2HG
chr1	30976	31097	+	ENSG00000243485	MIR1302-2HG
chr1	35721	36081	-	ENSG00000237613	FAM138A
chr1	35277	35481	-	ENSG00000237613	FAM138A
chr1	34554	35174	-	ENSG00000237613	FAM138A
chr1	65419	65433	+	ENSG00000186092	OR4F5
chr1	65520	65573	+	ENSG00000186092	OR4F5
chr1	69037	71585	+	ENSG00000186092	OR4F5

phenotype_group

X1	X2
7:102476270:102478811:clu_309_-:ENSG00000005075	ENSG00000005075
7:102476270:102478808:clu_309_-:ENSG00000005075	ENSG00000005075
X:47572961:47574002:clu_349_-:ENSG00000008056	ENSG00000008056
X:47572999:47574002:clu_349_-:ENSG00000008056	ENSG00000008056
8:27236905:27239971:clu_322_-:ENSG00000015592	ENSG00000015592
8:27239279:27239971:clu_322_-:ENSG00000015592	ENSG00000015592
8:27241262:27241677:clu_323_-:ENSG00000015592	ENSG00000015592
8:27241262:27242397:clu_323_-:ENSG00000015592	ENSG00000015592
8:27241757:27242397:clu_323_-:ENSG00000015592	ENSG00000015592
1:35558223:35559107:clu_4_+:ENSG00000020129	ENSG00000020129

Command Interface#

sos run pipeline/gene_annotation.ipynb -h

usage: sos run pipeline/gene_annotation.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  annotate_coord
  annotate_coord_biomart
  map_leafcutter_cluster_to_gene
  annotate_leafcutter_isoforms
  annotate_psichomics_isoforms

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --cwd output (as path)
                        Work directory & output directory
  --phenoFile VAL (as path, required)
                        Molecular phenotype matrix
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 5h
                        Wall clock time expected
  --mem 16G
                        Memory expected
  --numThreads 1 (as int)
                        Number of threads
  --container ''
  --entrypoint ''

Sections
  annotate_coord:
    Workflow Options:
      --sample-participant-lookup . (as path)
                        A file to map sample ID from expression to genotype,
                        must contain two columns, sample_id and participant_id,
                        mapping IDs in the expression files to IDs in the
                        genotype (these can be the same).
      --molecular-trait-type gene
                        Options: gene, protein, atac
      --coordinate-annotation VAL (as path, required)
                        gtf annotation for RNA-seq data; other types of
                        annotation index for protein and atac-seq
      --phenotype-id-column 'gene_id'
                        gene_id or gene_name for RNA-seq data, SOMAseqID for
                        proteins for example
      --auxiliary-id-mapping . (as path)
                        Optional mapping file (protein_name_index for protein,
                        etc.)
      --[no-]strip-id (default to False)
                        Optional processing flag
      --sep '\t'
  annotate_coord_biomart:
    Workflow Options:
      --ensembl-version VAL (as int, required)
  map_leafcutter_cluster_to_gene:
    Workflow Options:
      --intron-count VAL (as path, required)
                        Extract the code in case psichromatic needs to be
                        processed the same way PheoFile in this step is the
                        intron_count file
      --coordinate-annotation VAL (as path, required)
                        gtf annotation with exons (the full gene model), used to
                        map intron splice sites to genes
      --map-stra site
                        Defines the mapping strategy options: 'site' or
                        'region', with 'site' as the default. The 'site'
                        strategy maps introns to the start and end of each exon.
                        The 'region' strategy, to be recommended in leafcutter2,
                        maps each intron based on it overlapping more than
                        overlap_ratio of a gene's region.
      --overlap-ratio 0.8 (as float)
                        Define the overlap ratio as the proportion of the
                        cluster length that intersects with a gene, used to
                        determine mapping to the gene.
  annotate_leafcutter_isoforms:
    Workflow Options:
      --sample-participant-lookup . (as path)
      --coordinate-annotation VAL (as path, required)
                        gtf annotation with exons (the full gene model)
  annotate_psichomics_isoforms:
    Workflow Options:
      --sample-participant-lookup . (as path)
      --coordinate-annotation VAL (as path, required)
                        gtf annotation with exons (the full gene model)

Workflow implementation#

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# Work directory & output directory
parameter: cwd = path("output")
# Molecular phenotype matrix
parameter: phenoFile = path
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Wall clock time expected
parameter: walltime = "5h"
# Memory expected
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 1
parameter: container = ""
parameter: entrypoint= ""

Implementation using `pyqtl`#

Implementation based on GTEx pipeline.

Following step serves to annotate coord for gene expression file.

[annotate_coord]
# A file to map sample ID from expression to genotype, must contain two columns, sample_id and participant_id, mapping IDs in the expression files to IDs in the genotype (these can be the same).
parameter: sample_participant_lookup = path()
# Options: gene, protein, atac
parameter: molecular_trait_type = "gene" 
# gtf annotation for RNA-seq data; other types of annotation index for protein and atac-seq
parameter: coordinate_annotation = path
# gene_id or gene_name for RNA-seq data, SOMAseqID for proteins for example
parameter: phenotype_id_column = "gene_id"
# Optional mapping file (protein_name_index for protein, etc.)
parameter: auxiliary_id_mapping = path()
# Optional processing flag
parameter: strip_id = False
parameter: sep = "\t"

input: phenoFile, coordinate_annotation
output: f'{cwd:a}/{_input[0]:bn}.bed.gz', f'{cwd:a}/{_input[0]:bn}.region_list.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/gene_annotation.R \
        --step annotate_coord \
        --phenoFile "${_input[0]}" \
        --coordinate-annotation "${_input[1]}" \
        --sample-participant-lookup "${sample_participant_lookup}" \
        --molecular-trait-type "${molecular_trait_type}" \
        --phenotype-id-column "${phenotype_id_column}" \
        --auxiliary-id-mapping "${auxiliary_id_mapping}" \
        --sep "${sep}" \
        ${"--strip-id" if strip_id else ""} \
        --output-bed "${_output[0]}" \
        --output-region-list "${_output[1]}" \
        --gene-list-output "${cwd:a}/${_input[0]:bn}.gene_list.tsv"

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.

[annotate_coord_biomart]
parameter: ensembl_version=int
input: phenoFile
output: f'{cwd:a}/{_input:bn}.bed.gz',
        f'{cwd:a}/{_input:bn}.region_list.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, tags = f'{step_name}_{_output[0]:bn}'  
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/gene_annotation.R \
        --step annotate_coord_biomart \
        --phenoFile "${_input[0]}" \
        --ensembl-version ${ensembl_version} \
        --output-bed "${_output[0]}" \
        --output-region-list "${_output[1]}"

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.

[map_leafcutter_cluster_to_gene]
## Extract the code in case psichromatic needs to be processed the same way
## PheoFile in this step is the intron_count file
parameter: intron_count = path
# gtf annotation with exons (the full gene model), used to map intron splice sites to genes
parameter: coordinate_annotation = path
# Defines the mapping strategy options: 'site' or 'region', with 'site' as the default. 
# The 'site' strategy maps introns to the start and end of each exon. 
# The 'region' strategy, to be recommended in leafcutter2, maps each intron based on it overlapping more than overlap_ratio of a gene's region.
parameter: map_stra = "site" 
# Define the overlap ratio as the proportion of the cluster length that intersects with a gene, used to determine mapping to the gene.
parameter: overlap_ratio = 0.8
input: intron_count, coordinate_annotation
output: f'{cwd}/{_input[0]:b}.exon_list', f'{cwd}/{_input[0]:b}.leafcutter.clusters_to_genes.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, tags = f'{step_name}_{_output[0]:bn}'  
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/gene_annotation.R \
        --step map_leafcutter_cluster_to_gene \
        --intron-count "${_input[0]}" \
        --annotation-gtf "${_input[1]}" \
        --map-stra "${map_stra}" \
        --overlap-ratio ${overlap_ratio} \
        --output-exon-list "${_output[0]}" \
        --output-cluster-map "${_output[1]}"

[annotate_leafcutter_isoforms]
parameter: sample_participant_lookup = path()
# gtf annotation with exons (the full gene model)
parameter: coordinate_annotation = path
input: phenoFile, coordinate_annotation, output_from("map_leafcutter_cluster_to_gene")
output: f'{cwd:a}/{_input[0]:bn}.formated.bed.gz', f'{cwd:a}/{_input[0]:bn}.phenotype_group.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, tags = f'{step_name}_{_output[0]:bn}'  
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/gene_annotation.R \
        --step annotate_leafcutter_isoforms \
        --phenoFile "${_input[0]}" \
        --annotation-gtf "${_input[1]}" \
        --cluster-map "${_input[3]}" \
        --sample-participant-lookup "${sample_participant_lookup}" \
        --output-bed "${_output[0]}" \
        --output-phenotype-group "${_output[1]}"

Processing of psichomics output#

It occurs that the psichomatic by default grouped the isoforms by gene name, so only thing needs to be done is to extract this information and potentially renamed the gene symbol into ENSG ID

[annotate_psichomics_isoforms]
parameter: sample_participant_lookup = path()
# gtf annotation with exons (the full gene model)
parameter: coordinate_annotation = path
input: phenoFile, coordinate_annotation
output: f'{cwd:a}/{_input[0]:bn}.formated.bed.gz', f'{cwd:a}/{_input[0]:bn}.phenotype_group.txt'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime,  mem = mem, tags = f'{step_name}_{_output[0]:bn}'  
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout', container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/data_preprocessing/phenotype/gene_annotation.R \
        --step annotate_psichomics_isoforms \
        --phenoFile "${_input[0]}" \
        --annotation-gtf "${_input[1]}" \
        --sample-participant-lookup "${sample_participant_lookup}" \
        --output-bed "${_output[0]}" \
        --output-phenotype-group "${_output[1]}"

Gene Coordinate Annotation

Contents

Gene Coordinate Annotation#

Overview#

Input#

Output#

Minimal Working Example#

Gene Expression#

Proteins#

Leafcutter clusters to genes#

Leafcutter isoforms#

Psichomics isoforms#

Alternative coordinate annotation with biomaRt#

Command Interface#

Workflow implementation#

Implementation using pyqtl#

Implementation using biomaRt#

Annotation of leafcutter isoform#

Processing of psichomics output#

Implementation using `pyqtl`#