inst/doc/CAGEexp.R

In CAGEr: Analysis of CAGE (Cap Analysis of Gene Expression) sequencing data for precise mapping of transcription start sites and promoterome mining

## ----setup, echo = FALSE, results = "hide"------------------------------------
options(signif = 3, digits = 3)
knitr::opts_chunk$set(tidy = FALSE, cache = TRUE, autodep = TRUE, fig.height = 5.5,
                      message = FALSE, error = FALSE, warning = TRUE)
set.seed(0xdada)

## ----CAGEprotocol, echo=FALSE, fig.align="center", fig.cap="Overview of CAGE experiment"----
knitr::include_graphics("images/CAGEprotocol.svg")

## ----load_CAGEr---------------------------------------------------------------
library(CAGEr)

## ----specify_input_files------------------------------------------------------
inputFiles = list.files( system.file("extdata", package = "CAGEr")
                       , "ctss$"
                       , full.names = TRUE)

## ----create_CAGEexp-----------------------------------------------------------
ce <- CAGEexp( genomeName     = "BSgenome.Drerio.UCSC.danRer7"
             , inputFiles     = inputFiles
             , inputFilesType = "ctss"
             , sampleLabels   = sub( ".chr17.ctss", "", basename(inputFiles))
)

## ----display_created_object---------------------------------------------------
ce

## ----load_data----------------------------------------------------------------
getCTSS(ce)
ce

## -----------------------------------------------------------------------------
CTSStagCountSE(ce)
CTSScoordinatesGR(ce)
CTSStagCountDF(ce)

## -----------------------------------------------------------------------------
head(CTSScoordinates(ce))
head(CTSStagCountDf(ce))
head(CTSStagCount(ce))

## -----------------------------------------------------------------------------
sampleLabels(ce)

## -----------------------------------------------------------------------------
annotateCTSS(ce, exampleZv9_annot)

## -----------------------------------------------------------------------------
colData(ce)[,c("librarySizes", "promoter", "exon", "intron", "unknown")]
CTSScoordinatesGR(ce)

## -----------------------------------------------------------------------------
plotAnnot(ce, "counts")

## ----CorrelationScatterPlots, fig.cap="Correlation of raw CAGE tag counts per TSS"----
corr.m <- plotCorrelation2( ce, samples = "all"
                          , tagCountThreshold = 1, applyThresholdBoth = FALSE
                          , method = "pearson")

## -----------------------------------------------------------------------------
mergeSamples(ce, mergeIndex = c(3,2,4,4,1), 
			mergedSampleLabels = c("zf_unfertilized_egg", "zf_high", "zf_30p_dome", "zf_prim6"))
annotateCTSS(ce, exampleZv9_annot)

## -----------------------------------------------------------------------------
librarySizes(ce)

## ----ReverseCumulatives, fig.cap="Reverse cumulative distribution of CAGE tags"----
plotReverseCumulatives(ce, fitInRange = c(5, 1000), onePlot = TRUE)

## -----------------------------------------------------------------------------
normalizeTagCount(ce, method = "powerLaw", fitInRange = c(5, 1000), alpha = 1.2, T = 5*10^4)
ce[["tagCountMatrix"]]

## -----------------------------------------------------------------------------
clusterCTSS( object = ce
           , threshold = 1
           , thresholdIsTpm = TRUE
           , nrPassThreshold = 1
           , method = "distclu"
           , maxDist = 20
           , removeSingletons = TRUE
           , keepSingletonsAbove = 5)

## -----------------------------------------------------------------------------
head(tagClusters(ce, sample = "zf_unfertilized_egg"))

## ----CumulativeDistribution, echo=FALSE, fig.align="center", fig.cap="Cumulative distribution of CAGE signal and definition of interquantile width"----
knitr::include_graphics("images/CumulativeDistributionAndQuantiles.svg")

## -----------------------------------------------------------------------------
cumulativeCTSSdistribution(ce, clusters = "tagClusters", useMulticore = T)
quantilePositions(ce, clusters = "tagClusters", qLow = 0.1, qUp = 0.9)

## -----------------------------------------------------------------------------
tagClustersGR( ce, "zf_unfertilized_egg"
             , returnInterquantileWidth = TRUE,  qLow = 0.1, qUp = 0.9)

## -----------------------------------------------------------------------------
plotInterquantileWidth(ce, clusters = "tagClusters", tpmThreshold = 3, qLow = 0.1, qUp = 0.9)

## -----------------------------------------------------------------------------
aggregateTagClusters(ce, tpmThreshold = 5, qLow = 0.1, qUp = 0.9, maxDist = 100)
ce$outOfClusters / ce$librarySizes *100

## -----------------------------------------------------------------------------
consensusClustersGR(ce)

