knitr::opts_chunk$set(tidy = FALSE, cache = FALSE, dev = "png", message = FALSE, error = FALSE, warning = TRUE, fig.dpi = 96)
library(cinaR) data("atac_seq_consensus_bm")
Bed formatted consensus matrix (chr, start, end and samples)
dim(bed)
# bed formatted file head(bed[,1:4])
Create the contrasts you want to compare, here we create contrasts for 22 mice samples from different strains.
# create contrast vector which will be compared. contrasts<- c("B6", "B6", "B6", "B6", "B6", "NZO", "NZO", "NZO", "NZO", "NZO", "NZO", "B6", "B6", "B6", "B6", "B6", "NZO", "NZO", "NZO", "NZO", "NZO", "NZO")
cinaR
function directly computes the differentially accessible peaks.
# If reference genome is not set hg38 will be used! results <- cinaR(bed, contrasts, reference.genome = "mm10")
Now, you can access differential accessibility (DA) and enrichment results.
names(results)
Inside DA.results
, you have the consensus peaks (cp) and differentially accessible (DA) peaks.
If batch correction was run, then cp
will be a batch-corrected consensus matrix, otherwise it is the filtered and normalized version of original consensus peaks you provided.
names(results$DA.results)
There are many information cinaR
provides such as adjusted p value,
log fold-changes, gene names etc for each peak:
colnames(results$DA.results$DA.peaks$B6_NZO)
Here is an overview of those DA peaks:
head(results$DA.results$DA.peaks$B6_NZO[,1:5])
Since the comparison is
B6_NZO
, if fold-changes are positive it means they are more accesible in B6 compared to NZO and vice versa for negative values!
and here is a little overview for enrichment analyses results:
head(results$Enrichment.Results$B6_NZO[,c("module.name","overlapping.genes", "adj.p")])
You can easily get the PCA plots of the samples:
pca_plot(results, contrasts, show.names = F)
You can overlay different information onto PCA plots as well!
# Overlaid information overlaid.info <- c("B6-18mo", "B6-18mo", "B6-18mo", "B6-18mo", "B6-18mo", "NZO-18mo", "NZO-18mo", "NZO-18mo", "NZO-18mo", "NZO-18mo", "NZO-18mo", "B6-3mo", "B6-3mo", "B6-3mo", "B6-3mo", "B6-3mo", "NZO-3mo", "NZO-3mo", "NZO-3mo", "NZO-3mo", "NZO-3mo", "NZO-3mo") # Sample IDs sample.names <- c("S01783", "S01780", "S01781", "S01778", "S01779", "S03804", "S03805", "S03806", "S03807", "S03808", "S03809", "S04678", "S04679", "S04680", "S04681", "S04682", "S10918", "S10916", "S10919", "S10921", "S10917", "S10920")
pca_plot(results, overlaid.info, sample.names)
You can see the available comparisons using:
show_comparisons(results)
Then, plot the differential peaks for a selected contrast using:
heatmap_differential(results, comparison = "B6_NZO")
Also, you can configure your heatmaps using the additional arguments of pheatmap
function. For more information check out ?pheatmap
.
heatmap_differential(results, comparison = "B6_NZO", show_colnames = FALSE)
You can also plot most variable 100 peaks for all samples:
heatmap_var_peaks(results)
Plus, you can set the number of peaks to be used in these plots, and again you can change the additional arguments of pheatmap
function. For more information check out ?pheatmap
.
heatmap_var_peaks(results, heatmap.peak.count = 200, cluster_cols = F)
You can plot your enrichment results using:
dot_plot(results)
if it gets too crowded, you can get rid of the irrelevant pathways as well:
dot_plot(results, filter.pathways = T)
Note that you can further divide the resolution of contrasts, for instance this is also a valid vector
contrasts <- sapply(strsplit(colnames(bed), split = "-", fixed = T), function(x){paste(x[1:4], collapse = ".")})[4:25] unique(contrasts)
in this case, each of them will be compared to each other which will result in 28 different comparisons.
As default, cinaR
will use one vs one (OVA) scheme, comparing each contrast to others one by one. However, this can be changed to one vs all (OVA) which will compare each contrast to rest:
```{eval=FALSE}
cinaR(..., comparison.scheme = "OVO")
cinaR(..., comparison.scheme = "OVA")
## Running for bulk RNA-seq data To run `cinaR` with RNA-seq experiments, just provide the count matrix, and specify the experiment type: ```r cinaR(matrix = count.matrix, ..., experiment.type = "RNA-Seq")
Note that, count.matrix
should be in the form of $g \times (1+n)$ where $g$ is the number of genes and $n$ is the number of samples, and plus one is for gene names.
Note that currently
cinaR
can only handle gene symbols (e.g. FOSL2, FOXA) and ensembl stable IDs (e.g. ENSG00000010404) for both human and mice!
You can run the enrichment analyses with a custom gene set:
cinaR(..., geneset = new_geneset)
cinaRgenesets
Easiest way to do this is to use cinaRgenesets package. You can select your gene set of interest and just plug it into your pipeline.
You can also download different gene sets from MSigDB site. Note that you should use the human gene symbol versions.
You can use
read.gmt
function fromqusage
package to read these genesets into your current environment.
A geneset
must be a .gmt
formatted symbol file.
You can familiarize yourself with the format by checking out :
# default geneset to be used data("VP2008")
If you have gene and pathway names in a
data.frame
, you can usesplit
function to create your own.gmt
formatted gene sets e.g.split(df$genes, df$pathways)
.
For now, cinaR
supports 3 genomes for human and mice models:
hg38
hg19
mm10
You can set your it using reference.genome
argument.
If you suspect your data have unknown batch effects, you can use:
cinaR(..., batch.correction = T)
This option will run Surrogate Variable Analysis (SVA) and try to
adjust your data for unknown batch effects. If however, you already know the batches of the samples, you can simply set the batch.information
argument as well. This will not run the SVA but just add the batches to design matrix.
# runs SVA cinaR(..., batch.correction = T) # runs SVA with 2 surrogate variables cinaR(..., batch.correction = T, sv.number = 2) # adds only batch information to the design matrix! (does not run SVA) # batch.information should be number a vector where # the length of it equals to the number of samples. cinaR(..., batch.correction = T, batch.information = c(rep(0, 11), rep(1,11)))
Reminder - In our example data we have 22 samples
Sometimes, one might want to add additional covariates to adjust the design matrix further, which affects the down-stream analyses. For instance, ages or sexes of the samples could be additional covariates. To account for those:
# Ages of the samples could be not in biological interests but should be accounted for! cinaR(..., additional.covariates = c((18, 11), (3, 11))) # More than one covariate for instance, sex and age sex.info <- c("M", "F", "M", "F", "F", "F", "F", "F", "M", "M", "M", "F", "F", "M", "M", "M", "F", "F", "M", "M", "F", "M") age.info <- c((18, 11), (3, 11) covs <- data.frame(Sex = sex.info, Age = age.info) cinaR(..., additional.covariates = covs)
Setting save.DA.peaks = TRUE
in cinaR
function will create a DApeaks.xlsx
file
in the current directory. This file includes all the comparisons in different tabs.
Additionally, you can set the path/name of the file using
DA.peaks.path
argument after setting save.DA.peaks = TRUE
.
For instance,
results <- cinaR(..., save.DA.peaks = T, DA.peaks.path = "./Peaks_mice.xlsx")
will create an excel file with name Peaks_mice.xlsx
in the current directory.
Currently, cinaR
supports 4 different algorithms, namely;
If not set, it uses edgeR
for differential analyses.
You can change the used algorithm by simply setting DA.choice
argument.
For more information, ?cinaR
# new FDR threshold for DA peaks results <- cinaR(..., DA.fdr.threshold = 0.1) # filters out pathways results <- cinaR(..., enrichment.FDR.cutoff = 0.1) # does not run enrichment pipeline results <- cinaR(..., run.enrichment = FALSE) # creates the piechart from chIpSeeker package results <- cinaR(..., show.annotation.pie = TRUE) # change cut-off value for dot plots dot_plot(..., fdr.cutoff = 0.05)
Robinson MD, McCarthy DJ, Smyth GK (2010). “edgeR: a Bioconductor package for differential expression analysis of digital gene expression data.” Bioinformatics, 26(1), 139-140. doi: 10.1093/bioinformatics/btp616.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK (2015). “limma powers differential expression analyses for RNA-sequencing and microarray studies.” Nucleic Acids Research, 43(7), e47.
Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15:550. 10.1186/s13059-014-0550-8
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.