In ArnaudDroitLab/consensusSeekeR: Detection of consensus regions inside a group of experiences using genomic positions and genomic ranges

BiocStyle::markdown()
library(knitr)

Licensing and citing

This package and the underlying r Biocpkg("consensusSeekeR") code are distributed under the Artistic license 2.0. You are free to use and redistribute this software.

If you use this package for a publication, we would ask you to cite the following:

Samb R, Khadraoui K, Belleau P, et al. (2015) Using informative Multinomial-Dirichlet prior in a t-mixture with reversible jump estimation of nucleosome positions for genome-wide profiling. Statistical Applications in Genetics and Molecular Biology. Published online before print December 10, 2015. doi:10.1515/sagmb-2014-0098

Introduction

Genome data, such as genes, nucleosomes or single-nucleotide polymorphisms (SNPs) are linked to the genome by occupying either a range of positions or a single position on the sequence. Genomic related data integration is made possible by treating the data as ranges on the genome [@Lawrence2013]. Bioconductor has developed an infrastructure, including packages such as r Biocpkg("GenomicRanges"), r Biocpkg("IRanges") and r Biocpkg("GenomicFeatures"), which facilitate the integrative statistical analysis of range-based genomic data.

Ranges format is a convenient way for the analysis of diffent experimental genomic data. As an example, the peak calling step, in the analysis of ChIP-seq data, commonly generates NarrowPeak outputs. The NarrowPeak format, wich is used by the ENCODE project [@Dunham2012a], includes a peak position located inside a genomic range.

In genomic analysis where feature identification generate a position value surrounded by a genomic range, such as ChIP-Seq peaks and nucleosome postions, the replication of an experiment may result in slight differences between predicted values. Conciliation of the results can be difficult, especially when many replicates are done. One current approach used to identify consensus regions in a group of results consist of extracting the overlapping regions of the genomic ranges. This approach, when used on a large number of experiments, can miss, as a side effect, regions when one of the experiment has missing or slightly shift features. On the other hand, the use of the union of the regions can result in wide consensus ranges.

As an example, the Figure \@ref(fig:peakRegion) shows, using Integrative genomics viewer [@Interests2011], two ChIP-Seq peaks from ENCODE for the FOSL2 transcription factor (DCC accession: ENCFF002CFN). The data have been analyzed using MACS2 [@Zhang2008] with the default parameters and the q-value set to 0.05. The ChIP-Seq peak is a genomic feature that can be defined by a position value (the peak position) and a genome range (the enriched region). This example shows that the peak position is not necessarily at the center of the enriched region.

knitr::include_graphics('figures/peak_and_region.jpg')

The r Biocpkg("consensusSeekeR") package implements a novative way to identify consensus which use the features positions, instead of the most commonly used genomic ranges.

The consensusSeekeR package

The r Biocpkg("consensusSeekeR") package implements a novative way to identify consensus ranges in a group of experiments which generated position values surrounded by genomic ranges. The r Biocpkg("consensusSeekeR") package is characterized by its use of the position values, instead of the genomic ranges, to identify the consensus regions. The positions values have the double advantages of being, most of the time, the most important information from features and allowing creation of consensius regions of smaller ranges.

Using iterative steps on ordered features position values from all experiments, a window of fixed size (specified by user) with the current feature position as starting point is set. All features which reside inside the window are gathered to calculate a median feature position which is then used to recreate a new window. This time, the new window has twice the size fixed by user and its center is the median feature position. An update of the features located inside the window is done and the median feature position is recalculated. This step is repeated up to the moment that the set of features remains identical between two iterations. The final set of features positions is used to fix the central position of the consensus region. This final region must respect the minimum number of experiments with at least one feature inside it to be retained as a final consensus region. The minimum number of experiments is set by the user. At last, the consensus region can be extended or/and shrinked to fit the regions associated to the position values present inside. If new features positions are added during the consensus region resizing, the iterative steps are not reprocessed. It is possible that the extension step adds new features in the extended consensus region. However, those new features ranges won't be taken into account during the extension step.

Loading consensusSeekeR package

