Repository folder:

Link to GitHub

Relevant reading to prepare yourself to the course

• Reinke, Aaron W, Raymond Mak, Emily R Troemel, and Eric J Ben. 2017. “In Vivo Mapping of Tissue- and Subcellular-Specific Proteomes in Caenorhabditis Elegans.” Science Advances 3 (5): 1–12.
• Jänes, Jürgen, Yan Dong, Michael Schoof, Jacques Serizay, Alex Appert, Chiara Cerrato, Carson Woodbury, et al. 2018. “Chromatin accessibility dynamics across C. elegans development and ageing.” eLife, October. https://doi.org/10.7554/eLife.37344.
• Bioconductor (Huber et al. 2015) and GenomicRanges (Marc Carlson and Morgan 2019) documentations

Required sowftware

• R (3.5.2)
• org.Hs.eg.db_3.7.0
• clusterProfiler_3.10.1
• htmlwidgets_1.3
• rtracklayer_1.42.2
• ggplot2_3.1.1
• DOSE_3.8.2
• tidyverse_1.2.1
• gprofiler2_0.1.3
• biomaRt_2.38.0

Dependencies (which should have been installed with the main required packages)

• S4Vectors_0.20.1
• stringr_1.4.0
• readr_1.3.1
• AnnotationDbi_1.44.0
• GenomeInfoDb_1.18.2
• BiocGenerics_0.28.0
• dplyr_0.8.0.1
• tidyr_0.8.3
• Biobase_2.42.0
• IRanges_2.16.0
• forcats_0.4.0
• purrr_0.3.2
• tibble_2.1.1

• Prerequisites
  • Relevant reading to prepare yourself to the course
  • Required sowftware
  • Dependencies (which should have been installed with the main required packages)
• Day 1: Introduction
  • Research in aging
  • Aging and the worm
  • The question
  • The approach
  • About this vignette
• Day 2: Proteomics
  • Set up R environment
  • Load and inspect the data
  • Adding gene information
  • Separating into lists of tissue-specific proteins
  • Compare it to existing tissue-specific datasets from RNA-seq
  • Save your work!!
• Day 3: ATAC-seq
  • Set up R environment and load previous work
  • Investigate the ATAC-seq aging data
  • Investigate the genes associated to varying promoters during aging
  • [EXTRA] Using external visualization tools
  • Save your work!!
• Day 4: ChIP-seq
  • [EXTRA] Get the ChIP-seq data by yourself from modENCODE platform
  • Set up R environment and load previous work
  • Investigate the TFs associated to varying promoters during aging
  • Defining key TFs involved in aging?
  • [EXTRA] Finding motifs associated to each TFs
  • Save your work!!
• SessionInfo
• References

Research in aging

Aging is a time-dependent functional decline that affects most living organisms. It is characterized by a progressive loss of physiological integrity, leading to impaired function and increased vulnerability to death.

The accumulation of cellular damage is widely considered to be the general cause of aging (López-Otı́n et al. 2013), which explains why the studies from the last three decades have mostly focused on the cellular aspects of aging.
Currently, the aging hallmarks are:

• Genomic instability;
• Telomere attrition;
• Epigenetic alterations;
• Loss of proteostasis;
• Deregulated nutrient sensing;
• Mitochondrial dysfunction;
• Cellular senescence;
• Stem cell exhaustion;
• Altered intercellular communication.

However, aging should not only be seen as a cellular disfunction but as an organismal process. All the cells within an organism are undergoing aging, and not all of them are aging simultaneously. Overall, different tissues are aging differently (Son et al. 2019). Thus, aging is a decline affecting not only cellular components but also organismal biological processes in general.

Aging and the worm

Aging research was inaugurated more than 30 years ago, following the isolation of the first long-lived strains in C. elegans (Johnson 2013). In the past three decades, the nematode played a fundamental role in aging research. It is a powerful model organism when it comes to studying aging:

• It has a short life cycle (~ 21 days in normal conditions);
• Mutants can be easily generated and phenotypically screened for aging defects;
• Longevity assays are easy to perform and very informative;
• Its genome has been fully sequenced and presents fundamental resemblances with that of more complex metazoans.

One of the most used phenotypic assays in worm is the longevity assay. Longevity assays are useful to understand whether a gene is involved in aging (Chen et al. 2015). RNAi is used to knock-down a gene of interest and the resulting longevity is measured. If the longevity increases, the gene function normally reduce it.
In worms, RNAi can be used in large-scale screens. RNAi screens have been performed and have led to the identification of tens of genes involved in aging and regulation of longevity in general. It has been instrumental to decipher the molecular pathways involved in aging (Hamilton et al. 2005).

The question

