New computational tools for epigenome analysis

Outline

PEPATAC strengths

$ /pipelines/pepatac.py -h

usage: pepatac.py [-h] [-R] [-N] [-D] [-F] [-C CONFIG_FILE]
                  [-O PARENT_OUTPUT_FOLDER] [-M MEMORY_LIMIT]
                  [-P NUMBER_OF_CORES] -S SAMPLE_NAME -I INPUT_FILES
                  [INPUT_FILES ...] [-I2 [INPUT_FILES2 [INPUT_FILES2 ...]]] -G
                  GENOME_ASSEMBLY [-Q SINGLE_OR_PAIRED] [-gs GENOME_SIZE]
                  [--frip-ref-peaks FRIP_REF_PEAKS] [--TSS-name TSS_NAME]
                  [--anno-name ANNO_NAME] [--keep]
                  [--peak-caller {fseq,macs2}]
                  [--trimmer {trimmomatic,skewer}]
                  [--prealignments PREALIGNMENTS [PREALIGNMENTS ...]] [-V]

PEPATAC version 0.7.3

optional arguments:
  -h, --help            show this help message and exit
  -R, --recover         Overwrite locks to recover from previous failed run
  -N, --new-start       Overwrite all results to start a fresh run
  -D, --dirty           Don't auto-delete intermediate files
  -F, --force-follow    Always run 'follow' commands
  -C CONFIG_FILE, --config CONFIG_FILE
                        Pipeline configuration file (YAML). Relative paths are
                        with respect to the pipeline script.
  -O PARENT_OUTPUT_FOLDER, --output-parent PARENT_OUTPUT_FOLDER
                        Parent output directory of project
  -M MEMORY_LIMIT, --mem MEMORY_LIMIT
                        Memory limit (in Mb) for processes accepting such
  -P NUMBER_OF_CORES, --cores NUMBER_OF_CORES
                        Number of cores for parallelized processes
  -I2 [INPUT_FILES2 [INPUT_FILES2 ...]], --input2 [INPUT_FILES2 [INPUT_FILES2 ...]]
                        Secondary input files, such as read2
  -Q SINGLE_OR_PAIRED, --single-or-paired SINGLE_OR_PAIRED
                        Single- or paired-end sequencing protocol
  -gs GENOME_SIZE, --genome-size GENOME_SIZE
                        genome size for MACS2
  --frip-ref-peaks FRIP_REF_PEAKS
                        Reference peak set for calculating FRiP
  --TSS-name TSS_NAME   Name of TSS annotation file
  --anno-name ANNO_NAME
                        Name of reference bed file for calculating FRiF
  --keep                Keep prealignment BAM files
  --peak-caller {fseq,macs2}
                        Name of peak caller
  --trimmer {trimmomatic,pyadapt,skewer}
                        Name of read trimming program
  --prealignments PREALIGNMENTS [PREALIGNMENTS ...]
                        Space-delimited list of reference genomes to align to
                        before primary alignment.
  -V, --version         show program's version number and exit

required named arguments:
  -S SAMPLE_NAME, --sample-name SAMPLE_NAME
                        Name for sample to run
  -I INPUT_FILES [INPUT_FILES ...], --input INPUT_FILES [INPUT_FILES ...]
                        One or more primary input files
  -G GENOME_ASSEMBLY, --genome GENOME_ASSEMBLY
                        Identifier for genome assembly

PEP specification for sample metadata

pep_version: 2.0.0
sample_table: "path/to/sample_table.csv"

"sample_name", "protocol", "file"
"frog_1", "ATAC-seq", "frog1.fq.gz"
"frog_2", "ATAC-seq", "frog2.fq.gz"
"frog_3", "ATAC-seq", "frog3.fq.gz"
"frog_4", "ATAC-seq", "frog4.fq.gz"

MapReduce or Scatter/Gather

looper run config.yaml

looper runp config.yaml

PEPATAC strengths

PEPATAC strengths

Flexibility and Portability

# basic tools 
tools:  # absolute paths to required tools
  java: java
  python: python
  samtools: samtools
  bedtools: bedtools
  bowtie2: bowtie2
  fastqc: fastqc
  macs2: macs2
  picard: ${PICARD}
  skewer: skewer
  perl: perl
  # ucsc tools
  bedGraphToBigWig: bedGraphToBigWig
  wigToBigWig: wigToBigWig
  bigWigCat: bigWigCat
  bedSort: bedSort
  bedToBigBed: bedToBigBed
  # optional tools
  fseq: fseq  
  trimmo: ${TRIMMOMATIC}
  Rscript: Rscript 