As with any R package, the r Biocpkg("consensusSeekeR") package should first be loaded with the following command:

library(consensusSeekeR)

Inputs

Positions and Ranges

The main function of the r Biocpkg("consensusSeekeR") is findConsensusPeakRegions. The mains inputs of the findConsensusPeakRegions function are:

1. a GRanges of the feature positions of all experiments with a metadata field called name.
2. a GRanges of the feature ranges for all experiments with a metadata field called name.

Beware that the GRanges of the feature ranges is only mandatory if the expandToFitPeakRegion parameter and/or the shrinkToFitPeakRegion parameter are set to TRUE.

Both inputs must satify these conditions:

• All rows of each GRanges must be named after the experiment source. All entries from the same experiment must be assigned the same name.
• Each feature must have one entry in both GRanges. The metadata field name is used to associate the feature position to its range.

This is an example showing how a metadata field name can easily be created and row names can be assigned:

### Load dataset
data(A549_FOSL2_01_NarrowPeaks_partial)

### Remove dataset metadata field "name"
A549_FOSL2_01_NarrowPeaks_partial$name        <- NULL
A549_FOSL2_01_NarrowPeaks_partial$score       <- NULL
A549_FOSL2_01_NarrowPeaks_partial$qValue      <- NULL
A549_FOSL2_01_NarrowPeaks_partial$pValue      <- NULL
A549_FOSL2_01_NarrowPeaks_partial$signalValue <- NULL
A549_FOSL2_01_NarrowPeaks_partial$peak        <- NULL

### Initial dataset without metadata field
head(A549_FOSL2_01_NarrowPeaks_partial, n = 3)

### Adding a new metadata field "name" unique to each entry
A549_FOSL2_01_NarrowPeaks_partial$name <- paste0("FOSL2_01_entry_", 
                                1:length(A549_FOSL2_01_NarrowPeaks_partial))

### Assign the same row name to each entry
names(A549_FOSL2_01_NarrowPeaks_partial) <- rep("FOSL2_01", 
                                length(A549_FOSL2_01_NarrowPeaks_partial))

### Final dataset with metadata field 'name' and row names
head(A549_FOSL2_01_NarrowPeaks_partial, n = 3)

### Remove dataset
rm(A549_FOSL2_01_NarrowPeaks_partial)

Chromosomes information

The chromosomes information is mandatory. It ensures that the consensus regions do not exceed the length of the chromosomes.

The chromosomes information is contained inside a Seqinfo object. The information from some UCSC genomes can be fetched automatically using the r Biocpkg("GenomeInfoDb") package.

### Import library
library(GenomeInfoDb)

### Get the information for Human genome version 19
hg19Info <- Seqinfo(genome="hg19")

### Subset the object to keep only the analyzed chromosomes
hg19Subset <- hg19Info[c("chr1", "chr10", "chrX")]

A Seqinfo object can also be created using the chromosomes information specific to the analyzed genome.

### Create an Seqinfo Object
chrInfo <- Seqinfo(seqnames=c("chr1", "chr2", "chr3"),
            seqlengths=c(1000, 2000, 1500), isCircular=c(FALSE, FALSE, FALSE),
            genome="BioconductorAlien")

if (exists("chrInfo", inherits = FALSE)) rm(chrInfo)
if (exists("hg19Subset", inherits = FALSE)) rm(hg19Subset)
if (exists("hg19Info", inherits = FALSE)) rm(hg19Info)

Read NarrowPeak files

The NarrowPeak format is often used to provide called peaks of signal enrichment based on pooled, normalized data. The r Biocpkg("rtracklayer") package has functions which faciliates the loading of NarrowPeak files.

Since the main function of the r Biocpkg("consensusSeekeR") package needs 2 GRanges objects, some manipulations are needed to create one GRanges for the regions and one GRanges for the peaks.

### Load the needed packages
library(rtracklayer)
library(GenomicRanges)

### Demo file contained within the consensusSeekeR package
file_narrowPeak <- system.file("extdata",
        "A549_FOSL2_ENCSR000BQO_MZW_part_chr_1_and_12.narrowPeak", package = "consensusSeekeR")

