In sebastian-gregoricchio/Rseb: An R-package for NGS data managing and visualization
knitr::opts_chunk$set(collapse = TRUE, comment = ">", dev = "svg")
#knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file())
#options(tibble.print_min = 4L, tibble.print_max = 4L)
Introduction {#intro}
The concept behind the Rseb (R-package for Simplified End-to-end data Build-up) is to provide a toolkit that allows the automation of different type of tasks avoiding retyping of code and loss of time. Furthermore, the advantage is that, in most of the cases, the functions are built in R-language making it suitable for all the operating systems
From a more biological point of view, this package simplifies many downstream analyses of high-throughput data that otherwise would require many hours of code typing and case-to-case adaptation. Moreover, most of the functions aimed to visualize these kind of data are thought to provide a high level of possible customization with a large number of graphical parameters compared to the commonly used already available tools. Another advantage of this package is that it offers multiple methods, with a corresponding visualization, to quantify the difference of signal between samples, in a qualitatively and/or quantitatively way, without any additional coding.
The guide is divided in four parts: 1) the first one will explore some analyses and visualization of RNA-seq data; 2) the second one the representation and quantification of targeted sequencing data (ChIP-seq, ATAC-seq, etc.); 3) the third part is focused on the analyses of RNA expression by RT-qPCR; 4) while the last part is focused on some of the “general” tools available in this package.
Citation
If you use this package, please cite:
**Citation**
"HDAC1 and PRC2 mediate combinatorial control in SPI1/PU.1-dependent gene repression in murine erythroleukaemia.".
Gregoricchio S. *et al.*, *Nucleic Acids Research* (2022).
*doi*: [10.1093/nar/gkac613](https://doi.org/10.1093/nar/gkac613)
citation("Rseb")
RNA-seq data {#RNAseq}
A common analysis performed on RNA-seq data is the evaluation of the differentially expressed genes between two different conditions (e.g., untreated vs treated cells).
For this analysis it is common to use the R-package DESeq2^[Love, M.I., Huber, W., Anders, S., "Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2", Genome Biology 15(12):550 (2014)] which returns a table as the following one:
data("RNAseq", package = "Rseb")
RNAseq
data("RNAseq", package = "Rseb")
set.seed(floor(runif(1, 1, 1000)))
n = floor(runif(n = 10, min = 1, max = nrow(RNAseq)))
knitr::kable(RNAseq[n,], row.names = F, caption = "**DESeq2 results table example**")
Differentially expressed genes {#diffExp}
The differential genes are defined depending on the Fold Change (FC) of expression and the adjusted p-value (P~adj~).
The function DEstatus helps in this definition. It takes as input the fold change expression and the p-value adjusted, then a threshold for these two parameters can be set to define four status:
- UP-regulated (FC greater than threshold, significant p-value);
- DOWN-regulated (FC lower than threshold, significant p-value);
- UNRESPONSIVE/NoResp (FC within the range defined by the unresponsive FC threshold, not significant p-value);
- NULL (all the other genes).
All the labels and thresholds can be custom. We can proceed to add a column with the differential expression status to the original table.
require(dplyr)
RNAseq <-
RNAseq %>%
mutate(DE.status = Rseb::DE.status(log2FC = RNAseq$log2FC,
p.value.adjusted = RNAseq$padj,
FC_threshold = 2, # Linear value
FC_NoResp_left = 0.9, # Automatically 0.9 <= FC <= 1/0.9)
p.value_threshold = 0.05,
low.FC.status.label = "DOWN",
high.FC.status.label = "UP",
unresponsive.label = "UNRESPONSIVE",
null.label = "NULL"))
knitr::kable(RNAseq[n,], row.names = F, caption = "**RNA-seq table with differential expression status**")
It is possible now to use the DE.status column to count the number of genes per each group:
RNAseq.summary.table <-
RNAseq %>%
group_by(DE.status) %>%
summarise(N = n()) %>%
rbind(c("Total", nrow(RNAseq)))
knitr::kable(RNAseq.summary.table, row.names = F, caption = "<b>RNA-seq differential expression summary</b>",
format = "html", table.attr = "style='width:30%';")
Representation of the RNA-seq data {#RNAviz}
Two simple representations for RNA-seq data are the MA-plot and the volcano-plot. The first allows to visualize the Fold Change expression as function of basal expression of each gene, while the second always the FC but depending on the significance of the FC computation.
MA-plot {#MAplot}
The MA-plot helps to estimate the difference between two samples plotting the Fold Change expression of a gene as function of its expression among all the samples (all the replicates of both conditions compared, defined by the "baseMean" in the DESeq2 output table). Different colors could be used depending on the "DE.status" column that we just added to the RNA-seq table.
require(ggplot2)
MA.plot <-
ggplot(data = RNAseq,
aes(x = log2(baseMean),
y = log2FC,
col = DE.status)) +
geom_point(size = 2) +
scale_color_manual(values = c("#F8766D", "gray30", "#00A5CF", "#00BA38")) +
ggtitle("MA-plot") +
theme_classic()
MA.plot
Volcano-plot {#volcano}
To plot the Fold Change as function of the p-value it is possible to use the function volcano.
volcano.plot <-
Rseb::volcano(log2FC_data = RNAseq$log2FC,
padj_data = RNAseq$padj,
FC_t = 2,
p_t = 0.05,
FC_unresponsive_left = 0.9,
left_label = "DOWN",
unresponsive_label = "UNRESPONSIVE",
right_label = "UP",
null_label = "NULL",
title = "Volcano",
point_size = 2)
volcano.plot
This function allows to label the points corresponding to differentially expressed genes with the corresponding gene name if required.
volcano.plot.with.names <-
Rseb::volcano(log2FC_data = RNAseq$log2FC,
padj_data = RNAseq$padj,
names = RNAseq$geneName,
FC_t = 2,
p_t = 0.05,
FC_unresponsive_left = 0.9,
left_label = "DOWN",
unresponsive_label = "UNRESPONSIVE",
right_label = "UP",
null_label = "NULL",
title = "Volcano",
point_size = 2,
right_names = T,
names_size = 3)
volcano.plot.with.names
Sequencing signal representation {#sequencingSignal}
Transcription factors and histone marks occupancy (eg. ChIP/Cut&Tag-seq), chromatin condensation (eg. ATAC/MNase-seq), DNA-methylation sites (e.g. MeDIP) are some of the possible high-throughput sequencing signals that could be represented as tracks along the genome by software such as Integrative Genome Browser (IGV) as discussed further. This kind of visualization allows the representation of a single locus at the time, but it could be necessary to summarize the signal of multiple loci simultaneously. To reach this goal, density profiles are the most used method to represent the mean of signal among many loci and compare it for different samples/signals. On the other end, it could be useful to define whether two or more sets of regions show different levels of a signal comparing it region by region.
The general workflow to generate this type of graphs includes the generation of a score matrix for each position of each region and each signal (A). This matrix than can be reshaped in order to calculate the mean/median/sum of the signal among regions in a set (B, density profile) or for each single region (C, density summary. Furthermore, it is possible to compute and visualize the difference of signal between to samples (D, density difference area and/or their correlation). This process is schematized in the following figure and will be discussed in detail in the further paragraphs.
#magick::image_read(rsvg::rsvg_svg("url_image"))
knitr::include_graphics("https://sebastian-gregoricchio.github.io/Rseb/vignettes/images/Rseb_workflow.svg")
Score matrix computation {#matrix}
For the generation of the score matrix I refer to the function computeMatrix from deepTools^[Ramírez, Fidel, Devon P. Ryan, Björn Grüning, Vivek Bhardwaj, Fabian Kilpert, Andreas S. Richter, Steffen Heyne, Friederike Dündar, and Thomas Manke. deepTools2: A next Generation Web Server for Deep-Sequencing Data Analysis. Nucleic Acids Research (2016). doi:10.1093/nar/gkw257.]. This package provides a function called computeMatrix.deeptools that allows to automatically build the command line required to run the tool on the command line.
Alternatively, but based on a more time consuming algorithm, a fully R-based function, density.matrix, is available with this package (see box below). This function generates a list that can be directly used as input for others functions as described further. The parameters are more or less the same found in the computeMatrix.deeptools function and for this reason only the latter will be described in this guide.
**The `density.matrix` function**
This function is completely built in R-code but requires much more time for the execution. This drawback is to attribute to the way by which R handles .bigWig and .bed files. Indeed, a first step required to build a score matrix is to get the bigWig score for each position in each region. This process requires to import a bed file (either from a data.frame or a text file) as a `GRanges` object by the `rtracklayer` package (Bioconductor) and get the score position-by-position from a given bigWig. Unfortunately, `GRanges` objects are quite heavy to handle in R and therefore much more CPU performance and RAM memory are required.
Then, if the chosen modality of matrix computation is “reference-point” all the final regions will have the same length (defined by the upstream and downstream region) and it will be sufficient to split each region score vector (with length n equal to the number of bases of a given region) in chunks of length equal to the required bin size to then calculate the mean/median/sum for each chuck. This process will result in a final score vector of length n/binSize.
In the case the mode selected is instead “scale-regions” the process is more complicated. Indeed, to scale the regions in a perfect way it would require bringing every single region to the body size indicated by the users. This would require calculating the least common multiple (lcm) between the lengths of all the regions and the indicated body length and repeat each base score as much times as lcm/lengthRegion (“base-repetition”). The issue is the size order of the lcm number that makes quite hard the computation in a reasonable time. Hence, the escamotage used in the function is to bring, with the “base-repetition” method, each region to a length ≥ body length and then split in an almost-equal number of chunks. The number of chunks is equivalent to body length. This definition ends often in up to ~5 bases excluded by this splitting in chunks which will be integrated in the last chuck generated be-fore to compute the mean of the scores inside each chunk. At the end of this process to each region corresponds a list of scores of length equivalent to the user-defined body length. Now that all the regions have the same length it is possible to proceed to define the bins scores as described for the “reference-point” mode.
To generate the score matrix, the required arguments are the score/signal files (usually .bigWig/.bw files) and regions files (.bed files).
The function `build.bed` guides the user in the creation of bed files with the option to include header ([track line](https://genome.ucsc.edu/goldenPath/help/hgTracksHelp.html#TRACK)) useful for visualization on IGV or UCSC genome browser (eg. color shadow by region score, track color, track name, etc.). This function allows to export directly the bed file as well save it into a variable as data.frame. Furthermore, many internal controls are automatically performed such as check that END position comes after START position or that all the columns are filled if not provided by the user. Beside the quality controls, the bed could be automatically sorted depending on the genomic position and score^[The "classic" way to sort a bed file is to sort it by chromosome > start > end, but sometimes it is required a more extended sorting by score. Missed score sorting could lead to encounter errors in certain software.] (by the function `sort.bed` of this package) and duplicated rows are removed.
wzxhzdk:12
wzxhzdk:13
Sometimes it happens that in a bed file are present overlapping or adjacent regions that we would like to merge in order to obtain an unique one composed by the fusion of the original regions (*e.g.*, adjacent ChIP-seq peaks). To do that, in this package is available the `collapse.bed` function. This function allows to specify whether the merge should take into account the region strandness (the merge will be applied only if the regions are in the same strand), or to merge the regions in only one strand. Furthermore, the output file consists in a .bed^[Only up to 6 columns (BED6) will be kept: chr, start, end, name, score, strand] in which the region names will correspond to the concatenation of the names of the original regions that served to produced the overlapped region, and the scores (if available) will correspond to the median/mean/sum of the parental scores.
Once built and/or available the region and signal files it is possible to generate the score matrix. The key parameters are^[For more details on parameters visit [computeMatrix](https://deeptools.readthedocs.io/en/develop/content/tools/computeMatrix.html) manual page.]:
* `mode`:
- *"reference-point"* to center all the regions on a specific position
- *"scale-regions"* to re-size all the regions in order to make correspond the extremities of the regions;
* `scoreFileName`: a vector with the full path to the signal file(s);
* `regionsFileName`: a vector with the full path to the regions file(s);
* `outFileName`: the full path for the output matrix archive (.gz);
* `upstream`|`downstream`: to indicate the number of bases to consider before (upstream)/after (downstream) the reference point/region extremities depending on the `mode` used;
* `binSize`: the window size by which each region will be subdivided and for each of which the scores will be calculate;
* `computeMatrix.deeptools.command`: the path to the computeMatrix script (here an example for conda installed deepTools).
wzxhzdk:14
Score matrix from deepTools can be used directly (full path) to generate density plots or can be read by the function `read.computeMatrix.file` which returns a list composed of a data.frame containing the matrix metadata (bin size, regions size, sample_names, region_names, etc.) and the scores table. To plot the signals it can be used this list as well as the path of the source matrix.gz archive.
wzxhzdk:15
An example of deepTools matrix is available with the package:
wzxhzdk:16
wzxhzdk:17
wzxhzdk:18
wzxhzdk:19
wzxhzdk:20
Density profile {#densityProfile}
High-throughput sequencing data allow the detection of different types of signals at genome-wide scale as well at a single locus as described further. A common analysis involves the comparison of the average signal over different categories of regions (e.g., Transcription Staring Site (TSS) vs Enhancers, Activated vs Repressed genes, etc.). For this aim the average signal over regions of the same group can be plotted and compared sample-by-sample or group-by-group. To generate this kind of signal profile it can be used the plot.density.profile/plot.density.profile.smooth function which accepts as input a score matrix as described hereafter. Noteworthy, it is possible to narrow the final plot to a specific range by the option "subset.range" without the need to re-compute the matrix.
The function returns a list with the scores table, a list of each single plot and a single multiplot with the combination of all the single plots. The latter can be directly exported by passing the full output path to the option "multiplot.export.file", while the option "print.multiplot" allows to show the multiplot output without call/print the variable (useful to save time during the parameters test phase).
Here an example of the output list structure:
data("deeptools.matrix", package = "Rseb")
density.profile.by.group <-
Rseb::plot.density.profile(
matrix.file = deeptools.matrix,
signal.type = "mean",
error.type = "sem",
plot.by.group = T,
missing.data.as.zero = T,
y.identical.auto = T,
text.size = 10,
axis.line.width = 0.25,
line.width = 0.5,
plot.error = T,
write.reference.points = T,
plot.vertical.lines = F,
colors = c("steelblue", "mediumseagreen"),
n.row.multiplot = 2,
by.row = T,
y.lab = "Mean ChIP-seq signal \u00b1 SEM")
str(density.profile.by.group, max.level = 1, give.attr = F)
# Load the example matrix list
data("deeptools.matrix", package = "Rseb")
# Generate the density profile by group
density.profile.by.group <-
Rseb::plot.density.profile(
matrix.file = deeptools.matrix,
signal.type = "mean",
error.type = "sem",
plot.by.group = T,
missing.data.as.zero = T,
y.identical.auto = T,
text.size = 10,
axis.line.width = 0.25,
line.width = 0.5,
plot.error = T,
write.reference.points = T,
plot.vertical.lines = F,
colors = c("steelblue", "mediumseagreen"),
n.row.multiplot = 2,
by.row = T,
y.lab = "Mean ChIP-seq signal \u00b1 SEM")
# Generate the density profile by sample
density.profile.by.sample <-
Rseb::plot.density.profile(
matrix.file = deeptools.matrix,
signal.type = "mean",
error.type = "sem",
plot.by.group = FALSE,
missing.data.as.zero = TRUE,
y.identical.auto = TRUE,
text.size = 10,
axis.line.width = 0.25,
line.width = 0.5,
plot.error = TRUE,
write.reference.points = TRUE,
plot.vertical.lines = FALSE,
colors = rep(c("indianred", "darkgoldenrod2"), 2),
line.type = rep(c("solid", "dotted"), each = 2),
n.row.multiplot = 1,
legend.position = c(0.2, 0.75),
y.lab = "Mean ChIP-seq signal \u00b1 SEM")
cowplot::plot_grid(density.profile.by.group$multiplot,
density.profile.by.sample$multiplot,
nrow = 2, labels = "AUTO", rel_heights = c(2, 1))
Quantification of the signal differences
With the density profiles it is possible to qualitatively estimate the difference of signal between two or more groups of regions, or samples, over the same region type. Albeit sometimes the difference between two datasets is obvious, other cases might require a more precises quantification and statistical evaluation of the difference.
Therefore, below in this section, three different methods to infer and visualize the difference between two samples will be explored: the density summary, the area difference and the correlation plot.
All these methods import a score matrix generated, for instance, by deepTools or by the density.matrix function.
Noteworthy, in all cases it is possible to carryout the analyses narrowing down the range on which compute the scores by specifying the range to use in the "subset.range" parameter without the need to re-compute the matrix.
Statistics of the signal {#densityStat}
Both the two analyses that will be addressed in the next paragraphs, beside the graphics, provide a table with the p-values for the comparisons. Indeed, two types of comparisons could be carried out:
- region-by-region: in this case each score of one signal over a single region is compared with the score of an other signal on the same region in a paired-test.
- group-by-group/sample-by-sample: in this case all the scores on all the regions for one group/sample can be compared with those of an other group/sample in an un-paired test.
For both this two methods a t-test (parametric) or Wilcoxon-test (non parametric) could be automatically performed, whereas a test on the distributions cold be manually performed by Kruskal-Wallis test using the values stored in the data.table in the output.
Density summary {#densitySummary}
The first type of representation is generated by the function plot.density.summary which shows, through a violin plot, the distribution of the scores over each region in a given group/sample. An option allows to add a mean symbol with the corresponding SD/SEM. Furthermore, the statistical comparison could be added directly above each violin plot either as "stars" or numeric format. As well as for a boxplot, on each violin plot three lines indicate respectively the 25^th^, 50^th^ and 75^th^ percentile of the data distribution, respectively.
Since the original size of the regions could not be the same it, is preferable to use as score on which compute the analyses the sum of the score over each region (equivalent to the area under the curve of a density profile).
The output of this function is composed of a list of single plots, as well as a multiple plot (multiplot) with the distribution of the scores of the different samples per each group or alternatively the scores on different regions per each sample or group. The output list provides a summary plot with all the region groups for each sample, or alternatively all the samples for each region group. Finally, also a table with all the comparisons with the respective significance is available.
Here an example of the output list structure:
data("deeptools.matrix", package = "Rseb")
summary.plot.by.group <-
Rseb::plot.density.summary(
matrix.file = deeptools.matrix,
plot.by.group = TRUE,
missing.data.as.zero = TRUE,
signal.type = "sum",
stat.paired = FALSE,
stat.hide.ns = FALSE,
axis.line.width = 0.5,
mean.symbol.size = 0.2,
y.identical.auto = TRUE,
text.size = 10,
border.width = 0.25,
subset.range = c(-1000, 1000), #bp from 0 point
colors = c("steelblue", "mediumseagreen"),
legend.position = "right",
n.row.multiplot = 2,
by.row = TRUE)
str(summary.plot.by.group, max.level = 1, give.attr = F)
# Load the example matrix list
data("deeptools.matrix", package = "Rseb")
# Generate the density profile by group
summary.plot.by.group <-
Rseb::plot.density.summary(
matrix.file = deeptools.matrix,
plot.by.group = TRUE,
missing.data.as.zero = TRUE,
signal.type = "sum",
stat.paired = FALSE,
stat.hide.ns = FALSE,
axis.line.width = 0.5,
mean.symbol.size = 0.2,
y.identical.auto = TRUE,
text.size = 10,
border.width = 0.25,
subset.range = c(-1000, 1000), #bp from 0 point
colors = c("steelblue", "mediumseagreen"),
legend.position = "right",
n.row.multiplot = 2,
by.row = TRUE)
# Generate the density profile by sample
summary.plot.by.sample <-
Rseb::plot.density.summary(
matrix.file = deeptools.matrix,
plot.by.group = FALSE,
missing.data.as.zero = TRUE,
signal.type = "sum",
stat.paired = FALSE,
stat.hide.ns = FALSE,
axis.line.width = 0.5,
mean.symbol.size = 0.2,
y.identical.auto = TRUE,
text.size = 10,
border.width = 0.25,
subset.range = c(-1000, 1000), #bp from 0 point
colors = c("steelblue", "mediumseagreen"),
legend.position = "right",
n.row.multiplot = 2,
by.row = TRUE)
cowplot::plot_grid(summary.plot.by.group$multiplot,
summary.plot.by.sample$multiplot,
nrow = 1, labels = "AUTO", rel_widths = c(3, 2))
cowplot::plot_grid(summary.plot.by.group$summary.plot.samples,
summary.plot.by.group$summary.plot.regions,
nrow = 2, labels = "AUTO")
summary.plot.by.sample$means.comparisons
knitr::kable(summary.plot.by.sample$means.comparisons, row.names = F, caption = "**Example of comparisons table for the plots by sample**")
Density differences {#densityDifferences}
In this section two methods to represent the difference of signal between two conditions region-by-region will be described. Both methods, that involve the generation of a correlation/scatter plot and an area/line plot, can be computed by the function plot.density.differences. Also in this case the input of the function is a score matrix generated by computeMatrix.deeptools or density.matrix and it is possible to scale-down to a specific range the final plot by the option "subset.range" without the need to re-compute the matrix..
Signal correlation {#correlationPlot}
On method by which it is possible to visualize the difference between two samples is to correlate the signal for each given region in the two conditions, generating in this way a correlation/scatter plot. Dots corresponding to regions with no difference between two conditions will lie on the diagonal line (defined by the equation $y = x$), while regions with a signal difference will be placed in the half-quadrant of one or the other axis (the two conditions) depending on which condition has a stronger signal on the given region. To generate this kind of plots it is possible to use the function as described below.
# Load the example matrix list
data("deeptools.matrix", package = "Rseb")
# Generate the correlation.plot
correlation.plot <-
Rseb::plot.density.differences(matrix.file = deeptools.matrix,
inverted.comparisons = T,
missing.data.as.zero = T,
signal.type = "sum",
stat.paired = T,
points.size = 0.25,
axis.line.width = 0.25,
area.y.identical.auto = T,
text.size = 10,
correlation.show.equation = T,
correlation.correlation.line.width = 0.5,
correlation.correlation.line.color = "#ED8141",
subset.range = c(-1000, 1000), #bp from 0 point
colors = c("steelblue", "mediumseagreen", "purple"),
legend.position = c(0.8, 0.25),
error.type = "sem")
cowplot::plot_grid(plotlist = Rseb::uniform.y.axis(Rseb::uniform.x.axis(correlation.plot$correlation.plot.byGroup.list)),
nrow = 2, byrow = T)
Signal difference {#areaPlot}
As already mentioned above, another way to quantify the difference between two signals is literally the "difference", intended as subtraction, between the signal over a region (equivalent to the area under the curve in a density profile but for only one region) in one condition and another. This could be computed region-by-region by the function plot.density.differences, but in this case instead of the population shift the regions are ranked by the value of the area/signal subtraction. Then the closest to 0 ("no difference") rank position is calculated and its value subtracted to all the rank positions. In this way on the plot it is possible to know the number of regions with a negative (lower signal) or positive (higher signal) score. This representation allows firstly to know how many regions have an enriched signal (Sample 1) compared to those that have an enrichment of the other compared signal (Sample 2), secondly setting the axis to the same ranges it is possible to observe whether de differences between two samples for a region set compared to another are wide or rather narrowed around the 0 (no difference between the samples signal), beside eventual outliers that could distort the mean in the density profile.
# Load the example matrix list
data("deeptools.matrix", package = "Rseb")
# Generate the area.plot
areaDifference.plot <-
Rseb::plot.density.differences(matrix.file = deeptools.matrix,
inverted.comparisons = T,
missing.data.as.zero = T,
signal.type = "sum",
stat.paired = T,
points.size = 0.25,
axis.line.width = 0.25,
area.y.identical.auto = T,
text.size = 10,
subset.range = c(-1000, 1000), #bp from 0 point
colors = c("steelblue", "mediumseagreen", "purple"),
legend.position = c(0.2, 0.80),
error.type = "sem")
cowplot::plot_grid(plotlist = Rseb::uniform.y.axis(areaDifference.plot$area.plot.byGroup.list),
nrow = 2, byrow = T)
IGV script {#IGV}
The previous paragraphs described how to visualize the signal, or the difference of signal, over multiple regions. In this section instead I will focus on the generation of a script that can be run in the Integrative Genome Browser (IGV^[For more details on parameters visit the IGV batch parameters page.]) software in order to generate a large number of snapshots of multiple loci for the desired signal tracks (.bw files) and/or regions (.bed files) automatically.
The function IGVsnap accepts as input both a list of gene symbols (the genomic localization will be retrieved using biomaRt) or a list of loci in any format accepted by IGV. As output will be generated a text file with a script that can be run on IGV: IGV > Tools > Run Batch Script... > chose your batch script file generated by IGVsnap.
To chose the tracks that will be used there are two possibilities:
- specify an IGV session file (.xml) to be opened before the generation of the snapshots automatically when the script is run;
- run the batch script only after have opened all the required tracks on IGV.
# Define the gene list
gene_list <- c("Spi1", "Idh1", "Bcl2l11", "Mcl1", "Polr2a", "Hdac1")
# Generate the script
Rseb::IGVsnap(loci_vector = gene_list,
input_type = "genes",
biomart = "ensembl",
dataset = "mmusculus_gene_ensembl",
reference_genome = "mm10",
fivePrime = 1000,
threePrime = 1000,
snap_names = gene_list,
IGV_batch_file = "path/to/output/script.txt",
snap_directory = "path/to/directory/in/wich/the/images/will/be/generated",
snap_image_format = "png",
maxPanelHeight = 1500)
Real-Time qPCR analyses for RNA expression {#qPCR}
This tool is designed to analyse RNA expression RT-qPCR data. In particular it is based on the output of RT-qPCR machines of the Applied Biosystems. Indeed it is possible to use the function qPCR.rna.exp directly on results in excel file format generated by the qPCR devices (it could be required to indicate the position number of the excel sheet to use, using the parameter results.sheet.position). However, it is also possible to provide a data.frame with at least two columns: Target Name (e.g., primers), Sample Name (e.g. samples), CT [all must include space and capital letters as indicated].
A mandatory parameter is the list of "housekeeping" genes that will be used for the computation of the normalized expression. On the other hand, also fold changes (FC) over a reference sample are computed. The reference sample if it is not provided by the user will be automatically defined as the first one appearing in the list of all the samples. Notably, specific samples (e.g., blank the negative controls) could be excluded using the parameter exclude.samples.
In the list resulting form the analyses there is a plot showing the technical replicates detected in the data. Technical replicates that as a difference greater than a certain threshold could be excluded setting the parameter max.delta.reps (default value 0.5)
How to run the analyses
# Load the data
data_to_analyse <- Rseb::qPCR.results.rep1
print(data_to_analyse)
data_to_analyse = data.frame(Rseb::qPCR.results.rep1, check.names = F)
knitr::kable(head(data_to_analyse, n = 10), row.names = F, caption = '**Example of RT-qPCR data.frame for RNA expression analyses**')
# Run the analyses
analyses <- Rseb::qPCR.rna.exp(results.file = data_to_analyse,
housekeeping.genes = c("geneB", "housekeeping"),
reference.sample = "Ctrl",
max.delta.reps = 0.5)
Interpretation of the output
The structure of the analyzed data is the following:
cat(paste0(
" analyses\n",
" ¦--original.table\n",
" ¦--reshaped.table\n",
" ¦--reshaped.table.cleaned\n",
" ¦--reps.validation.plot\n",
" ¦--analyzed.data\n",
" ¦ ¦--geneB\n",
" ¦ ¦--housekeeping\n",
" ¦ °--mean_FC_housekeeping\n",
" ¦--expression.plots\n",
" ¦ ¦--geneB\n",
" ¦ °--housekeeping\n",
" °--foldChange.plots\n",
" ¦--geneB\n",
" ¦--housekeeping\n",
" °--mean_FC_housekeeping"))
The `original.table` is the input table, while the `reshaped.table` contains the replicates and their deltas. The `reshaped.table.cleaned` table instead is a filtered version without the discarded replicates.
knitr::kable(head(analyses$reshaped.table.cleaned, n = 15), row.names = F, caption = '**RT-qPCR analyses results: `reshaped.table.cleaned`**')
The replicate deltas can be inspected in the `reps.validation.plot`. The red background indicates that the corresponding replicate comparison did not pass the threshold test defined by `max.delta.reps`.
print(analyses$reps.validation.plot)
The `analyzed.data` branch is a list of tables, one for each housekeeping and an extra one for the mean of all housekeeping genes. Each table shows the replicate CT values, the replicate deltas, the mean of the (filtered) CTs, the calculated expression normalized to the corresponding houskeeping gene and the Fold Change (FC) of expression over the reference.sample.
knitr::kable(head(analyses$analyzed.data$housekeeping, n = 15), row.names = F, caption = '**RT-qPCR analyses results: `analyzed.data$housekeeping`**')
The results then are plotted either as normalized expression to each housekeeping gene (**A**) or the Fold Change expression of the `reference.sample` (**B**). Furthermore, also a mean of the fold changes obtained from the data normalized to each housekeeping genes is plotted as well (**C**).
cowplot::plot_grid(analyses$expression.plots$housekeeping,
cowplot::plot_grid(analyses$foldChange.plots$housekeeping,
analyses$foldChange.plots$mean.FC.housekeeping,
labels = c("B", "C")),
labels = c("A",""), nrow = 2)
Merge multiple experiments
In most of the cases, it is required to compute the average among multiple experiments. To achieve that, a list of qPCR.ran.exp objects can be submitted to the function qPCR.rna.mean.reps.
The function will automatically define which are the housekeeping genes and the reference sample using the list names or the column names of the tables in the qPCR.ran.exp objects. However, if the reference sample used are multiple and/or not shared among replicates, the first sample in the order is used as reference. The reference sample can be also redefined manually during the averaging: all the Fold Changes will be recomputed. In all cases a message indicating which reference sample has been used is printed at each time (it will be always the first sample in the fold change plots).
Also in this case, specific samples can be excluded from the analyses (e.g., blank negative control).
In the following example two different analyses are run starting from the same data (rep1 and rep2). Then the two replicates are merged by qPCR.rna.mean.reps.
# Run the analyses for the two replicates
rep1 <- Rseb::qPCR.rna.exp(results.file = Rseb::qPCR.results.rep1,
housekeeping.genes = c("geneB", "housekeeping"),
reference.sample = "Ctrl",
max.delta.reps = 0.5)
rep2 <- Rseb::qPCR.rna.exp(results.file = Rseb::qPCR.results.rep2,
housekeeping.genes = "housekeeping",
reference.sample = "sample_1",
max.delta.reps = 0.3)
# Compute the average of the two replicates
mean_reps <- Rseb::qPCR.rna.mean.reps(reps.list = list(rep1, rep2))
The structure of the output is similar to that of [`qPCR.ran.exp`](#qPCR):
cat(paste0(
" mean_reps\n",
" ¦--mean.reps.data.table\n",
" ¦ ¦--geneB\n",
" ¦ ¦--housekeeping\n",
" ¦ °--mean_FC.means\n",
" ¦--mean.exp.plots\n",
" ¦ ¦--geneB\n",
" ¦ °--housekeeping\n",
" °--mean.reps.FC.plots\n",
" ¦--geneB\n",
" ¦--housekeeping\n",
" °--mean_FC.means"))
However the analyzed data tables (`mean.reps.data.table`) are simplified showing only the means (exp = normalized expression, FC = Fold Change) and statistics without the single values of each CT in each replicate:
knitr::kable(head(mean_reps$mean.reps.data.table$housekeeping, n = 15), row.names = F, caption = '**RT-qPCR merged replicates: `mean.reps.data.table$housekeeping`**')
The plots are similar as well, but now also the number of values in each bar of gra[h are indicated since the number of samples and/or housekeeping genes could be different among the replicates:
cowplot::plot_grid(mean_reps$mean.reps.exp.plots$housekeeping,
mean_reps$mean.reps.FC.plots$mean_FC.means,
labels = c("AUTO"), nrow = 1)
Generic tools
This package provides functions that can simplify routine task such as axis equal re-scale among a list of plots or filter lines of a data.frame in a shell-like method, and many others.
Uniform plot axes {#uniformAxis}
In the previous paragraphs the functions uniform.x.axis and uniform.y.axis have been used to re-format and uniform the axis among a list of plots. These functions allow to re-set the maximum, minimum, ticks interval and number formatting of the plot axes.
uniformed.list.of.plots <- Rseb::uniform.x.axis(plot.list = list.of.plots,
x.min = 0,
x.max = NA,
ticks.each = 0.5,
digits = 1)
For this function could be useful to combine multiple lists of plots in a unique one using the function combine.lists of this package. The latter, from an input list of list, generates an output list result of the combination of all the lists contained in the provided input.
Filter data.frame rows by a specific pattern {#grepL}
The function grep in the shell allows to filter all the rows/lines of a table/file containing a specific pattern in any position. The function grepl.data.frame from this package mimics the behavior of the shell function in R and applies it to a data.frame. It returns a logic vector of length equal to the rows number of the data.frame in which TRUE indicates that the pattern have been found in the corresponding line. This vector than could be used to filter the data.frame rows such as using dplyr::filter() function.
In the following example is shown a dummy table of copy number variation (CNV) individuated in a cohort of patients for a list of genes. For a given gene in a certain patient if a CNV has been found it will be indicated the type of the CNV ("gain", "loss", "gain/loss") otherwise an NA is used.
# Loading the table
data("CNV.data", package = "Rseb")
print(CNV.data)
data("CNV.data", package = "Rseb")
knitr::kable(CNV.data[round(runif(n = 10, min = 1, max = nrow(CNV.data))),], row.names = F, caption = "**Example of CNV annotation table**")
Now it is possible to filter only the lines containing, in any position, the pattern "gain/loss" using the function grepl.data.frame.
CNV.gain.loss <- dplyr::filter(.data = CNV.data,
Rseb::grepl.data.frame(data.frame = CNV.data,
pattern = "gain/loss"))
print(CNV.gain.loss)
CNV.gain.loss = dplyr::filter(.data = CNV.data, Rseb::grepl.data.frame(data.frame = CNV.data, pattern = "gain/loss"))
knitr::kable(CNV.gain.loss, row.names = F, caption = '**CNV annotation table filtered for the pattern "gain/loss"**')
Package information {#info}
Documentation
With the package a detailed PDF manual with details for each function and respective parameters is available.
The R-package have been published on GitHub and a git-pages website is available as well. At both these sites it is possible to find the installation procedure, required dependencies, and the links for changeLog, manual and vignette.
Package history and releases
A list of all releases and respective description of changes applied could be found here.
Contact
For any suggestion, bug fixing, commentary please contact:
Sebastian Gregoricchio sebastian.gregoricchio@gmail.com
License
This package is under a GNU General Public License (version 3).
sebastian-gregoricchio/Rseb documentation built on June 14, 2024, 12:22 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(collapse = TRUE, comment = ">", dev = "svg") #knitr::opts_knit$set(root.dir = rprojroot::find_rstudio_root_file()) #options(tibble.print_min = 4L, tibble.print_max = 4L)
Introduction {#intro}
The concept behind the Rseb (R-package for Simplified End-to-end data Build-up) is to provide a toolkit that allows the automation of different type of tasks avoiding retyping of code and loss of time. Furthermore, the advantage is that, in most of the cases, the functions are built in R-language making it suitable for all the operating systems
From a more biological point of view, this package simplifies many downstream analyses of high-throughput data that otherwise would require many hours of code typing and case-to-case adaptation. Moreover, most of the functions aimed to visualize these kind of data are thought to provide a high level of possible customization with a large number of graphical parameters compared to the commonly used already available tools. Another advantage of this package is that it offers multiple methods, with a corresponding visualization, to quantify the difference of signal between samples, in a qualitatively and/or quantitatively way, without any additional coding.
The guide is divided in four parts: 1) the first one will explore some analyses and visualization of RNA-seq data; 2) the second one the representation and quantification of targeted sequencing data (ChIP-seq, ATAC-seq, etc.); 3) the third part is focused on the analyses of RNA expression by RT-qPCR; 4) while the last part is focused on some of the “general” tools available in this package.
Citation
If you use this package, please cite:
**Citation**
"HDAC1 and PRC2 mediate combinatorial control in SPI1/PU.1-dependent gene repression in murine erythroleukaemia.".
Gregoricchio S. *et al.*, *Nucleic Acids Research* (2022).
*doi*: [10.1093/nar/gkac613](https://doi.org/10.1093/nar/gkac613)
Gregoricchio S. *et al.*, *Nucleic Acids Research* (2022).
*doi*: [10.1093/nar/gkac613](https://doi.org/10.1093/nar/gkac613)
citation("Rseb")
RNA-seq data {#RNAseq}
A common analysis performed on RNA-seq data is the evaluation of the differentially expressed genes between two different conditions (e.g., untreated vs treated cells).
For this analysis it is common to use the R-package DESeq2^[Love, M.I., Huber, W., Anders, S., "Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2", Genome Biology 15(12):550 (2014)] which returns a table as the following one:
data("RNAseq", package = "Rseb") RNAseq
data("RNAseq", package = "Rseb") set.seed(floor(runif(1, 1, 1000))) n = floor(runif(n = 10, min = 1, max = nrow(RNAseq))) knitr::kable(RNAseq[n,], row.names = F, caption = "**DESeq2 results table example**")
Differentially expressed genes {#diffExp}
The differential genes are defined depending on the Fold Change (FC) of expression and the adjusted p-value (P~adj~).
The function
DEstatus helps in this definition. It takes as input the fold change expression and the p-value adjusted, then a threshold for these two parameters can be set to define four status:
- UP-regulated (FC greater than threshold, significant p-value);
- DOWN-regulated (FC lower than threshold, significant p-value);
- UNRESPONSIVE/NoResp (FC within the range defined by the unresponsive FC threshold, not significant p-value);
- NULL (all the other genes).
All the labels and thresholds can be custom. We can proceed to add a column with the differential expression status to the original table.
require(dplyr) RNAseq <- RNAseq %>% mutate(DE.status = Rseb::DE.status(log2FC = RNAseq$log2FC, p.value.adjusted = RNAseq$padj, FC_threshold = 2, # Linear value FC_NoResp_left = 0.9, # Automatically 0.9 <= FC <= 1/0.9) p.value_threshold = 0.05, low.FC.status.label = "DOWN", high.FC.status.label = "UP", unresponsive.label = "UNRESPONSIVE", null.label = "NULL"))
knitr::kable(RNAseq[n,], row.names = F, caption = "**RNA-seq table with differential expression status**")
It is possible now to use the DE.status column to count the number of genes per each group:
RNAseq.summary.table <- RNAseq %>% group_by(DE.status) %>% summarise(N = n()) %>% rbind(c("Total", nrow(RNAseq)))
knitr::kable(RNAseq.summary.table, row.names = F, caption = "<b>RNA-seq differential expression summary</b>", format = "html", table.attr = "style='width:30%';")
Representation of the RNA-seq data {#RNAviz}
Two simple representations for RNA-seq data are the MA-plot and the volcano-plot. The first allows to visualize the Fold Change expression as function of basal expression of each gene, while the second always the FC but depending on the significance of the FC computation.
MA-plot {#MAplot}
The MA-plot helps to estimate the difference between two samples plotting the Fold Change expression of a gene as function of its expression among all the samples (all the replicates of both conditions compared, defined by the "baseMean" in the DESeq2 output table). Different colors could be used depending on the "DE.status" column that we just added to the RNA-seq table.
require(ggplot2) MA.plot <- ggplot(data = RNAseq, aes(x = log2(baseMean), y = log2FC, col = DE.status)) + geom_point(size = 2) + scale_color_manual(values = c("#F8766D", "gray30", "#00A5CF", "#00BA38")) + ggtitle("MA-plot") + theme_classic() MA.plot
Volcano-plot {#volcano}
To plot the Fold Change as function of the p-value it is possible to use the function
volcano.
volcano.plot <- Rseb::volcano(log2FC_data = RNAseq$log2FC, padj_data = RNAseq$padj, FC_t = 2, p_t = 0.05, FC_unresponsive_left = 0.9, left_label = "DOWN", unresponsive_label = "UNRESPONSIVE", right_label = "UP", null_label = "NULL", title = "Volcano", point_size = 2) volcano.plot
This function allows to label the points corresponding to differentially expressed genes with the corresponding gene name if required.
volcano.plot.with.names <- Rseb::volcano(log2FC_data = RNAseq$log2FC, padj_data = RNAseq$padj, names = RNAseq$geneName, FC_t = 2, p_t = 0.05, FC_unresponsive_left = 0.9, left_label = "DOWN", unresponsive_label = "UNRESPONSIVE", right_label = "UP", null_label = "NULL", title = "Volcano", point_size = 2, right_names = T, names_size = 3) volcano.plot.with.names
Sequencing signal representation {#sequencingSignal}
Transcription factors and histone marks occupancy (eg. ChIP/Cut&Tag-seq), chromatin condensation (eg. ATAC/MNase-seq), DNA-methylation sites (e.g. MeDIP) are some of the possible high-throughput sequencing signals that could be represented as tracks along the genome by software such as Integrative Genome Browser (IGV) as discussed further. This kind of visualization allows the representation of a single locus at the time, but it could be necessary to summarize the signal of multiple loci simultaneously. To reach this goal, density profiles are the most used method to represent the mean of signal among many loci and compare it for different samples/signals. On the other end, it could be useful to define whether two or more sets of regions show different levels of a signal comparing it region by region.
The general workflow to generate this type of graphs includes the generation of a score matrix for each position of each region and each signal (A). This matrix than can be reshaped in order to calculate the mean/median/sum of the signal among regions in a set (B, density profile) or for each single region (C, density summary. Furthermore, it is possible to compute and visualize the difference of signal between to samples (D, density difference area and/or their correlation). This process is schematized in the following figure and will be discussed in detail in the further paragraphs.
#magick::image_read(rsvg::rsvg_svg("url_image")) knitr::include_graphics("https://sebastian-gregoricchio.github.io/Rseb/vignettes/images/Rseb_workflow.svg")
Score matrix computation {#matrix}
For the generation of the score matrix I refer to the function
computeMatrix from deepTools^[Ramírez, Fidel, Devon P. Ryan, Björn Grüning, Vivek Bhardwaj, Fabian Kilpert, Andreas S. Richter, Steffen Heyne, Friederike Dündar, and Thomas Manke. deepTools2: A next Generation Web Server for Deep-Sequencing Data Analysis. Nucleic Acids Research (2016). doi:10.1093/nar/gkw257.]. This package provides a function called computeMatrix.deeptools that allows to automatically build the command line required to run the tool on the command line.
Alternatively, but based on a more time consuming algorithm, a fully R-based function, density.matrix, is available with this package (see box below). This function generates a list that can be directly used as input for others functions as described further. The parameters are more or less the same found in the computeMatrix.deeptools function and for this reason only the latter will be described in this guide.
**The `density.matrix` function**
This function is completely built in R-code but requires much more time for the execution. This drawback is to attribute to the way by which R handles .bigWig and .bed files. Indeed, a first step required to build a score matrix is to get the bigWig score for each position in each region. This process requires to import a bed file (either from a data.frame or a text file) as a `GRanges` object by the `rtracklayer` package (Bioconductor) and get the score position-by-position from a given bigWig. Unfortunately, `GRanges` objects are quite heavy to handle in R and therefore much more CPU performance and RAM memory are required.
Then, if the chosen modality of matrix computation is “reference-point” all the final regions will have the same length (defined by the upstream and downstream region) and it will be sufficient to split each region score vector (with length n equal to the number of bases of a given region) in chunks of length equal to the required bin size to then calculate the mean/median/sum for each chuck. This process will result in a final score vector of length n/binSize.
In the case the mode selected is instead “scale-regions” the process is more complicated. Indeed, to scale the regions in a perfect way it would require bringing every single region to the body size indicated by the users. This would require calculating the least common multiple (lcm) between the lengths of all the regions and the indicated body length and repeat each base score as much times as lcm/lengthRegion (“base-repetition”). The issue is the size order of the lcm number that makes quite hard the computation in a reasonable time. Hence, the escamotage used in the function is to bring, with the “base-repetition” method, each region to a length ≥ body length and then split in an almost-equal number of chunks. The number of chunks is equivalent to body length. This definition ends often in up to ~5 bases excluded by this splitting in chunks which will be integrated in the last chuck generated be-fore to compute the mean of the scores inside each chunk. At the end of this process to each region corresponds a list of scores of length equivalent to the user-defined body length. Now that all the regions have the same length it is possible to proceed to define the bins scores as described for the “reference-point” mode.
Then, if the chosen modality of matrix computation is “reference-point” all the final regions will have the same length (defined by the upstream and downstream region) and it will be sufficient to split each region score vector (with length n equal to the number of bases of a given region) in chunks of length equal to the required bin size to then calculate the mean/median/sum for each chuck. This process will result in a final score vector of length n/binSize.
In the case the mode selected is instead “scale-regions” the process is more complicated. Indeed, to scale the regions in a perfect way it would require bringing every single region to the body size indicated by the users. This would require calculating the least common multiple (lcm) between the lengths of all the regions and the indicated body length and repeat each base score as much times as lcm/lengthRegion (“base-repetition”). The issue is the size order of the lcm number that makes quite hard the computation in a reasonable time. Hence, the escamotage used in the function is to bring, with the “base-repetition” method, each region to a length ≥ body length and then split in an almost-equal number of chunks. The number of chunks is equivalent to body length. This definition ends often in up to ~5 bases excluded by this splitting in chunks which will be integrated in the last chuck generated be-fore to compute the mean of the scores inside each chunk. At the end of this process to each region corresponds a list of scores of length equivalent to the user-defined body length. Now that all the regions have the same length it is possible to proceed to define the bins scores as described for the “reference-point” mode.
To generate the score matrix, the required arguments are the score/signal files (usually .bigWig/.bw files) and regions files (.bed files).
The function `build.bed` guides the user in the creation of bed files with the option to include header ([track line](https://genome.ucsc.edu/goldenPath/help/hgTracksHelp.html#TRACK)) useful for visualization on IGV or UCSC genome browser (eg. color shadow by region score, track color, track name, etc.). This function allows to export directly the bed file as well save it into a variable as data.frame. Furthermore, many internal controls are automatically performed such as check that END position comes after START position or that all the columns are filled if not provided by the user. Beside the quality controls, the bed could be automatically sorted depending on the genomic position and score^[The "classic" way to sort a bed file is to sort it by chromosome > start > end, but sometimes it is required a more extended sorting by score. Missed score sorting could lead to encounter errors in certain software.] (by the function `sort.bed` of this package) and duplicated rows are removed.
wzxhzdk:12
wzxhzdk:13
Sometimes it happens that in a bed file are present overlapping or adjacent regions that we would like to merge in order to obtain an unique one composed by the fusion of the original regions (*e.g.*, adjacent ChIP-seq peaks). To do that, in this package is available the `collapse.bed` function. This function allows to specify whether the merge should take into account the region strandness (the merge will be applied only if the regions are in the same strand), or to merge the regions in only one strand. Furthermore, the output file consists in a .bed^[Only up to 6 columns (BED6) will be kept: chr, start, end, name, score, strand] in which the region names will correspond to the concatenation of the names of the original regions that served to produced the overlapped region, and the scores (if available) will correspond to the median/mean/sum of the parental scores.
Once built and/or available the region and signal files it is possible to generate the score matrix. The key parameters are^[For more details on parameters visit [computeMatrix](https://deeptools.readthedocs.io/en/develop/content/tools/computeMatrix.html) manual page.]:
* `mode`:
- *"reference-point"* to center all the regions on a specific position
- *"scale-regions"* to re-size all the regions in order to make correspond the extremities of the regions;
* `scoreFileName`: a vector with the full path to the signal file(s);
* `regionsFileName`: a vector with the full path to the regions file(s);
* `outFileName`: the full path for the output matrix archive (.gz);
* `upstream`|`downstream`: to indicate the number of bases to consider before (upstream)/after (downstream) the reference point/region extremities depending on the `mode` used;
* `binSize`: the window size by which each region will be subdivided and for each of which the scores will be calculate;
* `computeMatrix.deeptools.command`: the path to the computeMatrix script (here an example for conda installed deepTools).
wzxhzdk:14
Score matrix from deepTools can be used directly (full path) to generate density plots or can be read by the function `read.computeMatrix.file` which returns a list composed of a data.frame containing the matrix metadata (bin size, regions size, sample_names, region_names, etc.) and the scores table. To plot the signals it can be used this list as well as the path of the source matrix.gz archive.
wzxhzdk:15
An example of deepTools matrix is available with the package:
wzxhzdk:16
wzxhzdk:17
wzxhzdk:18
wzxhzdk:19
wzxhzdk:20
Density profile {#densityProfile}
High-throughput sequencing data allow the detection of different types of signals at genome-wide scale as well at a single locus as described further. A common analysis involves the comparison of the average signal over different categories of regions (e.g., Transcription Staring Site (TSS) vs Enhancers, Activated vs Repressed genes, etc.). For this aim the average signal over regions of the same group can be plotted and compared sample-by-sample or group-by-group. To generate this kind of signal profile it can be used the
plot.density.profile/plot.density.profile.smooth function which accepts as input a score matrix as described hereafter. Noteworthy, it is possible to narrow the final plot to a specific range by the option "subset.range" without the need to re-compute the matrix.
The function returns a list with the scores table, a list of each single plot and a single multiplot with the combination of all the single plots. The latter can be directly exported by passing the full output path to the option "multiplot.export.file", while the option "print.multiplot" allows to show the multiplot output without call/print the variable (useful to save time during the parameters test phase).
Here an example of the output list structure:
data("deeptools.matrix", package = "Rseb") density.profile.by.group <- Rseb::plot.density.profile( matrix.file = deeptools.matrix, signal.type = "mean", error.type = "sem", plot.by.group = T, missing.data.as.zero = T, y.identical.auto = T, text.size = 10, axis.line.width = 0.25, line.width = 0.5, plot.error = T, write.reference.points = T, plot.vertical.lines = F, colors = c("steelblue", "mediumseagreen"), n.row.multiplot = 2, by.row = T, y.lab = "Mean ChIP-seq signal \u00b1 SEM") str(density.profile.by.group, max.level = 1, give.attr = F)
# Load the example matrix list data("deeptools.matrix", package = "Rseb") # Generate the density profile by group density.profile.by.group <- Rseb::plot.density.profile( matrix.file = deeptools.matrix, signal.type = "mean", error.type = "sem", plot.by.group = T, missing.data.as.zero = T, y.identical.auto = T, text.size = 10, axis.line.width = 0.25, line.width = 0.5, plot.error = T, write.reference.points = T, plot.vertical.lines = F, colors = c("steelblue", "mediumseagreen"), n.row.multiplot = 2, by.row = T, y.lab = "Mean ChIP-seq signal \u00b1 SEM") # Generate the density profile by sample density.profile.by.sample <- Rseb::plot.density.profile( matrix.file = deeptools.matrix, signal.type = "mean", error.type = "sem", plot.by.group = FALSE, missing.data.as.zero = TRUE, y.identical.auto = TRUE, text.size = 10, axis.line.width = 0.25, line.width = 0.5, plot.error = TRUE, write.reference.points = TRUE, plot.vertical.lines = FALSE, colors = rep(c("indianred", "darkgoldenrod2"), 2), line.type = rep(c("solid", "dotted"), each = 2), n.row.multiplot = 1, legend.position = c(0.2, 0.75), y.lab = "Mean ChIP-seq signal \u00b1 SEM") cowplot::plot_grid(density.profile.by.group$multiplot, density.profile.by.sample$multiplot, nrow = 2, labels = "AUTO", rel_heights = c(2, 1))
Quantification of the signal differences
With the density profiles it is possible to qualitatively estimate the difference of signal between two or more groups of regions, or samples, over the same region type. Albeit sometimes the difference between two datasets is obvious, other cases might require a more precises quantification and statistical evaluation of the difference.
Therefore, below in this section, three different methods to infer and visualize the difference between two samples will be explored: the density summary, the area difference and the correlation plot.
All these methods import a score matrix generated, for instance, by deepTools or by the density.matrix function.
Noteworthy, in all cases it is possible to carryout the analyses narrowing down the range on which compute the scores by specifying the range to use in the "subset.range" parameter without the need to re-compute the matrix.
Statistics of the signal {#densityStat}
Both the two analyses that will be addressed in the next paragraphs, beside the graphics, provide a table with the p-values for the comparisons. Indeed, two types of comparisons could be carried out:
- region-by-region: in this case each score of one signal over a single region is compared with the score of an other signal on the same region in a paired-test.
- group-by-group/sample-by-sample: in this case all the scores on all the regions for one group/sample can be compared with those of an other group/sample in an un-paired test.
For both this two methods a t-test (parametric) or Wilcoxon-test (non parametric) could be automatically performed, whereas a test on the distributions cold be manually performed by Kruskal-Wallis test using the values stored in the data.table in the output.
Density summary {#densitySummary}
The first type of representation is generated by the function
plot.density.summary which shows, through a violin plot, the distribution of the scores over each region in a given group/sample. An option allows to add a mean symbol with the corresponding SD/SEM. Furthermore, the statistical comparison could be added directly above each violin plot either as "stars" or numeric format. As well as for a boxplot, on each violin plot three lines indicate respectively the 25^th^, 50^th^ and 75^th^ percentile of the data distribution, respectively.
Since the original size of the regions could not be the same it, is preferable to use as score on which compute the analyses the sum of the score over each region (equivalent to the area under the curve of a density profile).
The output of this function is composed of a list of single plots, as well as a multiple plot (multiplot) with the distribution of the scores of the different samples per each group or alternatively the scores on different regions per each sample or group. The output list provides a summary plot with all the region groups for each sample, or alternatively all the samples for each region group. Finally, also a table with all the comparisons with the respective significance is available.
Here an example of the output list structure:
data("deeptools.matrix", package = "Rseb") summary.plot.by.group <- Rseb::plot.density.summary( matrix.file = deeptools.matrix, plot.by.group = TRUE, missing.data.as.zero = TRUE, signal.type = "sum", stat.paired = FALSE, stat.hide.ns = FALSE, axis.line.width = 0.5, mean.symbol.size = 0.2, y.identical.auto = TRUE, text.size = 10, border.width = 0.25, subset.range = c(-1000, 1000), #bp from 0 point colors = c("steelblue", "mediumseagreen"), legend.position = "right", n.row.multiplot = 2, by.row = TRUE) str(summary.plot.by.group, max.level = 1, give.attr = F)
# Load the example matrix list data("deeptools.matrix", package = "Rseb") # Generate the density profile by group summary.plot.by.group <- Rseb::plot.density.summary( matrix.file = deeptools.matrix, plot.by.group = TRUE, missing.data.as.zero = TRUE, signal.type = "sum", stat.paired = FALSE, stat.hide.ns = FALSE, axis.line.width = 0.5, mean.symbol.size = 0.2, y.identical.auto = TRUE, text.size = 10, border.width = 0.25, subset.range = c(-1000, 1000), #bp from 0 point colors = c("steelblue", "mediumseagreen"), legend.position = "right", n.row.multiplot = 2, by.row = TRUE) # Generate the density profile by sample summary.plot.by.sample <- Rseb::plot.density.summary( matrix.file = deeptools.matrix, plot.by.group = FALSE, missing.data.as.zero = TRUE, signal.type = "sum", stat.paired = FALSE, stat.hide.ns = FALSE, axis.line.width = 0.5, mean.symbol.size = 0.2, y.identical.auto = TRUE, text.size = 10, border.width = 0.25, subset.range = c(-1000, 1000), #bp from 0 point colors = c("steelblue", "mediumseagreen"), legend.position = "right", n.row.multiplot = 2, by.row = TRUE) cowplot::plot_grid(summary.plot.by.group$multiplot, summary.plot.by.sample$multiplot, nrow = 1, labels = "AUTO", rel_widths = c(3, 2))
cowplot::plot_grid(summary.plot.by.group$summary.plot.samples, summary.plot.by.group$summary.plot.regions, nrow = 2, labels = "AUTO")
summary.plot.by.sample$means.comparisons
knitr::kable(summary.plot.by.sample$means.comparisons, row.names = F, caption = "**Example of comparisons table for the plots by sample**")
Density differences {#densityDifferences}
In this section two methods to represent the difference of signal between two conditions region-by-region will be described. Both methods, that involve the generation of a correlation/scatter plot and an area/line plot, can be computed by the function
plot.density.differences. Also in this case the input of the function is a score matrix generated by computeMatrix.deeptools or density.matrix and it is possible to scale-down to a specific range the final plot by the option "subset.range" without the need to re-compute the matrix..
Signal correlation {#correlationPlot}
On method by which it is possible to visualize the difference between two samples is to correlate the signal for each given region in the two conditions, generating in this way a correlation/scatter plot. Dots corresponding to regions with no difference between two conditions will lie on the diagonal line (defined by the equation $y = x$), while regions with a signal difference will be placed in the half-quadrant of one or the other axis (the two conditions) depending on which condition has a stronger signal on the given region. To generate this kind of plots it is possible to use the function as described below.
# Load the example matrix list data("deeptools.matrix", package = "Rseb") # Generate the correlation.plot correlation.plot <- Rseb::plot.density.differences(matrix.file = deeptools.matrix, inverted.comparisons = T, missing.data.as.zero = T, signal.type = "sum", stat.paired = T, points.size = 0.25, axis.line.width = 0.25, area.y.identical.auto = T, text.size = 10, correlation.show.equation = T, correlation.correlation.line.width = 0.5, correlation.correlation.line.color = "#ED8141", subset.range = c(-1000, 1000), #bp from 0 point colors = c("steelblue", "mediumseagreen", "purple"), legend.position = c(0.8, 0.25), error.type = "sem") cowplot::plot_grid(plotlist = Rseb::uniform.y.axis(Rseb::uniform.x.axis(correlation.plot$correlation.plot.byGroup.list)), nrow = 2, byrow = T)
Signal difference {#areaPlot}
As already mentioned above, another way to quantify the difference between two signals is literally the "difference", intended as subtraction, between the signal over a region (equivalent to the area under the curve in a density profile but for only one region) in one condition and another. This could be computed region-by-region by the function
plot.density.differences, but in this case instead of the population shift the regions are ranked by the value of the area/signal subtraction. Then the closest to 0 ("no difference") rank position is calculated and its value subtracted to all the rank positions. In this way on the plot it is possible to know the number of regions with a negative (lower signal) or positive (higher signal) score. This representation allows firstly to know how many regions have an enriched signal (Sample 1) compared to those that have an enrichment of the other compared signal (Sample 2), secondly setting the axis to the same ranges it is possible to observe whether de differences between two samples for a region set compared to another are wide or rather narrowed around the 0 (no difference between the samples signal), beside eventual outliers that could distort the mean in the density profile.
# Load the example matrix list data("deeptools.matrix", package = "Rseb") # Generate the area.plot areaDifference.plot <- Rseb::plot.density.differences(matrix.file = deeptools.matrix, inverted.comparisons = T, missing.data.as.zero = T, signal.type = "sum", stat.paired = T, points.size = 0.25, axis.line.width = 0.25, area.y.identical.auto = T, text.size = 10, subset.range = c(-1000, 1000), #bp from 0 point colors = c("steelblue", "mediumseagreen", "purple"), legend.position = c(0.2, 0.80), error.type = "sem") cowplot::plot_grid(plotlist = Rseb::uniform.y.axis(areaDifference.plot$area.plot.byGroup.list), nrow = 2, byrow = T)
IGV script {#IGV}
The previous paragraphs described how to visualize the signal, or the difference of signal, over multiple regions. In this section instead I will focus on the generation of a script that can be run in the Integrative Genome Browser (IGV^[For more details on parameters visit the IGV batch parameters page.]) software in order to generate a large number of snapshots of multiple loci for the desired signal tracks (.bw files) and/or regions (.bed files) automatically.
The function IGVsnap accepts as input both a list of gene symbols (the genomic localization will be retrieved using biomaRt) or a list of loci in any format accepted by IGV. As output will be generated a text file with a script that can be run on IGV: IGV > Tools > Run Batch Script... > chose your batch script file generated by IGVsnap.
To chose the tracks that will be used there are two possibilities:
- specify an IGV session file (.xml) to be opened before the generation of the snapshots automatically when the script is run;
- run the batch script only after have opened all the required tracks on IGV.
# Define the gene list gene_list <- c("Spi1", "Idh1", "Bcl2l11", "Mcl1", "Polr2a", "Hdac1") # Generate the script Rseb::IGVsnap(loci_vector = gene_list, input_type = "genes", biomart = "ensembl", dataset = "mmusculus_gene_ensembl", reference_genome = "mm10", fivePrime = 1000, threePrime = 1000, snap_names = gene_list, IGV_batch_file = "path/to/output/script.txt", snap_directory = "path/to/directory/in/wich/the/images/will/be/generated", snap_image_format = "png", maxPanelHeight = 1500)
Real-Time qPCR analyses for RNA expression {#qPCR}
This tool is designed to analyse RNA expression RT-qPCR data. In particular it is based on the output of RT-qPCR machines of the Applied Biosystems. Indeed it is possible to use the function
qPCR.rna.exp directly on results in excel file format generated by the qPCR devices (it could be required to indicate the position number of the excel sheet to use, using the parameter results.sheet.position). However, it is also possible to provide a data.frame with at least two columns: Target Name (e.g., primers), Sample Name (e.g. samples), CT [all must include space and capital letters as indicated].
A mandatory parameter is the list of "housekeeping" genes that will be used for the computation of the normalized expression. On the other hand, also fold changes (FC) over a reference sample are computed. The reference sample if it is not provided by the user will be automatically defined as the first one appearing in the list of all the samples. Notably, specific samples (e.g., blank the negative controls) could be excluded using the parameter exclude.samples.
In the list resulting form the analyses there is a plot showing the technical replicates detected in the data. Technical replicates that as a difference greater than a certain threshold could be excluded setting the parameter max.delta.reps (default value 0.5)
How to run the analyses
# Load the data data_to_analyse <- Rseb::qPCR.results.rep1 print(data_to_analyse)
data_to_analyse = data.frame(Rseb::qPCR.results.rep1, check.names = F) knitr::kable(head(data_to_analyse, n = 10), row.names = F, caption = '**Example of RT-qPCR data.frame for RNA expression analyses**')
# Run the analyses analyses <- Rseb::qPCR.rna.exp(results.file = data_to_analyse, housekeeping.genes = c("geneB", "housekeeping"), reference.sample = "Ctrl", max.delta.reps = 0.5)
Interpretation of the output
The structure of the analyzed data is the following:
cat(paste0( " analyses\n", " ¦--original.table\n", " ¦--reshaped.table\n", " ¦--reshaped.table.cleaned\n", " ¦--reps.validation.plot\n", " ¦--analyzed.data\n", " ¦ ¦--geneB\n", " ¦ ¦--housekeeping\n", " ¦ °--mean_FC_housekeeping\n", " ¦--expression.plots\n", " ¦ ¦--geneB\n", " ¦ °--housekeeping\n", " °--foldChange.plots\n", " ¦--geneB\n", " ¦--housekeeping\n", " °--mean_FC_housekeeping"))
The `original.table` is the input table, while the `reshaped.table` contains the replicates and their deltas. The `reshaped.table.cleaned` table instead is a filtered version without the discarded replicates.
knitr::kable(head(analyses$reshaped.table.cleaned, n = 15), row.names = F, caption = '**RT-qPCR analyses results: `reshaped.table.cleaned`**')
The replicate deltas can be inspected in the `reps.validation.plot`. The red background indicates that the corresponding replicate comparison did not pass the threshold test defined by `max.delta.reps`.
print(analyses$reps.validation.plot)
The `analyzed.data` branch is a list of tables, one for each housekeeping and an extra one for the mean of all housekeeping genes. Each table shows the replicate CT values, the replicate deltas, the mean of the (filtered) CTs, the calculated expression normalized to the corresponding houskeeping gene and the Fold Change (FC) of expression over the reference.sample.
knitr::kable(head(analyses$analyzed.data$housekeeping, n = 15), row.names = F, caption = '**RT-qPCR analyses results: `analyzed.data$housekeeping`**')
The results then are plotted either as normalized expression to each housekeeping gene (**A**) or the Fold Change expression of the `reference.sample` (**B**). Furthermore, also a mean of the fold changes obtained from the data normalized to each housekeeping genes is plotted as well (**C**).
cowplot::plot_grid(analyses$expression.plots$housekeeping, cowplot::plot_grid(analyses$foldChange.plots$housekeeping, analyses$foldChange.plots$mean.FC.housekeeping, labels = c("B", "C")), labels = c("A",""), nrow = 2)
Merge multiple experiments
In most of the cases, it is required to compute the average among multiple experiments. To achieve that, a list of
qPCR.ran.exp objects can be submitted to the function qPCR.rna.mean.reps.
The function will automatically define which are the housekeeping genes and the reference sample using the list names or the column names of the tables in the qPCR.ran.exp objects. However, if the reference sample used are multiple and/or not shared among replicates, the first sample in the order is used as reference. The reference sample can be also redefined manually during the averaging: all the Fold Changes will be recomputed. In all cases a message indicating which reference sample has been used is printed at each time (it will be always the first sample in the fold change plots).
Also in this case, specific samples can be excluded from the analyses (e.g., blank negative control).
In the following example two different analyses are run starting from the same data (rep1 and rep2). Then the two replicates are merged by qPCR.rna.mean.reps.
# Run the analyses for the two replicates rep1 <- Rseb::qPCR.rna.exp(results.file = Rseb::qPCR.results.rep1, housekeeping.genes = c("geneB", "housekeeping"), reference.sample = "Ctrl", max.delta.reps = 0.5) rep2 <- Rseb::qPCR.rna.exp(results.file = Rseb::qPCR.results.rep2, housekeeping.genes = "housekeeping", reference.sample = "sample_1", max.delta.reps = 0.3) # Compute the average of the two replicates mean_reps <- Rseb::qPCR.rna.mean.reps(reps.list = list(rep1, rep2))
The structure of the output is similar to that of [`qPCR.ran.exp`](#qPCR):
cat(paste0( " mean_reps\n", " ¦--mean.reps.data.table\n", " ¦ ¦--geneB\n", " ¦ ¦--housekeeping\n", " ¦ °--mean_FC.means\n", " ¦--mean.exp.plots\n", " ¦ ¦--geneB\n", " ¦ °--housekeeping\n", " °--mean.reps.FC.plots\n", " ¦--geneB\n", " ¦--housekeeping\n", " °--mean_FC.means"))
However the analyzed data tables (`mean.reps.data.table`) are simplified showing only the means (exp = normalized expression, FC = Fold Change) and statistics without the single values of each CT in each replicate:
knitr::kable(head(mean_reps$mean.reps.data.table$housekeeping, n = 15), row.names = F, caption = '**RT-qPCR merged replicates: `mean.reps.data.table$housekeeping`**')
The plots are similar as well, but now also the number of values in each bar of gra[h are indicated since the number of samples and/or housekeeping genes could be different among the replicates:
cowplot::plot_grid(mean_reps$mean.reps.exp.plots$housekeeping, mean_reps$mean.reps.FC.plots$mean_FC.means, labels = c("AUTO"), nrow = 1)
Generic tools
This package provides functions that can simplify routine task such as axis equal re-scale among a list of plots or filter lines of a data.frame in a shell-like method, and many others.
Uniform plot axes {#uniformAxis}
In the previous paragraphs the functions
uniform.x.axis and uniform.y.axis have been used to re-format and uniform the axis among a list of plots. These functions allow to re-set the maximum, minimum, ticks interval and number formatting of the plot axes.
uniformed.list.of.plots <- Rseb::uniform.x.axis(plot.list = list.of.plots, x.min = 0, x.max = NA, ticks.each = 0.5, digits = 1)
For this function could be useful to combine multiple lists of plots in a unique one using the function combine.lists of this package. The latter, from an input list of list, generates an output list result of the combination of all the lists contained in the provided input.
Filter data.frame rows by a specific pattern {#grepL}
The function
grep in the shell allows to filter all the rows/lines of a table/file containing a specific pattern in any position. The function grepl.data.frame from this package mimics the behavior of the shell function in R and applies it to a data.frame. It returns a logic vector of length equal to the rows number of the data.frame in which TRUE indicates that the pattern have been found in the corresponding line. This vector than could be used to filter the data.frame rows such as using dplyr::filter() function.
In the following example is shown a dummy table of copy number variation (CNV) individuated in a cohort of patients for a list of genes. For a given gene in a certain patient if a CNV has been found it will be indicated the type of the CNV ("gain", "loss", "gain/loss") otherwise an NA is used.
# Loading the table data("CNV.data", package = "Rseb") print(CNV.data)
data("CNV.data", package = "Rseb") knitr::kable(CNV.data[round(runif(n = 10, min = 1, max = nrow(CNV.data))),], row.names = F, caption = "**Example of CNV annotation table**")
Now it is possible to filter only the lines containing, in any position, the pattern "gain/loss" using the function grepl.data.frame.
CNV.gain.loss <- dplyr::filter(.data = CNV.data, Rseb::grepl.data.frame(data.frame = CNV.data, pattern = "gain/loss")) print(CNV.gain.loss)
CNV.gain.loss = dplyr::filter(.data = CNV.data, Rseb::grepl.data.frame(data.frame = CNV.data, pattern = "gain/loss")) knitr::kable(CNV.gain.loss, row.names = F, caption = '**CNV annotation table filtered for the pattern "gain/loss"**')
Package information {#info}
Documentation
With the package a detailed PDF manual with details for each function and respective parameters is available.
The R-package have been published on GitHub and a git-pages website is available as well. At both these sites it is possible to find the installation procedure, required dependencies, and the links for changeLog, manual and vignette.
Package history and releases
A list of all releases and respective description of changes applied could be found here.
Contact
For any suggestion, bug fixing, commentary please contact:
Sebastian Gregoricchio sebastian.gregoricchio@gmail.com
License
This package is under a GNU General Public License (version 3).
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.