GREGOR enrichment analysis

GREGOR enrichment analysis#

Tests whether a set of trait-associated variants falls inside annotated regulatory features more often than a matched control set does.

Overview#

A set of trait-associated variants is more interpretable if you can say what kind of sequence they fall in. Enrichment testing asks whether they overlap a class of genomic feature - an annotation, a chromatin state, a set of regulatory elements - more often than chance allows, where chance has to account for the fact that variants are not exchangeable: they differ in minor allele frequency, in the number of LD proxies they carry, and in distance to the nearest gene. GREGOR builds matched control variants on exactly those properties - minor allele frequency, LD proxy count, distance to the nearest TSS, and local gene density - so the comparison is fair.

GREGOR’s own PValue column, defined in the GREGOR paper, sometimes comes back greater than 1, and a p-value says less about a result than the size of the enrichment does. This module therefore parses GREGOR’s intermediate output into a 2x2 table - inside versus outside the annotation peak, positive set versus matched negative set - and runs Fisher’s exact test on it to recover both a p-value and an odds ratio.

When to run it. After you have a variant set worth characterising - fine-mapped credible sets, or GWAS hits - and an annotation you want to test them against.

Input#

--index_snp_file: the positive variant set, one chr:pos per line with no header, on the same build as the reference database (hg19 for the standard ones). Example input/ld/protocol_example.index.snps.txt:
```
chr1:752566
chr1:776546
chr2:233567
```
--bed_file_index: a two-column tab-separated index naming each annotation and pointing at its bed file. Example input/enrichment/protocol_example.bed.file.index:
```
annotationA    <path>/annotationA.bed
annotationB    <path>/annotationB.bed
annotationC    <path>/annotationC.bed
```
--gregor_db: the GREGOR reference database directory, one SQLite CUBE table per chromosome. The full per-population reference is around 20 GB and is distributed under the University of Michigan license, so it is not bundled here and input/enrichment/protocol_example.gregor_ref is a placeholder. A downsized chr22 EUR slice is committed under tests/fixtures/gregor/EUR so the notebook can be smoke-tested end to end.
--pop: the reference population, EUR by default. It has to match the population the database was built for.
--r2_threshold and --ld_window_size: what counts as an LD buddy, 0.7 and 10000 by default. Both have to match the values the reference database was built with, 1 Mb for the standard r2 >= 0.7 databases.
--fisher1 and --fisher2, used by gregor_fisher_plot only: two enrichment-result tables to plot against each other, in the form the gregor workflow writes. Example input/enrichment/protocol_example.trait1_enrichment_results.txt:
```
Bed_File              odds   low    high   p_fisher
Promoter_UCSC.bed     1.377  0.751  2.011  0.0653
Enhancer_Hoffman.bed  1.886  1.15   2.307  0.4567
```
--cwd: the directory outputs are written to.

Output#

{cwd}/{index_snp_file}.gregor.conf - the GREGOR configuration file built from the supplied paths and parameters. The key lines of output/gregor/protocol_example.index.gregor.conf:

INDEX_SNP_FILE = input/ld/protocol_example.index.snps.txt
BED_FILE_INDEX = input/enrichment/protocol_example.bed.file.index
REF_DIR = input/enrichment/protocol_example.gregor_ref
R2THRESHOLD = 0.7 ## must be greater than 0.7
LDWINDOWSIZE = 10000 ## must be less than this window; these two values define LD buddies
POPULATION = EUR  ## define the population, you can specify EUR, AFR, AMR or ASN

{name}_gregor_output/StatisticSummaryFile.txt - GREGOR’s raw overlap statistics, one row per annotation bed file.
{name}_variant_counts.txt and {name}_enrichment_results.txt - the parsed 2x2 counts (inside versus outside the annotation peak, positive versus matched-negative set) and the per-annotation Fisher exact-test p-values and odds ratios.
{cwd}/{fisher1}_vs_{fisher2}_enrichment.pdf - the odds-ratio comparison plot across annotations.