### Information about the extra columns present in the file need
### to be specified
extraCols <- c(signalValue = "numeric", pValue = "numeric", qValue = "numeric", peak = "integer")

### Create genomic ranges for the regions
regions <- import(file_narrowPeak, format = "BED", extraCols = extraCols)

### Create genomic ranges for the peaks
peaks           <-  regions
ranges(peaks)   <-  IRanges(start = (start(regions) + regions$peak), width = rep(1, length(regions$peak)))

### First rows of each GRanges object
head(regions, n = 2)
head(peaks, n = 2)

rm(peaks)
rm(regions)

A simpler way is to use the readNarrowPeakFile function of the r Biocpkg("consensusSeekeR") package which generates both the peaks and the narrowPeak GRanges.

library(consensusSeekeR)

### Demo file contained within the consensusSeekeR package
file_narrowPeak <- system.file("extdata",
    "A549_FOSL2_ENCSR000BQO_MZW_part_chr_1_and_12.narrowPeak", package = "consensusSeekeR")

### Create genomic ranges for both the regions and the peaks
result <- readNarrowPeakFile(file_narrowPeak, extractRegions = TRUE, extractPeaks = TRUE)

### First rows of each GRanges object
head(result$narrowPeak, n = 2)
head(result$peak, n = 2)

rm(result)

Case study: nucleosome positioning

Global gene expression patterns are established and maintained by the concerted actions of Transcription Factors (TFs) and the proteins that constitute chromatin. The key structural element of chromatin is the nucleosome, which consists of an octameric histone core wrapped by 147 bps of DNA and connected to its neighbor by approximately 10-80 pbs of linker DNA [@Kornberg1999]. Nucleosome occupancy and positioning have been proved to be dynamic. It also has a major impact on expression, regulation, and evolution of eukaryotic genes [@Jiang2015].

Comparing nucleosome positioning results from different software

With the development of Next-generation sequencing, nucleosome positioning using MNase-Seq data or MNase- or sonicated- ChIP-Seq data combined with either single-end or paired-end sequencing have evolved as popular techniques. Software such as PING [@Woo2013] and NOrMAL [@Polishko2012], generates output which contains the positions of the predicted nucleosomes, which simply are one base pair positions on the reference genome. This position represents the center of the predicted nucleosome. A range of $\pm$ 73 bps is usually superposed to the predicted nucleosome to repesent the nucleosome occupancy.

First, the r Biocpkg("consensusSeekeR") package must be loaded.

library(consensusSeekeR)

The datasets, which are included in the r Biocpkg("consensusSeekeR") package, have to be loaded. Those include results obtained using syntethic reads distributed following a normal distribution with a variance of 20 from three different nucleosome positioning software: PING [@Woo2013], NOrMAL [@Polishko2012] and NucPosSimulator [@Schopflin2013]. The genomic ranges have been obtained by adding $\pm$ 73 bps to the detected positions.

### Loading dataset from NOrMAL
data(NOrMAL_nucleosome_positions) ; data(NOrMAL_nucleosome_ranges)

### Loading dataset from PING
data(PING_nucleosome_positions) ; data(PING_nucleosome_ranges)

### Loading dataset from NucPosSimulator
data(NucPosSimulator_nucleosome_positions) ; data(NucPosSimulator_nucleosome_ranges)

rownames(NOrMAL_nucleosome_positions) <- NULL
rownames(NOrMAL_nucleosome_ranges)    <- NULL

For the positions and ranges dataset from the same software, the name field is paired to ensure that each position can be associated to its range. The metadata field name must be unique to each feature for all datasets.

### Each entry in the positions dataset has an equivalent metadata 
### "name" entry in the ranges dataset
head(NOrMAL_nucleosome_positions, n = 2)

head(NOrMAL_nucleosome_ranges, n = 2)

To be able to identify all entries from the same software, each row of the dataset has to be assigned a name. All positions and ranges from the same source must be assigned identical row names. In this exemple, datasets are going to be identified by the name of their source software.