## -----------------------------------------------------------------------------
annotateConsensusClusters(ce, exampleZv9_annot)
cumulativeCTSSdistribution(ce, clusters = "consensusClusters", useMulticore = TRUE)
quantilePositions(ce, clusters = "consensusClusters", qLow = 0.1, qUp = 0.9, useMulticore = TRUE)

## -----------------------------------------------------------------------------
consensusClustersGR( ce, sample = "zf_unfertilized_egg"
		               , returnInterquantileWidth = TRUE,  qLow = 0.1, qUp = 0.9)

## -----------------------------------------------------------------------------
exportCTSStoBedGraph(ce, values = "normalized", format = "bedGraph", oneFile = TRUE)

## -----------------------------------------------------------------------------
exportCTSStoBedGraph(ce, values = "normalized", format = "BigWig")

## ----CTSSbedGraph, echo=FALSE, fig.cap="CAGE data bedGraph track visualized in the UCSC Genome Browser"----
knitr::include_graphics("images/CTSSbedGraph.svg")

## -----------------------------------------------------------------------------
exportToBed(object = ce, what = "tagClusters", qLow = 0.1, qUp = 0.9, oneFile = TRUE)

## ----tagClustersBed, echo=FALSE, fig.align="center", fig.cap="Tag clusters visualization in the genome browser"----
knitr::include_graphics("images/TagClustersBed.svg")

## -----------------------------------------------------------------------------
# Not supported yet for CAGEexp objects, sorry.
# getExpressionProfiles(ce, what = "consensusClusters", tpmThreshold = 10, 
# 		nrPassThreshold = 1, method = "som", xDim = 4, yDim = 2)

## -----------------------------------------------------------------------------
# Not supported yet for CAGEexp objects, sorry.
# plotExpressionProfiles(ce, what = "consensusClusters")

## ----ConsensusClustersExpressionProfiles, echo=FALSE, fig.align="center", fig.cap="Expression clusters"----
knitr::include_graphics("images/ConsensusClustersExpressionProfiles.svg")

## -----------------------------------------------------------------------------
# Not supported yet for CAGEexp objects, sorry.
# class3_1 <- extractExpressionClass(ce, 
# 		what = "consensusClusters", which = "3_1")
# head(class3_1)

## -----------------------------------------------------------------------------
# Not supported yet for CAGEexp objects, sorry.
# exportToBed(ce, what = "consensusClusters", 
# 		colorByExpressionProfile = TRUE)

## ----ConsensusClustersBed, echo=FALSE, fig.align="center", fig.cap="Consensus clusters colored by expression profile in the genome browser"----
knitr::include_graphics("images/ConsensusClustersBed.svg")

## -----------------------------------------------------------------------------
ce$group <- factor(c("a", "a", "b", "b"))
dds <- consensusClustersDESeq2(ce, ~group)

## -----------------------------------------------------------------------------
cumulativeCTSSdistribution(ce, clusters = "consensusClusters")

## -----------------------------------------------------------------------------
# Not supported yet for CAGEexp objects, sorry.
# scoreShift(ce, groupX = "zf_unfertilized_egg", groupY = "zf_prim6",
# 		testKS = TRUE, useTpmKS = FALSE)

## ----ShiftingScore, echo=FALSE, fig.cap="Calculation of shifting score"-------
knitr::include_graphics("images/ShiftingScore.svg")

## -----------------------------------------------------------------------------
# Not supported yet for CAGEexp objects, sorry.
# shifting.promoters <- getShiftingPromoters(ce, 
# 		tpmThreshold = 5, scoreThreshold = 0.6,
# 		fdrThreshold = 0.01)
# head(shifting.promoters)

## ----ShiftingPromoter, echo=FALSE, fig.cap="Example of shifting promoter"-----
knitr::include_graphics("images/ShiftingPromoter.svg")

## ----create_df----------------------------------------------------------------
TSS.df <- read.table(system.file( "extdata/Zf.unfertilized.egg.chr17.ctss"
                                , package = "CAGEr"))
# make sure the column names are as required
colnames(TSS.df) <- c("chr", "pos", "strand", "zf_unfertilized_egg")
# make sure the column classes are as required
TSS.df$chr <- as.character(TSS.df$chr)
TSS.df$pos <- as.integer(TSS.df$pos)
TSS.df$strand <- as.character(TSS.df$strand)
TSS.df$zf_unfertilized_egg <- as.integer(TSS.df$zf_unfertilized_egg)
head(TSS.df)

## ----coerce_to_df-------------------------------------------------------------
ce.coerced <- as(TSS.df, "CAGEexp")
ce.coerced

## ----sessionInfo--------------------------------------------------------------
sessionInfo()