Only the configuration file ships with the repository. Carrying the example past step 1 needs the full GREGOR reference database, so output/gregor holds the .conf and GREGOR’s generated job scripts but no summary, counts, or results tables.

Minimal Working Example#

Steps 1 and 2 need the GREGOR reference database and the GREGOR perl software; step 3 runs on the two example result tables alone.

GREGOR enrichment scan#

Configuration file only#

gregor_conf only writes a .gregor.conf from the supplied paths and parameters and never invokes GREGOR.

Timing: TBD (on toy dataset)

sos run pipeline/gregor.ipynb gregor_conf \
    --gregor_db input/enrichment/protocol_example.gregor_ref \
    --index_snp_file input/ld/protocol_example.index.snps.txt \
    --bed_file_index input/enrichment/protocol_example.bed.file.index \
    --pop EUR \
    --cwd output/gregor

The full chain#

gregor repeats the configuration step, then runs GREGOR and Fisher-tests its output. It expects the GREGOR perl tool as GREGOR on PATH and a reference database.

Timing: TBD (on toy dataset)

sos run pipeline/gregor.ipynb gregor \
    --gregor_db input/enrichment/protocol_example.gregor_ref \
    --index_snp_file input/ld/protocol_example.index.snps.txt \
    --bed_file_index input/enrichment/protocol_example.bed.file.index \
    --pop EUR \
    --cwd output/gregor

Plot an odds-ratio comparison#

gregor_fisher_plot is self-contained R and runs from two *_enrichment_results.txt tables alone, with no reference database or GREGOR install needed.

Timing: TBD (on toy dataset)

sos run pipeline/gregor.ipynb gregor_fisher_plot \
    --fisher1 input/enrichment/protocol_example.trait1_enrichment_results.txt \
    --fisher2 input/enrichment/protocol_example.trait2_enrichment_results.txt \
    --cwd output/gregor

Command Interface#

sos run pipeline/gregor.ipynb -h

usage: sos run pipeline/gregor.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:
  gregor_conf
  gregor
  gregor_fisher_plot

Global Workflow Options:
  --cwd output (as path)
                        working directory
  --container ''
                        Software container option

Sections
  gregor_conf, gregor_1: make configuration file for GREGOR
    Workflow Options:
      --gregor-db VAL (as path, required)
      --pop EUR
      --index-snp-file VAL (as path, required)
      --bed-file-index VAL (as path, required)
      --r2-threshold 0.7 (as float)
      --ld-window-size 10000 (as int)
      --min-neighbor 10 (as int)
      --job-number 10 (as int)
  gregor_2:             run GREGOR
  gregor_3:             Fisher test of enrichment
  gregor_fisher_plot:   enrichment plot with fisher test results
    Workflow Options:
      --fisher1 VAL (as path, required)
      --fisher2 VAL (as path, required)

Workflow implementation#

[global]
# working directory
parameter: cwd = path("output")
# Software container option
parameter: container = ""