### Assigning software name "NOrMAL" 
names(NOrMAL_nucleosome_positions) <- rep("NOrMAL", length(NOrMAL_nucleosome_positions))
names(NOrMAL_nucleosome_ranges) <- rep("NOrMAL", length(NOrMAL_nucleosome_ranges))

### Assigning experiment name "PING"
names(PING_nucleosome_positions) <- rep("PING", length(PING_nucleosome_positions))
names(PING_nucleosome_ranges) <- rep("PING", length(PING_nucleosome_ranges))

### Assigning experiment name "NucPosSimulator"
names(NucPosSimulator_nucleosome_positions) <- rep("NucPosSimulator", 
                                length(NucPosSimulator_nucleosome_positions))
names(NucPosSimulator_nucleosome_ranges) <- rep("NucPosSimulator", 
                                length(NucPosSimulator_nucleosome_ranges))

### Row names are unique to each software
head(NOrMAL_nucleosome_positions, n = 2)

head(PING_nucleosome_positions, n = 2)

head(NucPosSimulator_nucleosome_positions, n = 2)

The consensus regions for chromosome 1 only are calculated with a defaut region size of 50 bases pairs (2 * extendingSize) The regions are extended to include all nucleosome regions (expandToFitPeakRegion = TRUE and shrinkToFitPeakRegion = TRUE). To be retained as a consensus region, nucleosomes from at least 2 software must be present in the region (minNbrExp = 2). .

### Only choromsome 1 is going to be analyzed
chrList <- Seqinfo("chr1", 135534747, NA)

### Find consensus regions with both replicates inside it
results <- findConsensusPeakRegions(
            narrowPeaks = c(NOrMAL_nucleosome_ranges, PING_nucleosome_ranges,
                                NucPosSimulator_nucleosome_ranges),
            peaks = c(NOrMAL_nucleosome_positions, PING_nucleosome_positions,
                                NucPosSimulator_nucleosome_positions),
            chrInfo = chrList,
            extendingSize = 25,
            expandToFitPeakRegion = TRUE,
            shrinkToFitPeakRegion = TRUE,
            minNbrExp = 2,
            nbrThreads = 1)

The output of findConsensusPeakRegions function is a list containing an object call and an object consensusRanges. The object call contains the matched call while the object consensusRanges is a GRanges containing the consensus regions.

### Print the call 
results$call

### Print the 3 first consensus regions
head(results$consensusRanges, n = 3)

A total of r length(results$consensusRanges) consensus regions have been found. An exemple of the consensus regions (in dark blue) is shown in Figure \@ref(fig:nucleosomes) using Integrative genomics viewer [@Interests2011]:

knitr::include_graphics('figures/nucleosomes.jpg')

### Remove dataset from NOrMAL
rm(NOrMAL_nucleosome_positions)
rm(NOrMAL_nucleosome_ranges)

### Remove dataset from PING
rm(PING_nucleosome_positions)
rm(PING_nucleosome_ranges)

### Remove dataset from NucPosSimulator
rm(NucPosSimulator_nucleosome_positions)
rm(NucPosSimulator_nucleosome_ranges)

Case study: ChIP-Seq data

Next-generation DNA sequencing coupled with chromatin immunoprecipitation (ChIP-seq) has changed the ability to interrogate the genomic landscape of histone modifications, transcriptional cofactors and transcription-factors binding in living cells [@Mundade2014]. Consortium, such as ENCODE have developed and are constantly updating a set of standards and guidelines for ChIP-Seq experiments [@Landt2012].

ChIP-seq combines chromatin immunoprecipitation (ChIP) with massively parallel DNA sequencing. The obtained sequence reads are first mapped to the reference genome of the organism used in the experiments. Binding sites are then detected using software specialized in transcript factor binding sites identification, such as MACS2 [@Zhang2008] and PeakRanger [@Feng2011]. Peaks are defined as a single base pair position while statistically enriched regions are defined as genomic ranges.

ChIP-Seq replicates from one experiment

The Encyclopedia of DNA Elements (ENCODE) Consortium is an international collaboration of research groups funded by the National Human Genome Research Institute. The ENCODE website is a portal giving access to the data generated by the ENCODE Consortium. The amount of data gathered is extensive. Moreover, for some experiments, more than one ChIP-Seq replicate is often available.