# user configure 
resources:
  genomes: ${GENOMES}
  adapters: null  # Set to null to use default adapters

parameters:  # parameters passed to bioinformatic tools
  samtools:
    q: 10
  macs2: 
    f: BED
    q: 0.01
    shift: 0
  fseq:
    of: npf    # narrowPeak as output format
    l: 600     # feature length
    t: 4.0     # "threshold" (standard deviations)
    s: 1       # wiggle track step

PEPATAC strengths

Output

PEPATAC in practice

Problem

Who is the authoritative provider of the reference genome?

• NCBI?
• UCSC?
• Ensembl?

• hard, soft, or no repeat masking?
• are alternative scaffolds included?
• are haplotypes included?
• how are chromosomes named (chr1, 1, or NC_000001.11)?
• how is the assembly named (hg38, GRCh38, or GCF_000001405.39)?
• Are any decoy sequences included (like EBV)?

How refget works

Extending to sequence collections

• 1. An algorithm to create a deterministic, unique digest from a collection of sequences
• 2. A server capable of retrieving sequences given an identifier

First pass: Refgenie approach

Stolarczyk, Xue, and Sheffield (2021). NAR Genomics and Bioinformatics.

Limitations and discussion

• Should we include sequence topology in the digest?
• What other attributes could we include?
• Are there better delimiters?
• How do we construct the 'string-to-digest'?
• How do we handle order of sequences?
• How should the API respond to requests?

Project goal:

The project specifies:

How do we digest a sequence collection?

{
  "lengths": [
    4,
    4,
    8
  ],
  "names": [
    "chr1",
    "chr2",
    "chrX"
  ],
  "sequences": [
    "31fc6ca291a32fb9df82b85e5f077e31",
    "92c6a56c9e9459d8a42b96f7884710bc",
    "5f63cfaa3ef61f88c9635fb9d18ec945"
  ]
}

Digest algorithm

1. Canonicalize each attribute following RFC-8785 (JSON Canonicalization Scheme)
2. Digest each string (GA4GH digest: SHA512 truncated to 24 bits, converted to base64)
3. Canonicalize the entire object
4. Digest the canonicalized string

Advantages

• Accommodates new attributes with backwards-compatibility
• Additional layer of recursion to assess individual attributes
• Relies on existing JCS standard for string encoding

What gets digested?

Comparison function

• seqcol 1: 047c6e1eda552b50c5add59ff0995
• seqcol 2: 2230c535660fb4774114bfa966a62

How compatible are they?

Seqcol API demonstration

API endpoints

Conclusions

• Refget provides universal IDs for individual sequences
• Sequence collections extends this to reference genomes
• Using a deterministic algorithm, you can find the identifier
• A lookup service can retrieve the original sequence
• A comparison function allows fine-grained compatibility tests
• Please follow along: https://github.com/ga4gh/seqcol-spec

Word embeddings

Word2vec model

Word2vec model

Mikolov et al. (2013). arXiv:1301.3781v3.

Word context

Genomic context

Genomic Interval Embeddings

Evaluation

Evaluation 1: Classification performance

Evaluation 1: Classification performance

Evaluation 1: Classification performance

Evaluation 2: Perturbation similarity detection

Evaluation 3: Peak threshold robustness

Conclusion

• Regionset2vec uses an adapted word2vec model to train vectors for genomic regions
• Regionset2vec embeddings capture expected biological annotations
• Regionset2vec reflects known simulated perturbations
• Regionset2vec is robust to missing data
• NLP approaches can be adapted for applications in genomic interval analysis

Next steps

A high-performance server and API
for genomic interval data.

http://bedbase.org

• Data spans projects (e.g. all data on GEO; 40,000 accessions, 100,000+ BED files)
• Programmatic API for metadata, statistics, and data chunks
• Human browsing of statistical and biological attributes
• Aware of similarities among BED files
• Human-friendly search
• Shaped into 'non-redundant' sets for analysis

Thank You

Funding:



NIGMS R35-GM128636

nsheff

databio.org

nsheffield@virginia.edu