# make configuration file for GREGOR
[gregor_conf, gregor_1]
parameter: gregor_db = path
parameter: pop = 'EUR'
parameter: index_snp_file = path
parameter: bed_file_index = path
parameter:  r2_threshold = 0.7
parameter: ld_window_size = 10000
parameter: min_neighbor = 10
parameter: job_number = 10
input: index_snp_file, bed_file_index
output: f'{cwd:a}/{_input[0]:bnn}.gregor.conf'
report: output = f'{_output}', expand = True
    ##############################################################################
    # CHIPSEQ ENRICHMENT CONFIGURATION FILE
    # This configuration file contains run-time configuration of
    # CHIP_SEQ ENRICHMENT
    ###############################################################################
    ## KEY ELEMENTS TO CONFIGURE : NEED TO MODIFY
    ###############################################################################
    INDEX_SNP_FILE = {_input[0]}
    BED_FILE_INDEX = {_input[1]} 
    REF_DIR = {gregor_db}
    R2THRESHOLD = {r2_threshold} ## must be greater than 0.7
    LDWINDOWSIZE = {ld_window_size} ## must be less than this window; these two values define LD buddies
    OUT_DIR = {_output:nn}_gregor_output
    MIN_NEIGHBOR_NUM = {min_neighbor} ## define the minimum size of neighborhood
    BEDFILE_IS_SORTED = true  ## false, if the bed files are not sorted
    POPULATION = {pop}  ## define the population, you can specify EUR, AFR, AMR or ASN
    TOPNBEDFILES = 2 
    JOBNUMBER = {job_number}
    ###############################################################################
    #BATCHTYPE = mosix ##  submit jobs on MOSIX
    #BATCHOPTS = -E/tmp -i -m2000 -j10,11,12,13,14,15,16,17,18,19,120,122,123,124,125 sh -c
    ###############################################################################
    #BATCHTYPE = slurm   ##  submit jobs on SLURM
    #BATCHOPTS = --partition=broadwl --account=pi-mstephens --time=0:30:0
    ###############################################################################
    BATCHTYPE = local ##  run jobs on local machine

bash: expand = True, stderr = f'{_output}.stderr', stdout = f'{_output}.stdout'
    sed -i '/^$/d' {_output}

GREGOR is written in perl. If you don’t use containers, some libraries are required before one can run GREGOR:

sudo apt-get install libdbi-perl libswitch-perl libdbd-sqlite3-perl

With the docker image:

cd GREGOR_folder
docker run -v "$PWD:/usr/src/myapp" -it custom-perl
perl script/GREGOR --conf example/mvsusie_annotation.conf
perl script/GREGOR --conf example/susie_annotation.conf

# run GREGOR
[gregor_2]
output: f'{_input:nn}_gregor_output/StatisticSummaryFile.txt'
bash: expand = True, container = container, stderr = f'{_output}.stderr', stdout = f'{_output}.stdout'
    GREGOR --conf {_input} && touch {_output}