The software used to identify transcript factor binding sites generally generates a peak position and an enriched region for each binding site. However, it is quite unlikely that the exact peak position is exactly the same across replicates. Even more, there is not yet a consensus on how to analyze multiple-replicate ChIP-seq samples [@Yang2014].

The r Biocpkg("consensusSeekeR") package can be used to identify consensus regions for two or more replicates ChIP-Seq samples. The consensus regions are being found by using the peak positions.

The transcription factor binding for the CTCF transcription factor have been analyzed and 2 replicates are available in BAM files format on ENCODE website (DCCs: ENCFF000MYJ and ENCFF000MYN). The NarrowPeaks were generated using MACS2 [@Zhang2008] with the default parameters and the q-value set to 0.05.

To simplify this demo, only part of genome hg19, chr1:246000000-249250621 and chr10:10000000-12500000, have been retained in the datasets.

First, the r Biocpkg("consensusSeekeR") package must be loaded.

library(consensusSeekeR)

The datasets, which are included in the r Biocpkg("consensusSeekeR") package, have to be loaded.

### Loading datasets
data(A549_CTCF_MYN_NarrowPeaks_partial) ; data(A549_CTCF_MYN_Peaks_partial)
data(A549_CTCF_MYJ_NarrowPeaks_partial) ; data(A549_CTCF_MYJ_Peaks_partial)

To be able to identify data from the same source, each row of the dataset has to be assigned a source name. Beware that NarrowPeak and Peak datasets from the same source must be assigned identical names. In this exemple, datasets are replicates of the same experiment. So, the names "rep01" and "rep02" are going to be assigned to each dataset.

### Assigning experiment name "rep01" to the first replicate
names(A549_CTCF_MYJ_NarrowPeaks_partial) <- rep("rep01", length(A549_CTCF_MYJ_NarrowPeaks_partial))
names(A549_CTCF_MYJ_Peaks_partial)       <- rep("rep01", length(A549_CTCF_MYJ_Peaks_partial))

### Assigning experiment name "rep02" to the second replicate
names(A549_CTCF_MYN_NarrowPeaks_partial) <- rep("rep02", length(A549_CTCF_MYN_NarrowPeaks_partial))
names(A549_CTCF_MYN_Peaks_partial)       <- rep("rep02", length(A549_CTCF_MYN_Peaks_partial))

The consensus regions for chromosome 10 only are calculated with a defaut region size of 200 bases pairs (2 * extendingSize) The regions are extended to include all peaks regions (expandToFitPeakRegion = TRUE and shrinkToFitPeakRegion = TRUE). A peak from both replicates must be present withinin a region for it to be retained as a consensus region.

### Only choromsome 10 is going to be analyzed
chrList <- Seqinfo("chr10", 135534747, NA)

### Find consensus regions with both replicates inside it
results <- findConsensusPeakRegions(
            narrowPeaks = c(A549_CTCF_MYJ_NarrowPeaks_partial, A549_CTCF_MYN_NarrowPeaks_partial),
            peaks = c(A549_CTCF_MYJ_Peaks_partial, A549_CTCF_MYN_Peaks_partial),
            chrInfo = chrList,
            extendingSize = 100,
            expandToFitPeakRegion = TRUE,
            shrinkToFitPeakRegion = TRUE,
            minNbrExp = 2,
            nbrThreads = 1)

The output of findConsensusPeakRegions function is a list containing an object call and an object conesensusRanges. The object call contains the matched call while the object conesensusRanges is a GRanges containing the consensus regions.

### Print the call 
results$call

### Print the 3 first consensus regions
head(results$consensusRanges, n = 3)

A total of r length(results$consensusRanges) consensus regions have been found. An example of the consensus regions (in green) is shown in Figure \@ref(fig:ctcf), using Integrative genomics viewer [@Interests2011].

knitr::include_graphics('figures/ctcf_consensus.jpg')

### Remove datasets
rm(A549_CTCF_MYN_NarrowPeaks_partial)
rm(A549_CTCF_MYN_Peaks_partial)
rm(A549_CTCF_MYJ_NarrowPeaks_partial)
rm(A549_CTCF_MYJ_Peaks_partial)