If RNAi has been instrumental to generally study the process of aging in worm, the readouts are generally phenotypic readouts, e.g. does the worm live longer? Or does it rapidly stop moving after knock-down? These readouts are poorly suited to understand the aging-related role of a protein coded by a given gene in specific tissues. For instance, slo-1 is a gene coding for a potassium channel. Its mutation or knock-down accelerates aging (Chen et al. 2013). Because we have some information about the nature and function of the coded protein, it is rather easy to speculate in each tissue it has an aging-related function (hint: the neurons!).
However, if we do not know anything about the gene (e.g. in which tissue it is expressed, which transcription factors regulate its expression, the type of protein encoded by it, …), it becomes challenging to draw conclusions on the molecular mechanisms involving the protein encoded by the gene.
Rather than relying on RNAi screens, one could leverage the large amount of “omics” data generated in C. elegans to better characterize aging processes at the organismal level. The goal of this study is to analyze aging-related perturbations in worm in a comprehensive manner to (1) identify the tissues undergoing important remodeling during aging and (2) understand the molecular networks involved in aging. More specifically, we would like to fulfill the following aims:

1. Find tissue-specific proteins and genes deregulated during aging;
2. Find their associations with known diseases;
3. Find potential targets for drug engineering

The approach

RNAi has been extensively used to study aging in worm. However, next-generation techniques such as high-throughput sequencing or mass spectrometry now allow for more comprehensive multiomics studies of biological processes. We will rely on results from these techniques to:

1. Day 2: Annotate sets of tissue-specific proteins (using proteomics results) and compare it with sets of tissue-specific genes (using transcriptomics results);
2. Day 3: Find genomic regions regulating these genes (using epigenomics results);
3. Day 4: Identify networks of transcription factors regulating these genes (using epigenomics results).

About this project (and the vignette)

Over the next three days, we will be conducting analysis of several types of datasets in R. This vignette is to be used as a guideline and lists important questions or technical challenges. The R code is only here to suggest one solution to the question(s) that are raised. Before looking at it, it is recommended to investigate the data by yourself and try out your own commands.
More generally, you are welcome to investigate things your own way! If you think of relevant questions that are not raised in this vignette, which you could try and answer using the provided data, don’t hesitate to go off the beaten path. Mentors will be around to supervise everybody - even if the investigations deviate from the original direction - and any personal initiative is welcome and encouraged! Finally, it is important to stress out that everyone is not expected to investigate all the points covered in this vignette and to fully complete every single task as instructed. Again, this vignette is more of a general guideline; if you feel like you don’t want to dedicate time solving a particular point (e.g. you already know everything about investigating pathologies-related lists of genes), feel free to skip parts.

Set up R environment

You will first need to make sure that you are working in the right folder. You can also load important packages at this point.

    PROJECT_PATH <- getwd()
    require(tidyverse)
## Loading required package: tidyverse
## ── Attaching packages ────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.1       ✔ purrr   0.3.2  
## ✔ tibble  2.1.1       ✔ dplyr   0.8.0.1
## ✔ tidyr   0.8.3       ✔ stringr 1.4.0  
## ✔ readr   1.3.1       ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
    require(rtracklayer)
## Loading required package: rtracklayer
## Loading required package: GenomicRanges
## Loading required package: stats4
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind,
##     colMeans, colnames, colSums, dirname, do.call, duplicated,
##     eval, evalq, Filter, Find, get, grep, grepl, intersect,
##     is.unsorted, lapply, lengths, Map, mapply, match, mget, order,
##     paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind,
##     Reduce, rowMeans, rownames, rowSums, sapply, setdiff, sort,
##     table, tapply, union, unique, unsplit, which, which.max,
##     which.min
## Loading required package: S4Vectors
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## The following object is masked from 'package:tidyr':
## 
##     expand
## The following object is masked from 'package:base':
## 
##     expand.grid
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## The following object is masked from 'package:purrr':
## 
##     reduce
## Loading required package: GenomeInfoDb
    require(biomaRt)
## Loading required package: biomaRt

Load and inspect the data

The data has been fetched from the original paper and reformatted for the purpose of the course. Start by loading it from data/1602426_TableS4.csv.

    proteomics <- read.csv(
        file.path(PROJECT_PATH, '../data/1602426_TableS4_processed.csv'), 
        na.strings = "NA"
    )
    # Number of proteins in the file
    nrow(proteomics)
## [1] 3722
    # Number of proteins detected in each tissue
    colSums(proteomics[,38:45])
##         Epidermis_cytoplasm           Epidermis_nucleus 
##                        2442                        1278 
## Pharyngeal.muscle_cytoplasm   Pharyngeal.muscle_nucleus 
##                         837                         104 
##  Body.wall.muscle_cytoplasm    Body.wall.muscle_nucleus 
##                        1375                         609 
##         Intestine_cytoplasm           Intestine_nucleus 
##                        2449                         971

Adding gene information