# Fisher test of enrichment
[gregor_3]
output: f'{_input:ad}_variant_counts.txt', f'{_input:ad}_enrichment_results.txt'
bash: expand = '$[ ]', stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout'

    cd $[_input:ad]/  # The issue was that,{_input:ad}_variant_counts.txt is not in the same dir as the input  
    # Loop through each subdirectory
    for dir in */; do
        # Ensure that the directory is not empty
        if [[ -d "$dir" ]]; then
            # Calculate the total number of lines in neighbor.*.txt files
            total_lines=$(find "$dir" -name 'neighbor.*.txt' -exec cat {} + | wc -l)

            # Count the number of neighbor.*.txt files
            file_count=$(find "$dir" -name 'neighbor.*.txt' | wc -l)

            # Calculate the adjusted row number
            row_num=$((total_lines - file_count))

            # Check if PValue.txt exists
            if [[ -f "${dir}/PValue.txt" ]]; then
                # Append the row count to PValue.txt and save as PValue_new.txt
                (cat "${dir}/PValue.txt"; echo "N1_N2 = $row_num") > "${dir}/PValue_N1_N2.txt"
            fi
        fi
    done

    # Header for the new summary file
    echo -e "Bed_File\tInBed_Index_SNP\tExpectNum_of_InBed_SNP\tPValue\tN1_N2\tNp\tNp_Nn" > $[_output[0]]

    # Calculate Np and Np_Nn for the neighbor_SNP directory
    if [[ -d "neighbor_SNP" ]]; then
        if [[ -f "neighbor_SNP/index.snp.neighbors.txt" ]]; then
            Np=$(($(wc -l < "index_SNP/annotated.index.snp.txt") - 1))
        else
            Np=0
        fi
        Np_Nn_files=($(find "neighbor_SNP" -name 'neighbor.chr*txt'))
        Np_Nn=0
        for file in "${Np_Nn_files[@]}"; do
            Np_Nn=$(($Np_Nn + $(wc -l < "$file")))
        done
        Np_Nn=$(($Np_Nn - ${#Np_Nn_files[@]}))
    fi

    # Loop through each subdirectory
    for dir in */; do
        # Ensure that the directory is not empty
        if [[ -d "$dir" ]]; then
            # Check if PValue_new.txt exists
            if [[ -f "${dir}PValue_N1_N2.txt" ]]; then
                # Extract values from PValue_new.txt
                inBedIndexSNPNum=$(grep "inBedIndexSNPNum" "${dir}PValue_N1_N2.txt" | cut -d '=' -f2 | tr -d '[:space:]')
                expectedS=$(grep "expectedS" "${dir}PValue_N1_N2.txt" | cut -d '=' -f2 | tr -d '[:space:]')
                p3=$(grep "p3" "${dir}PValue_N1_N2.txt" | cut -d '=' -f2 | tr -d '[:space:]')
                N1_N2=$(grep "N1_N2" "${dir}PValue_N1_N2.txt" | cut -d '=' -f2 | tr -d '[:space:]')
                # Get the directory name as the Bed_File
                bedFile=$(basename "$dir")
                # Append the data to the summary file
                echo -e "$bedFile\t$inBedIndexSNPNum\t$expectedS\t$p3\t$N1_N2\t$Np\t$Np_Nn" >> $[_output[0]]
            fi
        fi
    done

R: expand = '${ }', stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout'
    res <- read.table(${_output[0]:r}, sep ='\t', header=T) 
    for(i in 1:nrow(res)){
        n1 = res$InBed_Index_SNP[i]
        n2 = res$N1_N2[i] - n1
        np = res$Np[i]
        nn = res$Np_Nn[i] - np

        # Construct the contingency matrix
        dat = matrix(c(n1, n2, np - n1, nn - n2), nrow = 2)

        # Perform Fisher's exact test
        test_res = fisher.test(dat, alternative = 'two.sided')

        # Store the results in 'res'
        res$odds[i] <- test_res$estimate
        res$low[i] <- test_res$conf.int[1]
        res$high[i] <- test_res$conf.int[2]
        res$p_fisher[i] <- test_res$p.value
    }

    res$odds <- as.numeric(res$odds)
    write.table(res, ${_output[1]:r}, quote = FALSE, row.names = FALSE)

# enrichment plot with fisher test results
[gregor_fisher_plot]
parameter: fisher1 = path
parameter: fisher2 = path
input: fisher1, fisher2
output: f'{cwd:a}/{_input[0]:dbn}_vs_{_input[1]:dbn}_enrichment.pdf'
R: expand = '${ }', stderr = f'{_output[0]}.stderr', stdout = f'{_output[0]}.stdout', container = container
    library(tidyverse)
    library(ggplot2)    
    group1 <- read.table(${_input[0]:r}, header=T)
    group2 <- read.table(${_input[1]:r}, header=T)
    group1$group <- ${_input[0]:dbnr}
    group2$group <- ${_input[1]:dbnr}
    res_all <- rbind(group1, group2)
    res_all$feature <- gsub(".bed","",res_all$Bed_File)
    res_all <- res_all %>%filter(!(str_detect(feature,"hg38")))%>%filter(!(p_fisher == 1))%>% na.omit 
    p <- res_all%>%
        arrange(odds)%>%ggplot()+geom_point(aes(x = odds, y = reorder(feature,-odds), color = group))+
            geom_vline(aes( xintercept = 0 ))+theme_bw()+theme(text = element_text(size = 20))+xlab("Odds Ratio")+
                geom_linerange(aes(xmin = low, xmax= high , y = feature ), color= 'grey40')+ylab("Functional Annotations")
    ggsave(plot = p, filename = ${_output:r}, height = 24, width = 20)

GREGOR enrichment analysis

Contents

GREGOR enrichment analysis#

Overview#

Input#

Output#

Minimal Working Example#

GREGOR enrichment scan#

Configuration file only#

The full chain#

Plot an odds-ratio comparison#

Command Interface#

Workflow implementation#