ChIP-Seq data from multiple experiments

The r Biocpkg("consensusSeekeR") package can also be used to identify consensus regions for two or more ChIP-Seq samples from multiple experiments. The peak positions are the feature used to identify the consensus regions.

The transcription factor binding for the NR3C1 transcription factor have been analyzed in more than one experiment. For each experiment, replicates have been analyzed together using the irreproducible discovery rate method [@Li2011]. Results are available in bed narrowPeak format on ENCODE website (DCCs: ENCFF002CFQ, ENCFF002CFR and ENCFF002CFS) [@Dunham2012a].

To simplify this demo, only part of genome hg19, chr2:40000000-50000000 and chr3:10000000-13000000, have been retained in the datasets.

First, the r Biocpkg("consensusSeekeR") package must be loaded.

library(consensusSeekeR)

The datasets, which are included in the r Biocpkg("consensusSeekeR") package, have to be loaded.

### Loading datasets
data(A549_NR3C1_CFQ_NarrowPeaks_partial) ; data(A549_NR3C1_CFQ_Peaks_partial)
data(A549_NR3C1_CFR_NarrowPeaks_partial) ; data(A549_NR3C1_CFR_Peaks_partial)
data(A549_NR3C1_CFS_NarrowPeaks_partial) ; data(A549_NR3C1_CFS_Peaks_partial)

To be able to identify data from the same source, each row of the dataset has to be assigned an experiment name. Beware that NarrowPeak and Peak datasets from the same source must be assigned identical names. In this exemple, datasets are coming from different experiments for the same transcription factor. So, the short name of each experiment "ENCFF002CFQ", "ENCFF002CFR" and "ENCFF002CFS" is going to be assigned to each dataset.

### Assign experiment name "ENCFF002CFQ" to the first experiment
names(A549_NR3C1_CFQ_NarrowPeaks_partial) <- rep("ENCFF002CFQ", 
                                length(A549_NR3C1_CFQ_NarrowPeaks_partial))
names(A549_NR3C1_CFQ_Peaks_partial) <- rep("ENCFF002CFQ", 
                                length(A549_NR3C1_CFQ_Peaks_partial))

### Assign experiment name "ENCFF002CFQ" to the second experiment
names(A549_NR3C1_CFR_NarrowPeaks_partial) <- rep("ENCFF002CFR", 
                                length(A549_NR3C1_CFR_NarrowPeaks_partial))
names(A549_NR3C1_CFR_Peaks_partial)       <- rep("ENCFF002CFR", 
                                length(A549_NR3C1_CFR_Peaks_partial))

### Assign experiment name "ENCFF002CFQ" to the third experiment
names(A549_NR3C1_CFS_NarrowPeaks_partial) <- rep("ENCFF002CFS", 
                                length(A549_NR3C1_CFS_NarrowPeaks_partial))
names(A549_NR3C1_CFS_Peaks_partial)       <- rep("ENCFF002CFS", 
                                length(A549_NR3C1_CFS_Peaks_partial))

In ENCODE bed narrowPeak format, entries don't have a specific metadata field called name. So, to be able to use the findConsensusPeakRegions function, specific names must manually be added to each entry.

### Assign specific name to each entry of to first experiment
### NarrowPeak name must fit Peaks name for same experiment
A549_NR3C1_CFQ_NarrowPeaks_partial$name <- paste0("NR3C1_CFQ_region_", 
                                1:length(A549_NR3C1_CFQ_NarrowPeaks_partial))
A549_NR3C1_CFQ_Peaks_partial$name       <- paste0("NR3C1_CFQ_region_",
                                1:length(A549_NR3C1_CFQ_NarrowPeaks_partial))

### Assign specific name to each entry of to second experiment
### NarrowPeak name must fit Peaks name for same experiment
A549_NR3C1_CFR_NarrowPeaks_partial$name <- paste0("NR3C1_CFR_region_", 
                                1:length(A549_NR3C1_CFR_NarrowPeaks_partial))