This file only contains the protein names. Eventually, we will be looking at genes as well. The simplest approach is to convert everything to unique gene IDs. For this, we will rely on BiomaRt.

    ### You can get IDs of genes and associated proteins using biomaRt
    # ensembl <- useDataset("celegans_gene_ensembl", mart = useMart("ensembl"))
    # ids <- getBM(
    #     attributes = c("uniprotswissprot", "uniprotsptrembl", "wormbase_gene"),
    #     mart = ensembl
    # )
    # Otherwise, you can also get the ids table from the data folder: 
    ids <- read.table('../data/proteins_genes_IDs.txt', header = TRUE, sep = '\t')
    # Merging SwissProt and Trembl IDs into one column will make your life easier
    ids$protID <- paste0(ids$uniprotswissprot, ids$uniprotsptrembl)
    # Then a simple match function can be used to translate the protein IDs into gene IDs
    proteomics$WormBaseID <- ids$wormbase_gene[match(proteomics$UniProtKB, ids$protID)]
    # Few proteins do not have any gene ID (why?). We will ignore these ones for the rest of the study.
    proteomics <- proteomics[!is.na(proteomics$WormBaseID),]

Separating into lists of tissue-specific proteins

Let’s figure out which proteins are specifically present in each tissue.

    # Let's retrieve the list of proteins specifically enriched in a single tissue
    list_prots <- lapply(
        levels(proteomics$Tissue.specific), 
        function(TISSUE) {
            proteomics$UniProtKB %>%
                '['(proteomics$Tissue.specific == TISSUE & !is.na(proteomics$Tissue.specific)) %>% 
                as.character()
        }
    ) %>% setNames(levels(proteomics$Tissue.specific))
    lengths(list_prots)
##  Body wall muscle specific         Epidermis specific 
##                         55                        125 
##         Intestine specific Pharyngeal muscle specific 
##                        130                         21
    # Because we will need it later, you can also get the same list but with gene IDs
    list_genes <- lapply(
        levels(proteomics$Tissue.specific), 
        function(TISSUE) {
            proteomics$WormBaseID %>%
                '['(proteomics$Tissue.specific == TISSUE & !is.na(proteomics$Tissue.specific)) %>% 
                as.character()
        }
    ) %>% setNames(levels(proteomics$Tissue.specific))
    # A list of most of the transcription factors found in C. elegans genome can be found in data/TFs.txt
    tfs <- readLines('../data/TFs.txt')
    length(tfs)
## [1] 824
    # Let's see which proteins enriched in tissues are transcription factors
    tfs_in_tissue_specific_prots <- lapply(list_prots, function(prots) {
        prots[proteomics$WormBaseID[match(prots, proteomics$UniProtKB)] %in% tfs]
    })
    # Let's see if transcription factors are enriched 
    n_prots <- nrow(ids)
    n_tfs <- length(tfs)
    n_tissue_spe_prots <- sum(lengths(list_prots))
    n_tfs_in_tissue_specific_prots <- sum(lengths(tfs_in_tissue_specific_prots))
    contingency_matrix <- matrix(
        cbind(
            c(n_tfs_in_tissue_specific_prots, n_tfs - n_tfs_in_tissue_specific_prots), 
            c(n_tissue_spe_prots - n_tfs_in_tissue_specific_prots, n_prots - n_tfs - n_tissue_spe_prots + n_tfs_in_tissue_specific_prots)
        ), 
        nrow = 2
    )
    fisher.test(contingency_matrix)
## 
##  Fisher's Exact Test for Count Data
## 
## data:  contingency_matrix
## p-value = 0.6
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.3354 1.5753
## sample estimates:
## odds ratio 
##     0.7862

Compare it to existing tissue-specific datasets from RNA-seq

Let’s see if this set of proteins overlap with tissue-specific gene annotation. For this, we can use the data in data/gene_annotations.gff3. This is a .gff3 file, a format used to add specific information to gene annotations. This annotation file has not been published yet.

    # Let's import the gene annotation information
    genes <- import('../data/gene_annotations.gff3')
    names(genes) <- genes$ID
    # It will be easier if we convert the tissue-specific gene annotation to a vector
    genes$which.tissues <- factor(genes$which.tissues, levels = c(
        'Germline', 'Sperm', 'EarlyGermline', 'Neurons', 'Muscle', 'Hypod.', 'Intest.', 
        'Germline_Neurons', 'Germline_Muscle', 'Germline_Hypod.', 'Germline_Intest.', 'Neurons_Muscle', 'Neurons_Hypod.', 'Neurons_Intest.', 'Muscle_Hypod.', 'Muscle_Intest.', 'Hypod._Intest.',
        'Germline_Neurons_Muscle', 'Germline_Neurons_Hypod.', 'Germline_Neurons_Intest.', 'Germline_Muscle_Hypod.', 'Germline_Muscle_Intest.', 'Germline_Hypod._Intest.', 'Neurons_Muscle_Hypod.', 'Neurons_Muscle_Intest.', 'Neurons_Hypod._Intest.', 'Muscle_Hypod._Intest.', 
        'Germline_Neurons_Muscle_Hypod.', 'Germline_Neurons_Muscle_Intest.', 'Germline_Neurons_Hypod._Intest.', 'Germline_Muscle_Hypod._Intest.', 'Neurons_Muscle_Hypod._Intest.', 
        'Germline_Neurons_Muscle_Hypod._Intest.', 
        'Soma', 'Ubiq.', 'Ubiq.-Reg', 'Unclassified', 'Low', 'non-prot-cod'
    ))
    # We can compare the two different tissue-specific annotations
    lapply(list_genes, function(g) {
        table(genes[g]$which.tissues)
    })
