Pathway Analysis

Tests groups of genes for over-representation in KEGG pathways and Gene Ontology terms.

Overview

A list of genes is hard to interpret on its own. Over-representation testing asks whether a group contains more members of some pathway or ontology term than chance would give, turning a list of identifiers into statements about biology.

This module tests each group against the KEGG database and all three Gene Ontology branches using clusterProfiler. Input ENSEMBL gene IDs are first converted to ENTREZ IDs via org.Hs.eg.db; genes without a valid ENTREZ mapping are dropped, and group associations are preserved through the conversion. KEGG over-representation then runs through enrichKEGG() and GO through enrichGO() for Biological Process, Cellular Component and Molecular Function, both applying a hypergeometric test with Benjamini-Hochberg FDR correction.

Every group is analysed in the same run, and results from all analysis types are harmonised into one column schema - metadata columns for analysis type and ontology category, plus a category to subcategory to term classification - then combined into a single data frame saved as a compressed RDS. That shared schema is what lets groups such as cell types, conditions or treatments be compared side by side.

When to run it. After a set of genes of interest has been grouped, to ask which pathways and ontology terms those groups are enriched for.

Input

• --genes-file: a TSV of grouped genes with two required columns, group naming the category a gene belongs to and gene_id giving its ENSEMBL identifier. Every group present is tested. Example input/enrichment/protocol_example.pathway_genes.tsv:

  group  gene_id
  AD     ENSG00000139618
  AD     ENSG00000091831
  TL1    ENSG00000196839
  TL1    ENSG0000008123
  

• --name: the stem of the output file. Required.

• --organism: the KEGG organism code, hsa (human) by default.

• --pvalue-cutoff: the enrichment p-value threshold, 1 by default, which keeps every term so filtering can be done downstream.

• --cwd: the directory outputs are written to.

• --modular-script-dir: the code/script directory holding the analysis scripts.

• --numThreads, --job-size, --walltime and --mem: job resources.

Output

• {cwd}/pathway_analysis/{name}.combined_pathway_results.rds - one compressed RDS holding the KEGG and GO results for every gene group in a single data frame, 2842 rows and 17 columns on the example. SoS always inserts the step name, so results land in a pathway_analysis/ subdirectory of whatever --cwd is given. Example output/gsea/pathway_analysis/protocol_example.combined_pathway_results.rds:

  'data.frame':  2842 obs. of  17 variables:
   $ ID            : chr  "hsa01524" "hsa04115" "hsa05210" "hsa05162" ...
   $ Description   : chr  "Platinum drug resistance" "p53 signaling pathway" "Colorectal cancer" "Measles" ...
   $ GeneRatio     : chr  "7/9" "6/9" "6/9" "6/9" ...
   $ BgRatio       : chr  "75/9416" "75/9416" "87/9416" "136/9416" ...
   $ RichFactor    : num  0.0933 0.08 0.069 0.0441 0.0438 ...
   $ FoldEnrichment: num  97.6 83.7 72.2 46.2 45.8 ...
  

The 17 columns are the same for every analysis type, which is what makes groups directly comparable:

• ID and Description: the pathway or GO term identifier and its name.

• GeneRatio and BgRatio: input genes found in the term, against all database genes in that term.

• RichFactor and FoldEnrichment: the ratio of the two, raw and fold-scaled.

• pvalue, p.adjust and qvalue: the raw hypergeometric p-value, the Benjamini-Hochberg adjusted value, and the q-value for FDR control.

• group: which gene group the row came from.

• analysis_type and ont_category: KEGG or GO, and PATHWAY, BP, CC or MF.

• category and subcategory: the high-level and detailed classification of the term. Both are NA for KEGG rows in this example.

Minimal Working Example

KEGG enrichment queries the online KEGG API, so this step needs network access. GO enrichment uses the local org.Hs.eg.db and works offline.

Timing: ~1-2 min (on toy dataset)

sos run pipeline/gsea.ipynb pathway_analysis \
    --genes_file input/enrichment/protocol_example.pathway_genes.tsv \
    --name protocol_example \
    --pvalue_cutoff 1 --organism hsa \
    --cwd output/gsea

Command Interface

sos run pipeline/gsea.ipynb -h

usage: sos run code/SoS/enrichment/gsea.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:
  pathway_analysis

Global Workflow Options:
  --cwd output (as path)
                        Path to the work directory of the analysis.
  --modular-script-dir code/script (as path)
  --genes-file . (as path)
  --name VAL (as str, required)
  --numThreads 8 (as int)
                        Number of threads
  --job-size 1 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 12h
  --mem 16G

Sections
  pathway_analysis:
    Workflow Options:
      --pvalue-cutoff 1 (as int)
      --organism hsa

Workflow implementation

[global]
# Path to the work directory of the analysis.
parameter: cwd = path('output')
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir

parameter: genes_file = path()  # TSV file with columns: group, gene_id
parameter: name = str
# Number of threads
parameter: numThreads = 8
# For cluster jobs, number commands to run per job
parameter: job_size = 1
parameter: walltime = '12h'
parameter: mem = '16G'

[pathway_analysis]
parameter: pvalue_cutoff = 1
parameter: organism = 'hsa'
output: pathway_results = f'{cwd:a}/{step_name}/{name}.combined_pathway_results.rds'

task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bnn}'
bash: expand = '${ }', stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout'
    Rscript ${modular_script_dir}/enrichment/gsea.R \
        --genes-file ${genes_file} \
        --organism ${organism} \
        --pvalue-cutoff ${pvalue_cutoff} \
        --output ${_output['pathway_results']}