A549_NR3C1_CFR_Peaks_partial$name       <- paste0("NR3C1_CFR_region_", 
                                1:length(A549_NR3C1_CFR_Peaks_partial))

### Assign specific name to each entry of to third experiment
### NarrowPeak name must fit Peaks name for same experiment
A549_NR3C1_CFS_NarrowPeaks_partial$name <- paste0("NR3C1_CFS_region_", 
                                1:length(A549_NR3C1_CFS_NarrowPeaks_partial))
A549_NR3C1_CFS_Peaks_partial$name       <- paste0("NR3C1_CFS_region_", 
                                1:length(A549_NR3C1_CFS_Peaks_partial))

The consensus regions for chromosome 2 only are calculated with a defaut region size of 400 bases pairs (2 * extendingSize) The regions are not extended to include all peaks regions but are shrinked when exceeding peaks regions (expandToFitPeakRegion = FALSE and shrinkToFitPeakRegion = TRUE). A peak from 2 out of 3 experiments must be present in a region for it to be retained as a consensus region.

### Only choromsome 2 is going to be analyzed
chrList <- Seqinfo("chr2", 243199373, NA)

### Find consensus regions with both replicates inside it
results <- findConsensusPeakRegions(
            narrowPeaks = c(A549_NR3C1_CFQ_NarrowPeaks_partial, 
                        A549_NR3C1_CFR_NarrowPeaks_partial,
                        A549_NR3C1_CFS_NarrowPeaks_partial),
            peaks = c(A549_NR3C1_CFQ_Peaks_partial, 
                        A549_NR3C1_CFR_Peaks_partial,
                        A549_NR3C1_CFS_Peaks_partial),
            chrInfo = chrList,
            extendingSize = 200,
            expandToFitPeakRegion = FALSE,
            shrinkToFitPeakRegion = TRUE,
            minNbrExp = 2,
            nbrThreads = 1)

The output of findConsensusPeakRegions function is a list containing an object call and an object consensusRanges. The object call contains the matched call while the object consensusRanges is a GRanges containing the consensus regions.

### Print the call 
results$call

### Print the first 3 consensus regions
head(results$consensusRanges, n = 3)

A total of r length(results$consensusRanges) consensus regions have been found. A example of the consensus regions (in green) is shown in Figure \@ref(fig:NR3C1) using Integrative genomics viewer [@Interests2011].

knitr::include_graphics('figures/NR3C1_consensus.jpg')

### Remove datasets
rm(A549_NR3C1_CFQ_NarrowPeaks_partial)
rm(A549_NR3C1_CFQ_Peaks_partial)
rm(A549_NR3C1_CFR_NarrowPeaks_partial)
rm(A549_NR3C1_CFR_Peaks_partial)
rm(A549_NR3C1_CFS_NarrowPeaks_partial)
rm(A549_NR3C1_CFS_Peaks_partial)

Parameters

Effect of the shrinkToFitPeakRegion parameter

The shrinkToFitPeakRegion allows the resizing of the consensus region to fit the minimum regions of the included features when those values are included inside the initial consensus region. When the extendingSize parameter is large, the effect can be quite visible on the final consensus regions. For exemple, Figure \@ref(fig:shrink) shows the same region, from the NR3C1 example with extendingSize of 200, when the shrinkToFitPeakRegion is set to TRUE (green color) and to FALSE (orange color) using Integrative genomics viewer [@Interests2011].

knitr::include_graphics('figures/shrink.jpg')

Effect of the expandToFitPeakRegion parameter

The expandToFitPeakRegion allows the resizing of the consensus region to fit the maximum of the included features when those values are outside the initial consensus region. When the extendingSize parameter is small, the effect can be quite visible on the final consensus regions. For exemple, the following figure shows the same region, from the CTCF example with extendingSize of 100, when the shrinkToFitPeakRegion is set to FALSE (orange color) and to TRUE (green color).

knitr::include_graphics('figures/extendNew.jpg')

Effect of the extendingSize parameter

The value of the extendingSize parameter can affect the final number of consensus regions. While small extendingSize value can miss some regions, large extendingSize value can gather consensus regions. Testing a range of extendingSize parameters can be an option worth considering.