## $`Body wall muscle specific`
## 
##                               Germline 
##                                      0 
##                                  Sperm 
##                                      0 
##                          EarlyGermline 
##                                      0 
##                                Neurons 
##                                      0 
##                                 Muscle 
##                                     34 
##                                 Hypod. 
##                                      0 
##                                Intest. 
##                                      0 
##                       Germline_Neurons 
##                                      0 
##                        Germline_Muscle 
##                                      0 
##                        Germline_Hypod. 
##                                      0 
##                       Germline_Intest. 
##                                      0 
##                         Neurons_Muscle 
##                                      3 
##                         Neurons_Hypod. 
##                                      0 
##                        Neurons_Intest. 
##                                      0 
##                          Muscle_Hypod. 
##                                      0 
##                         Muscle_Intest. 
##                                      0 
##                         Hypod._Intest. 
##                                      0 
##                Germline_Neurons_Muscle 
##                                      0 
##                Germline_Neurons_Hypod. 
##                                      0 
##               Germline_Neurons_Intest. 
##                                      0 
##                 Germline_Muscle_Hypod. 
##                                      0 
##                Germline_Muscle_Intest. 
##                                      0 
##                Germline_Hypod._Intest. 
##                                      0 
##                  Neurons_Muscle_Hypod. 
##                                      0 
##                 Neurons_Muscle_Intest. 
##                                      2 
##                 Neurons_Hypod._Intest. 
##                                      0 
##                  Muscle_Hypod._Intest. 
##                                      0 
##         Germline_Neurons_Muscle_Hypod. 
##                                      0 
##        Germline_Neurons_Muscle_Intest. 
##                                      0 
##        Germline_Neurons_Hypod._Intest. 
##                                      0 
##         Germline_Muscle_Hypod._Intest. 
##                                      0 
##          Neurons_Muscle_Hypod._Intest. 
##                                      0 
## Germline_Neurons_Muscle_Hypod._Intest. 
##                                      0 
##                                   Soma 
##                                      3 
##                                  Ubiq. 
##                                      4 
##                              Ubiq.-Reg 
##                                      0 
##                           Unclassified 
##                                      4 
##                                    Low 
##                                      0 
##                           non-prot-cod 
##                                      0 
## 
## $`Epidermis specific`
## 
##                               Germline 
##                                      0 
##                                  Sperm 
##                                      0 
##                          EarlyGermline 
##                                      0 
##                                Neurons 
##                                      0 
##                                 Muscle 
##                                      0 
##                                 Hypod. 
##                                     41 
##                                Intest. 
##                                      1 
##                       Germline_Neurons 
##                                      0 
##                        Germline_Muscle 
##                                      0 
##                        Germline_Hypod. 
##                                      1 
##                       Germline_Intest. 
##                                      0 
##                         Neurons_Muscle 
##                                      0 
##                         Neurons_Hypod. 
##                                      0 
##                        Neurons_Intest. 
##                                      1 
##                          Muscle_Hypod. 
##                                      0 
##                         Muscle_Intest. 
##                                      0 
##                         Hypod._Intest. 
##                                      0 
##                Germline_Neurons_Muscle 
##                                      0 
##                Germline_Neurons_Hypod. 
##                                      0 
##               Germline_Neurons_Intest. 
##                                      0 
##                 Germline_Muscle_Hypod. 
##                                      0 
##                Germline_Muscle_Intest. 
##                                      0 
##                Germline_Hypod._Intest. 
##                                      0 
##                  Neurons_Muscle_Hypod. 
##                                      1 
##                 Neurons_Muscle_Intest. 
##                                      0 
##                 Neurons_Hypod._Intest. 
##                                      0 
##                  Muscle_Hypod._Intest. 
##                                      0 
##         Germline_Neurons_Muscle_Hypod. 
##                                      5 
##        Germline_Neurons_Muscle_Intest. 
##                                      0 
##        Germline_Neurons_Hypod._Intest. 
##                                      0 
##         Germline_Muscle_Hypod._Intest. 
##                                      0 
##          Neurons_Muscle_Hypod._Intest. 
##                                      0 
## Germline_Neurons_Muscle_Hypod._Intest. 
##                                      0 
##                                   Soma 
##                                      7 
##                                  Ubiq. 
##                                     19 
##                              Ubiq.-Reg 
##                                      0 
##                           Unclassified 
##                                     10 
##                                    Low 
##                                      4 
##                           non-prot-cod 
##                                      0 
## 
## $`Intestine specific`
## 
##                               Germline 
##                                      2 
##                                  Sperm 
##                                      0 
##                          EarlyGermline 
##                                      0 
##                                Neurons 
##                                      1 
##                                 Muscle 
##                                      0 
##                                 Hypod. 
##                                      2 
##                                Intest. 
##                                     39 
##                       Germline_Neurons 
##                                      0 
##                        Germline_Muscle 
##                                      0 
##                        Germline_Hypod. 
##                                      0 
##                       Germline_Intest. 
##                                      0 
##                         Neurons_Muscle 
##                                      0 
##                         Neurons_Hypod. 
##                                      0 
##                        Neurons_Intest. 
##                                      0 
##                          Muscle_Hypod. 
##                                      0 
##                         Muscle_Intest. 
##                                      2 
##                         Hypod._Intest. 
##                                      2 
##                Germline_Neurons_Muscle 
##                                      0 
##                Germline_Neurons_Hypod. 
##                                      0 
##               Germline_Neurons_Intest. 
##                                      0 
##                 Germline_Muscle_Hypod. 
##                                      0 
##                Germline_Muscle_Intest. 
##                                      0 
##                Germline_Hypod._Intest. 
##                                      0 
##                  Neurons_Muscle_Hypod. 
##                                      0 
##                 Neurons_Muscle_Intest. 
##                                     11 
##                 Neurons_Hypod._Intest. 
##                                      0 
##                  Muscle_Hypod._Intest. 
##                                      1 
##         Germline_Neurons_Muscle_Hypod. 
##                                      0 
##        Germline_Neurons_Muscle_Intest. 
##                                      0 
##        Germline_Neurons_Hypod._Intest. 
##                                      0 
##         Germline_Muscle_Hypod._Intest. 
##                                      0 
##          Neurons_Muscle_Hypod._Intest. 
##                                      0 
## Germline_Neurons_Muscle_Hypod._Intest. 
##                                      0 
##                                   Soma 
##                                      5 
##                                  Ubiq. 
##                                     15 
##                              Ubiq.-Reg 
##                                      0 
##                           Unclassified 
##                                     18 
##                                    Low 
##                                      3 
##                           non-prot-cod 
##                                      0 
## 
## $`Pharyngeal muscle specific`
## 
##                               Germline 
##                                      0 
##                                  Sperm 
##                                      0 
##                          EarlyGermline 
##                                      0 
##                                Neurons 
##                                      0 
##                                 Muscle 
##                                      2 
##                                 Hypod. 
##                                      0 
##                                Intest. 
##                                      0 
##                       Germline_Neurons 
##                                      0 
##                        Germline_Muscle 
##                                      0 
##                        Germline_Hypod. 
##                                      0 
##                       Germline_Intest. 
##                                      0 
##                         Neurons_Muscle 
##                                      1 
##                         Neurons_Hypod. 
##                                      0 
##                        Neurons_Intest. 
##                                      0 
##                          Muscle_Hypod. 
##                                      0 
##                         Muscle_Intest. 
##                                      0 
##                         Hypod._Intest. 
##                                      0 
##                Germline_Neurons_Muscle 
##                                      1 
##                Germline_Neurons_Hypod. 
##                                      0 
##               Germline_Neurons_Intest. 
##                                      0 
##                 Germline_Muscle_Hypod. 
##                                      0 
##                Germline_Muscle_Intest. 
##                                      0 
##                Germline_Hypod._Intest. 
##                                      0 
##                  Neurons_Muscle_Hypod. 
##                                      0 
##                 Neurons_Muscle_Intest. 
##                                      3 
##                 Neurons_Hypod._Intest. 
##                                      0 
##                  Muscle_Hypod._Intest. 
##                                      0 
##         Germline_Neurons_Muscle_Hypod. 
##                                      0 
##        Germline_Neurons_Muscle_Intest. 
##                                      0 
##        Germline_Neurons_Hypod._Intest. 
##                                      0 
##         Germline_Muscle_Hypod._Intest. 
##                                      0 
##          Neurons_Muscle_Hypod._Intest. 
##                                      0 
## Germline_Neurons_Muscle_Hypod._Intest. 
##                                      0 
##                                   Soma 
##                                      1 
##                                  Ubiq. 
##                                      0 
##                              Ubiq.-Reg 
##                                      0 
##                           Unclassified 
##                                      4 
##                                    Low 
##                                      7 
##                           non-prot-cod 
##                                      0
    # We can also plot the two different tissue-specific annotations as a heatmap
    df <- data.frame(
        transcriptome_annotations = genes[unlist(list_genes)]$which.tissues, 
        proteomics_annotations = proteomics[match(unlist(list_genes), proteomics$WormBaseID),]$Tissue.specific
    ) %>% 
        table() %>% 
        as.data.frame() 
    p <- ggplot(df, aes(y = 1, x = transcriptome_annotations)) +
        geom_tile(aes(fill = Freq)) + 
        geom_text(aes(label = Freq)) +
        scale_fill_gradientn(colours = c('white', 'orange', 'red')) + 
        theme_bw() + 
        facet_wrap(~proteomics_annotations) + 
        labs(y = '', x = 'Transcriptomics-based annotations') + 
        theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank()) + 
        theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
    ggsave('../results/comparison-transcriptomics-proteomics.pdf', width = 15, height = 4)

