Cogito 1.0.0
Biological studies often consist of multiple conditions which are examined with different laboratory set ups like RNA-sequencing or ChIP-sequencing. To get an overview about the whole resulting data set, Cogito provides an automated and complete report about all samples and basic comparisons between all different samples. This overview can be used as documentation about the data set or as starting point for further custom analysis.
Cogito is not meant to do detailed genetic or epigenetic analysis. But provides a reproducible and clear report about data sets consisting of samples of different conditions and mesuared with different laboratory base technologies.
The package can be installed via Bioconductor
and loaded into R by typing
BiocManager::install("Cogito")
library("Cogito")
into the R console.
The workflow of Cogito consists of two main functions:
aggregateRanges(ranges, configfile, organism, name, verbose)
summarizeRanges(aggregated.ranges, verbose)
The process is demonstrated using a small murine example. The data set from
King et al.1 4 was downloaded from the NCBI GEO database under
accession number
GSE77004 and
preprocessed2 for details on preprocessing of the data have a look at the
manual pages of the single data sets to GRanges
/GRangesList
objects and
restricted to chromosome 5 to save space and processing time.
This example data set is loaded with the package and consists of three data objects.
The first GRanges
object contains gene expression values from RNA-sequencing
with 1195 ranges and 11
columns with the corresponding expression values. Each column represents an
experimental condition.
head(MurEpi.RNA.small[, 1:3])
## GRanges object with 6 ranges and 3 metadata columns:
## seqnames ranges strand | J1 TKO TKO3a1_1
## <Rle> <IRanges> <Rle> | <numeric> <numeric> <numeric>
## [1] chr5 101234211-101249984 - | 7.904 7.142 8.420
## [2] chr5 73039035-73062317 + | 0.101 0.344 0.431
## [3] chr5 6769030-6876523 - | 0.000 0.000 0.011
## [4] chr5 115303673-115333668 + | 0.150 0.000 0.022
## [5] chr5 3657004-3680325 + | 0.000 0.000 0.000
## [6] chr5 3583978-3596585 - | 8.752 6.195 8.677
## -------
## seqinfo: 35 sequences (1 circular) from mm9 genome
The second GRanges
object contains information about the methylation status
from RRBS. The object has 147644 ranges and
14 columns with the corresponding methylation
status. Where there is not information about the methylation status the value
is NA
. Each column represents an experimental condition.
head(MurEpi.RRBS.small[, 1:3])
## GRanges object with 6 ranges and 3 metadata columns:
## seqnames ranges strand | J1_1.meth J1_2.meth DKO.meth
## <Rle> <IRanges> <Rle> | <ordered> <ordered> <ordered>
## [1] chr5 3021138-3021139 * | high medium low
## [2] chr5 3021182-3021183 * | high medium low
## [3] chr5 3021357-3021358 * | high NA low
## [4] chr5 3021446-3021447 * | high NA medium
## [5] chr5 3025756-3025757 * | high high medium
## [6] chr5 3025763-3025764 * | high high low
## -------
## seqinfo: 35 sequences (1 circular) from mm9 genome
The third object is a GRangesList
object containing
36 GRanges
object. Each of this objects
represents an experimental condition and contains between
10 and
2468 ranges with one column of
attached value each. This value is the peak score resulting from the peak
calling with homer findPeaks3 4.
head(MurEpi.ChIP.small[[1]])
## GRanges object with 6 ranges and 1 metadata column:
## seqnames ranges strand | findPeaksScore
## <Rle> <IRanges> <Rle> | <numeric>
## [1] chr5 125869822-125870015 * | 282
## [2] chr5 114221367-114221560 * | 273
## [3] chr5 115729732-115729925 * | 269
## [4] chr5 122734457-122734650 * | 265
## [5] chr5 124875261-124875454 * | 275
## [6] chr5 117565637-117565830 * | 263
## -------
## seqinfo: 35 sequences (1 circular) from mm9 genome
The data set consists of the following samples:
RNA-seq | RRBS | ChIP-seq | |
---|---|---|---|
J1 | 1 | 2 | 5 |
TKO | 1 | 2 | 5 |
TKO3a1 | 2 | 2 | 4 |
TKO3a2 | 2 | 2 | 5 |
TKO3b1 | 1 | 2 | 4 |
DKO | 1 | 1 | 4 |
DKO3a1 | 1 | 1 | 3 |
DKO3a2 | 1 | 1 | 3 |
DKO3b1 | 1 | 1 | 3 |
The function aggregateRanges
is called with the example data set and the
corresponding genomic information.
mm9 <- TxDb.Mmusculus.UCSC.mm9.knownGene::TxDb.Mmusculus.UCSC.mm9.knownGene
example.dataset <- list(ChIP = MurEpi.ChIP.small,
RNA = MurEpi.RNA.small,
RRBS = MurEpi.RRBS.small)
aggregated.ranges <- aggregateRanges(ranges = example.dataset,
organism = mm9,
name = "murine.example")
## 84 genes were dropped because they have exons located on both strands
## of the same reference sequence or on more than one reference sequence,
## so cannot be represented by a single genomic range.
## Use 'single.strand.genes.only=FALSE' to get all the genes in a
## GRangesList object, or use suppressMessages() to suppress this message.
The function aggregates all provided data to the genes of the genome, given in
the parameter organism
: all single IRanges
are assigned to the closest
gene within a predefined default maximal distance of 100.000bp. This
parameter can be changed in the configuration file if necessary. Consequently,
the columns of attached values of each provided IRanges
object are assigned
to the corresponding gene. The result is one single GRanges
object, each row
representing a gene and each column a sample of the provided input data.
The function aggregateRanges
returns a list
object containing this single
GRanges
objects, i.e. a GRanges
object with rows representing genes with
columns of attached values mcols
where each column contains values from a
specific experimental condition (such as wildtype) and a specific underlying
base technology (such as RNA-seq expression).
head(aggregated.ranges$genes[, c(1, 2:3, 13:14, 27:28)])
## GRanges object with 6 ranges and 7 metadata columns:
## seqnames ranges strand | gene_id RNA.J1 RNA.TKO
## <Rle> <IRanges> <Rle> | <character> <numeric> <numeric>
## 12571 chr5 3343893-3522225 + | 12571 5.647 3.996
## 68152 chr5 3543833-3570546 + | 68152 NA NA
## 77036 chr5 3571716-3584341 + | 77036 3.176 1.396
## 71382 chr5 3596066-3637101 + | 71382 5.631 5.225
## 75010 chr5 3657004-3680325 + | 75010 0.000 0.000
## 79264 chr5 3803165-3844515 + | 79264 NA NA
## RRBS.J1_1.meth RRBS.J1_2.meth ChIP.H3K4me3_J1.findPeaksScore
## <ordered> <ordered> <numeric>
## 12571 low low 111
## 68152 low low 199
## 77036 high high NA
## 71382 low low 212
## 75010 high medium NA
## 79264 low low 199
## ChIP.H3K4me3_TKO.findPeaksScore
## <numeric>
## 12571 83
## 68152 164
## 77036 NA
## 71382 186
## 75010 NA
## 79264 163
## -------
## seqinfo: 35 sequences (1 circular) from mm9 genome
The return value of the function aggregateRanges
as well contains information
about the supposed underlying base technologies and conditions.
lapply(aggregated.ranges$config$technologies, head, 3)
## $ChIP
## [1] "ChIP.H3K4me3_J1.findPeaksScore" "ChIP.H3K4me3_TKO.findPeaksScore"
## [3] "ChIP.H3K4me3_TKO3a2.findPeaksScore"
##
## $RNA
## [1] "RNA.J1" "RNA.TKO" "RNA.TKO3a1_1"
##
## $RRBS
## [1] "RRBS.J1_1.meth" "RRBS.J1_2.meth" "RRBS.DKO.meth"
head(lapply(aggregated.ranges$config$conditions, head, 3), 3)
## $J1
## [1] "RNA.J1" "RRBS.J1_1.meth" "RRBS.J1_2.meth"
##
## $TKO
## [1] "RNA.TKO" "RRBS.TKO_1.meth" "RRBS.TKO_2.meth"
##
## $TKO3a1
## [1] "RNA.TKO3a1_1" "RNA.TKO3a1_2" "RRBS.TKO3a1_1.meth"
To advance in the workflow, these two list
s should be carefully checked and
corrected if necessary. Also it is possible to add user defined groups of
conditions to include them in the following analysis.
summarizeRanges(aggregated.ranges = aggregated.ranges)
The function summarizeRanges
does not has a return value but produces three
output files:
The main result is the report about the given data in a pdf. The report is divided into four chapters. The introduction holds general information about the underlying data set as the numbers of conditions (wildtype etc.) and used base technologies (ChIP-seq, RNA-seq, etc.).
The first chapter contains summaries about each single column of attached values, consisting of a location and a dispersion parameter as well as a plot for a visual impression. There are 61 samples in total included in the presented murine example data set.
The second chapter summarizes groups of columns of attached values which have the same condition and the same base technology, for example if there are two RNA-sequencing replicates of the condition wildtpye present. In the given murine example data set there are among others 2 samples of methylation status examined with RRBS in the condition J1, 5 samples of condition J1 of ChIP-seq experiments and 2 RNA-seq samples of condition TKO3a1, see table 1. Consequentely, this results in 16 plots in total.
The third chapter summarizes groups of columns of attached values which are examined by the same base technology but do not necessarily have the same underlying condition. This is visualized in an appropriate plot. The presented murine example has three involved base technologies: RNA-seq, ChIP-seq and RRBS.
Finally the fourth chapter compares every column of attached value to every
other column of attached value. The comparison is visualized with an appropriate
plot and a correlation test is made. Because this section may be long, the
report concentrates on comparisons where a significant correlation is found
(corrected p-value < 0.01).
In the exemplary murine data set there are
61 columns of attached values, so this
results in (61*60)/2=1830
comparisons, but of not all of them show a
significant correlation.
Besides the pdf report the function summarizeRanges
also produces a rmd file
with which the user may customize the report or take it as a starting point for
further analysis, as well as a RData file which holds all used processed data.
King AD, Huang K, Rubbi L, Liu S, Wang CY, Wang Y, Pellegrini M, Fan G. Reversible Regulation of Promoter and Enhancer Histone Landscape by DNA Methylation in Mouse Embryonic Stem Cells. Cell Rep. 2016 Sep 27;17(1):289-302. doi: 10.1016/j.celrep.2016.08.083
Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, Glass CK. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010 May 28;38(4):576-89. doi: 10.1016/j.molcel.2010.05.004
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] Cogito_1.0.0 entropy_1.3.1 GenomicFeatures_1.46.0
## [4] AnnotationDbi_1.56.0 Biobase_2.54.0 jsonlite_1.7.2
## [7] GenomicRanges_1.46.0 GenomeInfoDb_1.30.0 IRanges_2.28.0
## [10] S4Vectors_0.32.0 BiocGenerics_0.40.0 BiocStyle_2.22.0
##
## loaded via a namespace (and not attached):
## [1] MatrixGenerics_1.6.0
## [2] httr_1.4.2
## [3] sass_0.4.0
## [4] bit64_4.0.5
## [5] bslib_0.3.1
## [6] assertthat_0.2.1
## [7] BiocManager_1.30.16
## [8] BiocFileCache_2.2.0
## [9] blob_1.2.2
## [10] GenomeInfoDbData_1.2.7
## [11] Rsamtools_2.10.0
## [12] yaml_2.2.1
## [13] progress_1.2.2
## [14] lattice_0.20-45
## [15] pillar_1.6.4
## [16] RSQLite_2.2.8
## [17] glue_1.4.2
## [18] digest_0.6.28
## [19] XVector_0.34.0
## [20] colorspace_2.0-2
## [21] Matrix_1.3-4
## [22] htmltools_0.5.2
## [23] XML_3.99-0.8
## [24] pkgconfig_2.0.3
## [25] biomaRt_2.50.0
## [26] bookdown_0.24
## [27] zlibbioc_1.40.0
## [28] purrr_0.3.4
## [29] scales_1.1.1
## [30] BiocParallel_1.28.0
## [31] tibble_3.1.5
## [32] KEGGREST_1.34.0
## [33] ggplot2_3.3.5
## [34] generics_0.1.1
## [35] ellipsis_0.3.2
## [36] cachem_1.0.6
## [37] SummarizedExperiment_1.24.0
## [38] magrittr_2.0.1
## [39] crayon_1.4.1
## [40] memoise_2.0.0
## [41] evaluate_0.14
## [42] fansi_0.5.0
## [43] xml2_1.3.2
## [44] tools_4.1.1
## [45] prettyunits_1.1.1
## [46] hms_1.1.1
## [47] matrixStats_0.61.0
## [48] BiocIO_1.4.0
## [49] lifecycle_1.0.1
## [50] stringr_1.4.0
## [51] munsell_0.5.0
## [52] DelayedArray_0.20.0
## [53] Biostrings_2.62.0
## [54] compiler_4.1.1
## [55] jquerylib_0.1.4
## [56] rlang_0.4.12
## [57] grid_4.1.1
## [58] RCurl_1.98-1.5
## [59] rjson_0.2.20
## [60] rappdirs_0.3.3
## [61] bitops_1.0-7
## [62] rmarkdown_2.11
## [63] gtable_0.3.0
## [64] restfulr_0.0.13
## [65] DBI_1.1.1
## [66] curl_4.3.2
## [67] R6_2.5.1
## [68] GenomicAlignments_1.30.0
## [69] knitr_1.36
## [70] dplyr_1.0.7
## [71] rtracklayer_1.54.0
## [72] fastmap_1.1.0
## [73] bit_4.0.4
## [74] utf8_1.2.2
## [75] filelock_1.0.2
## [76] stringi_1.7.5
## [77] parallel_4.1.1
## [78] Rcpp_1.0.7
## [79] TxDb.Mmusculus.UCSC.mm9.knownGene_3.2.2
## [80] vctrs_0.3.8
## [81] png_0.1-7
## [82] dbplyr_2.1.1
## [83] tidyselect_1.1.1
## [84] xfun_0.27