As an example, the number of consensus regions obtained with different values of extendingSize is calculated.

``` {r sizeEffect, collapse=TRUE, eval=FALSE}

Set different values for the extendingSize parameter

size <- c(1, 10, 50, 100, 300, 500, 750, 1000)

Only chrompsome 10 is going to be analyzed

chrList <- Seqinfo("chr10", 135534747, NA)

Find consensus regions using all the size values

resultsBySize <- lapply(size, FUN = function(size) findConsensusPeakRegions( narrowPeaks = c(A549_CTCF_MYJ_NarrowPeaks_partial, A549_CTCF_MYN_NarrowPeaks_partial), peaks = c(A549_CTCF_MYJ_Peaks_partial, A549_CTCF_MYN_Peaks_partial), chrInfo = chrList, extendingSize = size, expandToFitPeakRegion = TRUE, shrinkToFitPeakRegion = TRUE, minNbrExp = 2, nbrThreads = 1))

Extract the number of consensus regions obtained for each extendingSize

nbrRegions <- mapply(resultsBySize, FUN = function(x) return(length(x$consensusRanges)))

A graph can be used to visualize the variation of the number of consensus 
regions in function of the `extendingSize` 
parameter (see Figure \@ref(fig:sizeEffectG)).

```r

library(ggplot2)

data <- data.frame(extendingSize = size, nbrRegions = nbrRegions)

ggplot(data, aes(extendingSize, nbrRegions)) + scale_x_log10("Extending size") + 
        stat_smooth(se = FALSE, method = "loess", size=1.4) +
        ylab("Number of consensus regions") +
        ggtitle(paste0("Number of consensus regions in function of the extendingSize"))

knitr::include_graphics('figures/sizeEffect.jpg')

Parallelizing consensusSeekeR

Due to the size of the analyzed genomes, the findConsensusPeakRegions function can take a while to process. However, a job can be separated by chromosome and run in parallel. This takes advantage of multiple processors and reduce the total execution time. The number of threads used can be set with the nbrThreads parameter in the findConsensusPeakRegions function.

### Load data
data(A549_FOSL2_01_NarrowPeaks_partial) ; data(A549_FOSL2_01_Peaks_partial)
data(A549_FOXA1_01_NarrowPeaks_partial) ; data(A549_FOXA1_01_Peaks_partial)

### Assigning name "FOSL2" 
names(A549_FOSL2_01_NarrowPeaks_partial) <- rep("FOSL2", 
                                    length(A549_FOSL2_01_NarrowPeaks_partial))
names(A549_FOSL2_01_Peaks_partial)       <- rep("FOSL2", 
                                    length(A549_FOSL2_01_Peaks_partial))

### Assigning name "FOXA1"
names(A549_FOXA1_01_NarrowPeaks_partial) <- rep("FOXA1", 
                                    length(A549_FOXA1_01_NarrowPeaks_partial))
names(A549_FOXA1_01_Peaks_partial)       <- rep("FOXA1", 
                                    length(A549_FOXA1_01_Peaks_partial))

### Two chromosomes to analyse
chrList <- Seqinfo(paste0("chr", c(1,10)), c(249250621, 135534747), NA)

### Find consensus regions using 2 threads
results <- findConsensusPeakRegions(
            narrowPeaks = c(A549_FOSL2_01_NarrowPeaks_partial, 
                            A549_FOXA1_01_Peaks_partial),
            peaks = c(A549_FOSL2_01_Peaks_partial, 
                            A549_FOXA1_01_NarrowPeaks_partial),
            chrInfo = chrList, extendingSize = 200, minNbrExp = 2, 
            expandToFitPeakRegion = FALSE, shrinkToFitPeakRegion = FALSE,
            nbrThreads = 4)

Acknowledgment

We thank Imène Boudaoud for her advice on the vignette content.

Session info

Here is the output of sessionInfo() on the system on which this document was compiled:

sessionInfo()

References

ArnaudDroitLab/consensusSeekeR documentation built on Feb. 6, 2020, 6:45 p.m.