Save your work!!

That’s it for today! Don’t forget to save your progress.

    save.image('../results/Day2.RData')

Set up R environment and load previous work

Once again, first make sure that you are working in the right folder. You can also load important packages at this point.

    PROJECT_PATH <- getwd()
    require(tidyverse)
    require(rtracklayer)
    require(htmlwidgets)
## Loading required package: htmlwidgets
    require(gprofiler2)
## Loading required package: gprofiler2
    require(DOSE)
## Loading required package: DOSE
## 
## DOSE v3.8.2  For help: https://guangchuangyu.github.io/DOSE
## 
## If you use DOSE in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609
    require(clusterProfiler)
## Loading required package: clusterProfiler
## clusterProfiler v3.10.1  For help: https://guangchuangyu.github.io/software/clusterProfiler
## 
## If you use clusterProfiler in published research, please cite:
## Guangchuang Yu, Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.
## 
## Attaching package: 'clusterProfiler'
## The following object is masked from 'package:purrr':
## 
##     simplify
    require(org.Hs.eg.db)
## Loading required package: org.Hs.eg.db
## Loading required package: AnnotationDbi
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## Attaching package: 'AnnotationDbi'
## The following object is masked from 'package:dplyr':
## 
##     select
## 
    # The work done in Day 2 has been stored in an RData object.
    load('../results/Day2.RData')

Investigate the ATAC-seq aging data

The data has been fetched from the original paper and reformatted for the purpose of the course. Start by loading it from data/ATAC-seq_aging.gff3.

    REs <- import('../data/ATAC-seq_aging.gff3')
    proms <- REs[REs$ageing_prom_cluster_label != '.']
    ATAC_aging <- mcols(proms)[, 8:12] %>% 
        data.frame() %>% 
        data.matrix() %>% 
        as.data.frame() 
    # How many promoters does each cluster contain?
    table(proms$ageing_prom_cluster_label)
## 
##       I I+H [1] I+H [2]    Mix1    Mix2    Mix3    Mix4    Mix5  stable 
##     304     252     409     210      98     273     112     142   11796
    # Comparison of accessibility changes across clusters
    norm_changes_df <- apply(ATAC_aging, 1, function(ROW) {log2(ROW/mean(ROW))}) %>% 
        t() %>% 
        data.frame()
    df <- data.frame(
        locus = proms$locus,
        norm_changes_df,
        cluster = proms$ageing_prom_cluster_label
    ) %>% 
        gather(stage, expr, -cluster, -locus) %>%
        mutate(stage = factor(stage, levels = colnames(norm_changes_df))) %>%
        filter(cluster != '.')
    p <- ggplot(df, aes(x = stage, y = expr, color = cluster, group = locus)) + 
        geom_line(alpha = 0.05) + 
        theme_bw() + 
        facet_wrap(~cluster) +
        theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank()) + 
        theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) + 
        theme(legend.position = 'none')
    ggsave('../results/accessibility-changes_clusters.pdf', height = 6, width = 6.4)
## Warning: Removed 30 rows containing missing values (geom_path).

    # Check the association between proms and proteins
    list_genes_aging <- lapply(levels(df$cluster), function(cluster) {
        proms <- proms[proms$ageing_prom_cluster_label == cluster & proms$regulatory_class %in% c('bidirect-promoter', 'fwd-promoter', 'rev-promoter')]
        genes_cluster <- proms$WormBaseID %>% as.character() %>% strsplit(',') %>% unlist()
        return(genes_cluster)
    }) %>% setNames(levels(df$cluster))
    # We can also plot something
    df <- data.frame(
        cluster_annotations = unlist(lapply(seq_along(list_genes_aging), function(K) {rep(names(list_genes_aging)[K], length(list_genes_aging[[K]]))})),
        proteomics_annotations = proteomics[match(unlist(list_genes_aging), proteomics$WormBaseID),]$Tissue.specific
    ) %>% 
        table() %>% 
        as.data.frame() %>% 
        filter(cluster_annotations != 'stable')
    p <- ggplot(df, aes(y = 1, x = cluster_annotations)) +
        geom_tile(aes(fill = Freq)) + 
        geom_text(aes(label = Freq)) +
        scale_fill_gradientn(colours = c('white', 'orange', 'red')) + 
        theme_bw() + 
        facet_wrap(~proteomics_annotations) + 
        labs(y = '', x = 'Ageing cluster') + 
        theme(axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank()) + 
        theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
    ggsave('../results/comparison-ATAC-clusters-proteomics.pdf', width = 8, height = 4)

Investigate the genes associated to varying promoters during aging

Now that we know which genomic loci are varying during aging, we can retrieve the associated genes and investigate their function and their associated with diseases.

    # What is the function of genes going up (down) during aging?
    genes_down <- unique(unlist(list_genes_aging[1:4]))
    genes_up <- unique(unlist(list_genes_aging[5:8]))
    gos <- gost(
        list(
            up = genes_up, 
            down = genes_down
        ),
        organism = 'celegans', 
        multi_query = TRUE, 
        user_threshold = 0.01, 
        source = c('GO:BP', 'GO:MF', 'GO:CC')
    )
    p <- gostplot(gos, capped = TRUE, interactive = TRUE)

    # What are the homologs of these genes in mammals? 
    ensembl <- useDataset("celegans_gene_ensembl", mart = useMart("ensembl"))
    #
    human_homologues_down <- getBM(
        filters = c("ensembl_gene_id"),
        values = genes_down, 
        attributes = c('ensembl_gene_id', 'hsapiens_homolog_ensembl_gene'),
        mart = ensembl
    )
## Batch submitting query [================>---------] 67% eta: 6s Batch
## submitting query [==========================] 100% eta: 0s
    human_homologues_down_ENTREZ <- human_homologues_down$hsapiens_homolog_ensembl_gene %>%
        '['(nchar(.) > 0) %>%
        bitr(fromType = "ENSEMBL", toType = c("ENTREZID"), OrgDb = org.Hs.eg.db) %>%
        '['(,"ENTREZID")
## 'select()' returned 1:many mapping between keys and columns
## Warning in bitr(., fromType = "ENSEMBL", toType = c("ENTREZID"), OrgDb =
## org.Hs.eg.db): 1.75% of input gene IDs are fail to map...
    #
    human_homologues_up <- getBM(
        filters = c("ensembl_gene_id"),
        values = genes_up, 
        attributes = c('ensembl_gene_id', 'hsapiens_homolog_ensembl_gene'),
        mart = ensembl
    )
## Batch submitting query [==========================] 100% eta: 0s
    human_homologues_up_ENTREZ <- human_homologues_up$hsapiens_homolog_ensembl_gene %>%
        '['(nchar(.) > 0) %>%
        bitr(fromType = "ENSEMBL", toType = c("ENTREZID"), OrgDb = org.Hs.eg.db) %>%
        '['(,"ENTREZID")
## 'select()' returned 1:many mapping between keys and columns
## Warning in bitr(., fromType = "ENSEMBL", toType = c("ENTREZID"), OrgDb =
## org.Hs.eg.db): 1.96% of input gene IDs are fail to map...

    diseases_enrichment_down <- enrichDO(
        gene = human_homologues_down_ENTREZ, 
        ont = "DO", 
        pvalueCutoff = 0.01,
        pAdjustMethod = "BH",
        readable = TRUE
    )
    head(diseases_enrichment_down@result[,1:7])
##                    ID                       Description GeneRatio  BgRatio
## DOID:2978   DOID:2978   carbohydrate metabolic disorder     9/401  27/8007
## DOID:655     DOID:655      inherited metabolic disorder    31/401 331/8007
## DOID:589     DOID:589       congenital hemolytic anemia     5/401  16/8007
## DOID:9993   DOID:9993                      hypoglycemia     4/401  10/8007
## DOID:14504 DOID:14504              Niemann-Pick disease     4/401  12/8007
## DOID:4330   DOID:4330 non-langerhans-cell histiocytosis     4/401  14/8007
##               pvalue p.adjust   qvalue
## DOID:2978  3.794e-06 0.002443 0.002443
## DOID:655   5.323e-04 0.164617 0.164617
## DOID:589   8.479e-04 0.164617 0.164617
## DOID:9993  1.022e-03 0.164617 0.164617
## DOID:14504 2.225e-03 0.286527 0.286527
## DOID:4330  4.154e-03 0.432135 0.432135

[EXTRA] Using external visualization tools

Data visualization softwares are key to investigate. If you have some spare time, you can try to use Seqplots software (Stempor and Ahringer 2016) to get more insights about the promoters clustered by their accessibility changes. The software requires:

• A genome to be installed (ce11, the latest genome version of C. elegans, should already be installed);
• Track files (typically bigWig files): for instance, genome-wide accessibility tracks generated during aging (see the data/ folder…)
• Some feature tracks (typically bed files): these files generally contain sets of genomic loci. For instance, you could generate 9 different bed files, each containing the promoters from a single cluster.

    # Export bed files of promoters from each cluster to visualize the accessibility profiles over each set using SeqPlots
    for (cluster in levels(df$cluster)) {
        export(granges(proms[proms$ageing_prom_cluster_label == cluster]), paste0('ATAC_cluster-', cluster, '.bed'))
    }
    seqplots::run()

Once the bed files are generated, you can load the tracks and the bed files in Seqplots and try the different functionalities of the software.

Save your work!!

Don’t forget to save your progress!

    save.image('../results/Day3.RData')

[EXTRA] Get the ChIP-seq data by yourself from modENCODE platform

modENCODE/modERN data is accessible on the ENCODE project web portal, which hosts many more datasets from different model systems. Raw data and processed files can be downloaded directly, or a list of files can be retrieved and downloaded programatically. If you have some spare time, you can try to download the original processed files from there and start to wrangle them by yourself.

Set up R environment and load previous work

Once again, first make sure that you are working in the right folder. You can also load important packages at this point.

    PROJECT_PATH <- getwd()
    require(tidyverse)
    require(rtracklayer)
    # The work done in Day 3 has been stored in an RData object.
    load('../results/Day3.RData')

Investigate the TFs associated to varying promoters during aging

Because getting the right list of files, tyding the downloaded files, QC-ing and removing bad quality datasets, merging replicates, etc… always require a lot of data wrangling, a clean set of .bed files (one per transcription factor) is provided in the data/ folder to avoid these cumbersome steps. Each bed file contains the transcription factor ChIP-seq “peaks”, corresponding to the genomic loci bound by the transcription factor. Some TFs are binding to thousands of loci, some others are bound to a few dozens of them. Moreover, the quality of the ChIP-seq datasets can vary between samples.

    # Read all the TF CHIP-seq bed files
    TFs <- list.files('../data/bed/') %>% gsub('.bed', '', .)
    TFs_peaks <- lapply(TFs, function(TF) {
        import(paste0('../data/bed/', TF, '.bed'))
    }) %>% setNames(TFs)
    # Compile a matrix summarizing which TFs are binding to which promoters
    TFs_matrix <- matrix(
        FALSE, 
        ncol = length(TFs), 
        nrow = length(proms)
    )
    colnames(TFs_matrix) <- TFs ; row.names(TFs_matrix) <- names(proms)
    for (TF in TFs) {
        hits <- subjectHits(findOverlaps(TFs_peaks[[TF]], proms))
        TFs_matrix[hits, TF] <- TRUE
    }
    # Integrate this to the promoters GRanges object
    proms$TFs <- apply(TFs_matrix, 1, function(ROW) {paste(colnames(TFs_matrix)[ROW], collapse = ',')})
    proms$TF_nb <- rowSums(TFs_matrix)
    proms$is_HOT <- ifelse(proms$TF_nb >= 18, TRUE, FALSE)
    cold_proms <- proms[!proms$is_HOT]

    clusters <- c("I", "I+H [1]", "I+H [2]", "Mix1", "Mix2", "Mix3", "Mix4", "Mix5")
    fisher_matrix <- data.frame(matrix(0, nrow = length(clusters), ncol = length(TFs))) ; colnames(fisher_matrix) = TFs ; row.names(fisher_matrix) = clusters
    pval_matrix <- data.frame(matrix(0, nrow = length(clusters), ncol = length(TFs))) ; colnames(pval_matrix) = TFs ; row.names(pval_matrix) = clusters
    pct_matrix <- data.frame(matrix(0, nrow = length(clusters), ncol = length(TFs))) ; colnames(pct_matrix) = TFs ; row.names(pct_matrix) = clusters
    for (CLUSTER in clusters) {
        # Get the cold promoters that are in the cluster CLUSTER
        cold_proms_in_cluster <- cold_proms[cold_proms$ageing_prom_cluster_label == CLUSTER]
        for (TF in TFs) {
            # Get the cold promoters bound by the transcription factor TF
            cold_proms_withTF <- sum(grepl(TF, cold_proms$TFs))
            # Get the cold promoters that are in the cluster CLUSTER and bound by the transcription factor TF
            cold_proms_in_cluster_withTF <- sum(grepl(TF, cold_proms_in_cluster$TFs))
            # Build a contingency matrix
            contingency_matrix <- matrix(
                c(
                    cold_proms_in_cluster_withTF, 
                    length(cold_proms_in_cluster) - cold_proms_in_cluster_withTF, 
                    cold_proms_withTF - cold_proms_in_cluster_withTF, 
                    length(cold_proms) - cold_proms_withTF - length(cold_proms_in_cluster) + cold_proms_in_cluster_withTF
                ), 
                nrow = 2
            )
            # Run a Fisher test and populate the result data.frames
            fisher <- fisher.test(contingency_matrix)
            fisher_matrix[CLUSTER, TF] <- fisher$estimate
            pval_matrix[CLUSTER, TF] <- fisher$p.val
            pct_matrix[CLUSTER, TF] <- cold_proms_in_cluster_withTF/length(cold_proms_in_cluster)*100
        }
    }
    # Filter the results to only see the significant ones
    fisher_matrix_curated <- t(fisher_matrix)
    fisher_matrix_curated[t(pval_matrix) > 0.01] <- 0 # Discard non-significant enrichments
    fisher_matrix_curated <- fisher_matrix_curated[rowSums(fisher_matrix_curated > 2) > 1,] # Discard TFs that don't have more than 2-fold enrichment in at least 1 cluster
    df <- fisher_matrix_curated %>% 
        data.frame() %>%
        setNames(clusters) %>% 
        rownames_to_column('TF') %>%
        gather('cluster', 'Odd_ratio', -TF)
    p <- ggplot(df, aes(y = TF, x = cluster)) +
        geom_tile(aes(fill = Odd_ratio)) + 
        scale_fill_gradientn(colours = c('white', 'orange', 'darkred')) + 
        theme_bw() + 
        labs(y = 'TFs', x = 'Promoter clusters in aging') + 
        theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
    ggsave('../results/TFs-enrichments_promoter-clusters.pdf', width = 5, height = 7)

Defining key TFs involved in aging?

[EXTRA] Finding motifs associated to each TFs

Most transcription factors recognize specific sequence motifs to bind to DNA. Using MEME web interface and the bed files located in data/, you can investigate which motifs are enriched in transcription factor binding sites (provided in the .bed files).

Save your work!!

Don’t forget to save your progress!

    save.image('../results/Day4.RData')